Qualitative features of a nonlinear, nonlocal, agent-based PDE model with applications to homelessness

In this paper, we develop a continuum model for the movement of agents on a lattice, taking into account location desirability, local and far-range migration, and localized entry and exit rates. Specifically, our motivation is to qualitatively describe the homeless population in Los Angeles. The model takes the form of a fully nonlinear, nonlocal, non-degenerate parabolic partial differential equation. We derive the model and prove useful properties of smooth solutions, including uniqueness and [Formula: see text]-stability under certain hypotheses. We also illustrate numerical solutions to the model and find that a simple model can be qualitatively similar in behavior to observed homeless encampments.


Introduction
Homelessness is a growing problem for society, with large cities like Los Angeles having homeless populations of nearly 40,000 people 32 at present. To date, correlations have been established between homeless populations in different cities; factors such as the cost of rent, number of affordable housing units, and poverty rates 15,19 are relevant. But outside of a few anecdotal examples in the social science litera-ture, 49 little is known about the daily lives of homeless people, how they became homeless, and how long they may remain homeless. There is also the question of how the homeless population changes over time and space, which is essential for planning social support services to counteract this problem. A few studies address the dynamics of the homeless population, 19 but there are still many unanswered questions. Within Los Angeles, the homeless population is a big concern. Complaints from Los Angeles residents regarding homeless encampments (see Fig. 1) are among the most common reason residents call the city. 42 The encampments can form quickly and require some careful consideration in terms of clean up. 41 But the encampments may also serve as places where homeless people themselves live preferentially, feeling safer. 22 Human interactions with their environment can be studied both from a social science perspective and through mathematical and statistical methodologies. Models of human social interaction have studied in many contexts, for example voting, 27 lobbying; 18 leader-follower dynamics; 1 and evolutions of opinions, 20,50 including through agent-based models of continuous-valued opinions within a finite group of peers. 36 Models relating humans and the environment have also been put forth, such as those pertaining to climate change. 44 Pedestrian travel has been studied experimentally and modeled with mathematics, 38 even taking into account personal objectives of minimal travel time and avoiding high density areas; 24 and statistics can be employed to describe large crowds at special events, including the fatal pressures that can build up in the center, asphyxiating standing people. 33 Mathematical models for riots have also been developed. 3 On a larger scale, agent-based models have been used to study residential burglary, taking into account historydependent target attractiveness, 43,45 truncated travel distances, 39 and agents burglarizing according to independent Poisson clocks. 48 Mathematical models of traffic Qualitative features of a nonlinear, nonlocal, agent-based PDE model 3 flow have also been developed. 16 More broadly, urban crime has received a lot of attention in models that incorporate police into the models; 13,21,25,40,47,51 models that predict locations of offenders; 37 rigorous proofs of solution properties of such equations; 5,35 and the study of emergent patterns and behavior. 4,14,26 Species aggregation has also been thoroughly studied, 7 with many authors modeling such phenomena by nonlocal differential equations. 2,10,46 On even larger scales, we can see the natural emergence of pack formation in predator-prey systems, for sufficiently strong intraspecies competition; 6 descriptions of human territorial conquest; 17 and ecology models for a species in an environment with patchy resources. 34 But to our knowledge, little has been done directly to model homelessness, which is a growing societal concern.
In our study, we address the dynamics of the homeless population on small geographic scales. We develop an agent-based model for the homeless population living on the streets, which, in the continuum limit, yields a partial differential equation for the population density. Our data sources for the homeless population are a collection of annual point-in-time estimates broken down by census tract from the Los Angeles Homeless Services Authority. [28][29][30][31] We remark that accurate counting of the homeless population is challenging at best, with some authors suggesting the counts could be significant underestimates. 8 Our work focuses upon building a qualitative model so absolute precision is not necessary. This paper makes three main contributions to the scientific community: first, the derivation of a new continuum partial differential equation to describe homeless populations; second, rigorous analysis of basic properties of solutions; and finally, in using the model to yield qualitatively consistent behaviors observed in the true homeless population. The paper structure is as follows: in Sec. 2, we derive the model; in Sec. 3, we prove useful properties of solutions; in Sec. 4, we demonstrate how the model can qualitatively describe encampment formation; and finally, in Sec. 5, we summarize our work and discuss future research.

Model Formulation
In describing the homeless population, we consider several processes taking place at each location (see Fig. 2): • Entry into population -due to local features of the area such as cost of living, etc. individuals may become homeless; • Exit from population -through localized social services, family support, death, or other processes, individuals may cease to be homeless; • Diffusive movement -a homeless individual who does not like their current area due to lack of resources or other factors may venture to a surrounding neighborhood; and • Nonlocal movement -they may intentionally travel to another part of the city, such as through public transit. The features of the cityscape vary by location so some areas are more residential, some are more industrial, etc.

Derivation
We use an Eulerian frame of reference discretized by a regular lattice of points I in R d , with spacing of δx in each coordinate direction and volume δx d . We also discretize time into time steps of δt with t k = kδt, k = 0, 1, 2, . . . . For modeling homeless populations, our main focus is d = 2, but d = 1 could describe populations of people (or even other species) in a long/narrow geometry and d = 3 could describe populations that can freely move in three dimensions such as birds or fish. We concern ourselves with the evolution of N k i , the number of individuals occupying site i ∈ I at time t k . Intrinsic to each site i at time t k , we denote an "attractiveness" A k i ∈ (0, 1). This could be influenced by the resources available at the location. In general, we denote subscripts for the spatial location index and superscripts for the time index. Over each very small time interval δt, we assume the following: • individuals who do not leave the population will remain at their current location with probability A k i ; • individuals who do not leave the population and choose to leave their current lattice point will make deliberate travel, such as with public transit in the case of homelessness, at a Poisson rate Γ k i ; and • if an individual remains in the population, does not stay at their location, and does not travel deliberately from their site to another, they will travel to each neighboring site with equal probability.
We incorporate all terms, allowing them to potentially be "in balance" so that transportation and movement may be comparable to arrivals/exits. In reality some of these effects may be negligible in describing the actual homeless population but we opt for generality.
At each site i, up to O(δt) the following probabilities describe an individual: Pr(stay at i and stay in Pr(leave i but not travel far and stay in I) (2.8) so that on average at site i starting at time t k , over each δt: • E k i δt people enter the population; δt people travel deliberately to another point; and We make one further assumption that there exists a transition matrix T k = (T k ji ) (j,i)∈I 2 describing the probability an agent deliberately travels from site j to site i. We denote j ∼ i to signify that i = j are neighbors and j → i to signify that i = j and there is deliberate (potentially far-range) movement from j to i. We also 6 M. R. Lindstrom & A. L. Bertozzi assume that each lattice point has n = 2d neighbors. Then The terms of the form j∼i • k j − n• k i are the discrete second-order centered difference Laplacian at position i and time t k multiplied by δx 2 . Denote the discrete Laplacian by ∆. Dividing the equation by δt, we have Under a diffusive scaling, such that δx 2 nδt = D = O(1), in taking the limit, we furnish the partial differential equation for a spatial domain Ω where ρ is a spatial density of agents; η is an entry rate per unit area; ω is an exit rate; τ governs transitions from y to x at time t; γ is a travel rate term; and a is a continuous attractiveness field. We can express (2.9) more nicely by defining the unattractiveness u = 1 − a (2.10) Qualitative features of a nonlinear, nonlocal, agent-based PDE model 7 so that over a domain Ω subject to the normalization of τ ≥ 0 with Ω τ (y, x, t)dx = 1. (2.12) We remark that a hyperbolic scaling with δx/δt = O(1) would not allow for diffusion in the system and we would lose out on the local travel term altogether, only having deliberate travel terms. This is equivalent to setting D = 0 in (2.11). We do remark, though, that sometimes hyperbolic scalings are necessary when the reality of finite propagation speeds cannot be well approximated through a diffusive limit. 11 For this preliminary work, we focus on a diffusive approximation.

Nondimensionalization
We begin with the PDE where the bars are scales and the asterisk variables are dimensionless. We also define Ω * = 1 x d Ω. Noting that if y =xy * then dy =x d dy * , we have We can adopt scalings relevant homeless count data: given the homeless counts are done annually, we chooset to be 1 year so an O(1) timescale roughly represents an interval between counts. We also choosex as a characteristic length within our region of interest (later we use the geometric mean of the length/width of the spatial domain), allowing the entire domain to have approximately unit length. There is also a natural scaling forρ as the initial average population density over the region. If the dimensionless density is ≈ 1 then it is roughly average. These choices suggest the derived scales (to make more dimensionless ratios equal to 1) Then for a constant δ = Dt/x 2 , and after removing the asterisks, we have In general, the unattractiveness u may depend upon ρ and other localized features. In Sec. 2.3.2, we choose a particular form for u (Eq. (2.16)). Because of this potential density-dependence, the system of equations is nonlinear. We can interpret (2.13)-(2.15) as consisting of a nonlinear diffusion operator where ρu is diffusing; a local source term η; a localized exit rate ω; a nonlocal operator I[γρu](·) describing intentional travel, and a local loss due to intentional travel γρu. Note that the nonlocal operator is not a convolution.

Notation
To denote solution spaces, we may explicitly label a function's argument and specify the continuity/differentiability assumed with a subscript with that label. We denote to be a function f that is twice continuously differentiable in x and once in t. We use similar notation for continuity in higher derivatives. The notation C represents continuity in all arguments, possibly with superscripts to denote the number of derivatives, e.g. C ∞ .

Assumptions
We denote our spatial domain to be Ω and assume it is bounded and simply connected. For simplicity in the analysis, we will choose Ω to be the torus Periodic boundary conditions actually exist in cities e.g. with perimeter roads, however here the choice is more for convenience of the mathematics. Also cities like Los Angeles can have repeating patches of residential (both affluent and disadvantaged) neighborhoods and commercial regions, which can give the impression of a repeating pattern. For numerical studies, we sometimes use no-flux boundaries in R 2 , identifying a flux of −δ∇(ρu).
We denote the following: • ρ(x, t) : Ω × R ≥0 → R ≥0 , the population density (people per unit area); • θ(x, t) : Ω × R ≥0 → R p , a features vector of size p, which can vary over space; Qualitative features of a nonlinear, nonlocal, agent-based PDE model 9 Fig. 3. Plots of unattractiveness versus density withρ = 1. The unattractiveness increases with density for all fixed κ. It is always at least u + − u − and never exceeds u + . As κ increases, the unattractiveness decreases.
the (dimensionless) unattractiveness for 0 < u − < u + < 1, 0 <ρ (see Fig. 3 for qualitative depiction); • τ (y, x, t) : Ω × Ω × R ≥0 → R ≥0 , the travel term (per unit area probability), such that , the entry rate (people per unit area per unit time); For simplicity later on, we will often write κ or κ(x, t) instead of κ(θ(x, t)), etc. The fact that unattractiveness is density-dependent is a hypothesis based on the observation that no region can supply unlimited resources. The potential for including more general τ to include preferential travel to locations with a higher population density is briefly discussed in Sec. 3.3 as this could be a possible mechanism. 22 For simpler analysis, however, we do not consider such τ .

Model exploration
To build a basic understanding of the model features, we study the solutions numerically and vary the parameters. We do this in one dimension on the torus with length 1. The numerical method employed is explained in Appendix A. We wish to qualitatively understand the effects of the different model features. As part of this work, we locally perturb some features sometimes making use of the bump function with a wide range of x -values where it is nearly constant (owing to the large power of 20). We remark that Υ has a support of length 2, which is larger than the unit torus, but by shifting and rescaling its argument, we ensure our use of Υ does not violate the properties of the torus. We perform a series of numerical experiments upon the model, as seen in Fig. 4 with explicit parameter listings in Table 1. In particular, we study the following: • baseline: We study how the population density evolves from an initial distribution when all the parameter functions are constants, observing a steady approach of the population density to a spatially constant solution at steady-state. • enhanced local entry: From the baseline, we increase the entry rate η in a region and at steady-state find the population is largest where there is more entry. • enhanced local exit: From the baseline, we increase the exit rate ω in a region and at steady-state find the population is reduced where there is a greater exit rate. • enhanced local to far migration: From the baseline, we increase the travel rate γ in a region and at steady-state find the population has diminished due to a higher rate of moving away from the area. • biased transfer: From the baseline, we choose τ to be biased to relocate individuals to one particular region; at steady-state, the population density is higher at this destination. • exponential decay transfer: From the baseline, we choose τ such that from y, the probability density in moving to x decays exponentially with |x − y|. Compared to the baseline, the change is quite small, but we observe the population decreases slightly slower in more concentrated regions since the people are not travelling as far.
• textured relative capacity: From the baseline, we vary the relative capacity κ in an oscillatory fashion; we find that at steady-state, the population density fluctuates with this varying relative capacity. • enhanced diffusion: From the baseline, we raise the diffusion coefficient and find the behavior can quickly be dominated by diffusion, resulting in flat density profiles. • more density sensitivity: From the baseline, we make the unattractiveness more sensitive to density variations by decreasingρ in (2.16). We observe the population disperses more rapidly in regions where the population density was initially highest. • narrow unattractiveness range: From the baseline, we modify u − and u + in (2.16) so that the unattractiveness is almost constant, regardless of ρ or κ. We also perturb κ as in the "textured relative capacity" experiment. Here we find the steady state population is much less influenced by the relative capacity. • region becomes uninhabitable: From the baseline, a region Q becomes uninhabitable. This is modeled by reducing the relative capacity to zero inside Q and modifying so that transfer into the region is zero and out of the region is sizeable. In addition the travel rate is increased inside of Q and the entry rate is zero within Q. We find the population chooses to relocate to more habitable areas.

Properties of Smooth Periodic Solutions
We study the solutions to

2)
Qualitative features of a nonlinear, nonlocal, agent-based PDE model 13 Here, we work in a spatial domain Ω = T d , the d-dimensional torus. For T > 0, we denote )} with t-derivatives understood to be right-, respectively, left-derivatives at t = 0 and t = T . In general, we will be assuming the existence of solutions to (3.1)-(3.4) within the space U T . If we write U ∞ we refer to a solution that exists globally for all time. At all times, we assume that hypotheses 1 and 2 are satisfied.
As Eqs. (3.1)-(3.4) are parabolic and nondegenerate (u is never zero), we anticipate global existence of smooth solutions, without the formation of shocks or phenomena such as finite-time blowup. However, such proofs are beyond this paper. A reader interested in the question of existence of solutions could refer to proofs of local existence to similar or related models such as for chemotaxis 9,12,23 or residential burglary. 43 Remark 3.1. The solutions we consider are continuous on the torus T d , which is compact, so we can use sup and max, respectively, inf and min, interchangeably.

Useful properties of ρu
The term ρu appears many times in our analysis and we make a list of some useful properties. Note that ρu = u + ρ − κM (ρ) where we define (3.7) Thus, for ρ ≥ 0, We also have that (3.14) In one spatial dimension, we also have Owing to (3.11), hypothesis 2, and 0 ≤ κ ≤ 1, we have where the subscripts in the Lipschitz spaces denote the bounding constants.

Results
We have a series of results of solution properties below. All proofs are provided in the paper, but the following lemma and subsequent three propositions are proved in Appendix A.2.
with F monotonically nondecreasing with respect to q, and Θ[q](x, t) is an operator depending on q, x ∈ T d , and t ∈ R ≥0 such that if q achieves a global maximum over T d at x then Θ[q](x, t) ≤ 0, respectively, if q achieves a global minimum over    Proof. From subtracting the respective PDEs, we have that By multiplying the equation by (ρ 1 − ρ 2 ), integrating over T d , and integrating by parts once, we have d dt where we used the fact that (ρ 1 − ρ 2 ) and γ(ρ 1 u(κ, ρ 1 ) − ρ 2 u(κ, ρ 2 )) will have the same sign (because ρu is positive and monotonically increasing in ρ). Now, by (3.16), Working with the derivative terms, we have 20) for some K that depends on max t∈[0,T ] ∇ρ 2 (·, t) L ∞ (T d ) . It is finite as solutions are twice continuously differentiable in T d so the gradient cannot blow up. The last inequality stems from the final (and positive) term in the integrand being Lipschitz: since ∇ρ 2 is bounded on T d (ρ 2 is C 2 in x) and using (3.17)-(3.18), the term being squared in the numerator is Lipschitz. Also, the denominator is bounded below by 4(u + − u − ). Whence, by combining inequalities (3.19) and (3.20) and by a standard application of Grönwall's inequality, the result is proven.
Qualitative features of a nonlinear, nonlocal, agent-based PDE model 17 Proof. The proof is trivial by using the arguments that prove 3.1: replace ρ 1 with and ρ 2 by ρ.
Remark 3.3. We are assuming global existence here in order for the limit to make sense.
Proof of theorem 3.2. Our argument will come from two parts: first, we will establish a uniform upper bound for the solution over T d that converges toρ as t ↑ ∞. We will then show that the solution cannot remain smaller thanρ except on sets whose measure vanishes. With κ constant in space, ρu = F (ρ, t) is monotonically increasing with ρ for each fixed t. Then (3.1)-(3.4) become We wish to change variables to remove as many terms from (3.21) as possible so as to apply the maximum principle supplied by Lemma 3.1. For smooth ξ 1 and ξ 2 , write ρ = e ξ1(t) q(x, t)+ξ 2 (t) so thatξ The dot signifies a time derivative. We choosė so that in simpler terms Note that γe −ξ1 (F − F ) is an operator that satisfies the conditions of Θ in Lemma 3.1. By Lemma 3.1, the global maximum for q in U T is achieved at t = 0. By choice of initial conditions (3.22) 2 and (3.23) 2 , ρ(·, 0) = q(·, 0) and thus By recalling η, ω are constant in space and solving the ODEs (3.22)-(3.23), we can precisely state that This proves directly that the set of points upon which ρ exceedsρ + must have vanishing measure as t ↑ ∞ for all > 0. Having established an upper bound, we now work on the second part of the proof. Integrating (3.21) over T d we have (as in proof of Proposition 3.2 in Appendix A.2) that Let > 0 and define S (t) = {x ∈ T d |ρ(x, t) ≤ρ(t) − }. We show that in measure, ρ →ρ(t).
Observe that by Eqs. (3.25) and (3.24) Rearranging the first line and the last line, using |T d | = 1, and recalling there is an inequality produces Since we assume here that ω is strictly bounded below by 0 so that ∞ 0 ω(s)ds diverges, the right-hand side decays to zero and thus In one dimension, we can prove a stronger result.
Remark 3.4. Note that we require three continuous derivatives in space. Also, we potentially allow ω = 0 here.

Qualitative features of a nonlinear, nonlocal, agent-based PDE model 19
Proof of theorem 3.3. We will prove that the gradient ρ x tends to zero uniformly. We can rewrite (2.13) in one spatial dimension and with the imposed hypotheses as We begin by taking an x−partial-derivative using (3.15) with κ, ω, and γ constant in space to obtain since the nonlocal operator results in a constant in space. By interchanging x-and t-derivatives since the mixed partial derivatives are assumed continuous (and hence equal), we rewrite the equation with ψ = ρ x Now, we multiply by ψ so that Let q = 1 2 ψ 2 . If we can prove q → 0 uniformly over T we are done. Note that, we have Note that q ≥ 0 and if q(x, t) = r(x, t)e ξt then r t − δ(u + − κM (ρ))r xx + δ(u + − κM (ρ))ψ 2 x re −ξt + 6κM (ρ)ψ x r − κM (ρ)ψe −ξt ψ 2 x = −(ξ + 2(ω + γ(u + − κM (ρ))))r. If the supremum of r is zero, there is nothing more to prove. Otherwise, at any nonzero local maximum for r, ψψ x = 0 and since ψ = 0 in such a case, we have ψ x = 0. If ξ is chosen so that ξ + 2(ω + u + − κM (ρ)) ≥ 0 then by Lemma 3.1, the global maximum for r is achieved at t = 0. Many ξ can be chosen but to establish decay, we choose ξ = −2m if ω ≥ m > 0 with γ reaching zero and otherwise we pick ξ = −2m(u + − u − ). This establishes 1 2 ψ 2 ≤ e −mt max T 1 2 ψ(·, 0) 2 so the gradient of ρ tends to zero uniformly. Combining this fact with Theorem 3.2 proves the result.
We examine the consequences of Theorem 3.3 in Fig. 5 by plotting the decay of the sup-norm of the gradient of a solution subject to time-dependent but spatially constant forcing. We now turn to the problem of parameters and solutions that are timeindependent and their perturbations. alongside (2.14) and (2.15). Then, for sufficiently small (in · L ∞ (T d ) ), periodic ρ 0 ∈ L 2 (T d ), the solution to (2.13) with ρ(·, 0) = ρ 0 +ρ 0 upholds ρ → ρ 0 in L 2 (T d ) as t → ∞.

Generalizations
Some of the proofs provided here could readily be generalized to include other features. For example, the proof of positivity, and the bounds on the L 1 and L ∞ norms of solutions depended upon the positivity, normalization, and boundedness of τ . But more complicated τ choices could also be used, such as those of the form τ (y, x, t, ρ/ ρ L 1 (Ω) ) to furnish the same results, allowing effects of aggregation due to localized peaks in the population density. Other proofs may require more control over τ and further properties of solutions ρ.

Qualitative Features of Real Data
After nondimensionalization and smoothing, the homeless population densities for four consecutive years are plotted in Fig. 6. Within the real data, we remark that while some areas have high homeless population densities consistently from year to year, encampments can form over the course of a year and a new bump in the density appears.
The PDE model has many parameters, most of which can vary over space and time. As such, there is a danger of overfitting. In order to illustrate qualitatively Fig. 6. Plots of the homeless population densities for a subset of Los Angeles. The population has been nondimensionalized by a scale of 68 mi −2 and the length scale is 12 mi. The data have been smoothed: population densities were interpolated onto a regular mesh from LAHSA data and then locally averaged over a radius of 1 mi. Here, x and y denote spatial position and ρ is population density.
consistent behavior in our model with the real data of encampment formation, we restrict ourselves to describing the formation of a new encampment (a new local maximum). To achieve this, we consider a scenario where all of η, ω, and γ are constant in space and time, the travel term τ is given by τ (y, x, t) = κ(x, t) Ω κ(x , t)dx where κ steadily increases in a region over a dimensionless time interval of l. Over the domain we assume no-flux boundary conditions. This results in the numerical solutions being consistent with real encampment formation (see Fig. 7).

Conclusions and Future Work
We formulated a continuum PDE model to describe the evolution of the homeless population density. We have proven that smooth periodic solutions to the model equations enjoy a maximum principle, a positive population density, uniqueness, a flattening phenomena with spatially uniform forcing, and L 2 -stability for constant transfer rates and symmetric travel kernels. The model is well behaved and, for suit-able choices of parameters, it can produce population changes that are qualitatively consistent with real homeless population data. From this preliminary model, many new avenues of research open up. It would be worth understanding how to empirically model the various parameters/functions in (3.1)-(3.4) so that the model can be used in a quantitative capacity and in identifying what can be done to combat homelessness. From the viewpoint of mathematical analysis, a rigorous proof of existence of solutions is worth pursuing. It would also be interesting to study this model or its natural variants as applied to other contexts. = δ∆( (3) u( (2) )) for (3) , also with semi-implicit boundary conditions; (4) solve = δ∆(ρ k+1 u( (4) )) for ρ k+1 , also with semi-implicit boundary conditions.
The Laplacian term solves the equation with the appropriate boundary conditions so, we split it into three. This ensures the nonlocal step and stepping forward with the sink will be done with the correct boundary conditions and so that the output also has the proper boundary conditions. The nonlocal operator is done explicitly to avoid inverting a large, dense matrix in the case of an implicit scheme. It is combined with the positive source term η. Positivity can be preserved when done explicitly for those terms. The sink term −ωρ is dealt with implicitly to preserve positivity.
The quadrature for the nonlocal movement is delicate. In second-order quadrature, values of γρu in the integral are weighted by either 1/2 or 1/4 at the boundary. Proper mass balance is achieved by weighting the corresponding sinks of −γρu by either 1/2 or 1/4.