Beyond Normal: Learning Spatial Density Models of Node Mobility

There is an accumulating body of work of using machine learning (ML) in network performance models [1, 2, 3, 4]. ML models are considered to be a natural choice to master the complexity of expressing functions that quantitatively capture the multiple interactions existing in the performance exhibited by a networked system consisting of many components. One of the major advantages of ML models, like deep neural networks (DNN), is their ability to generalize out of a set of (training set) examples [5], which, alongside their ability to relate data to underlying latent-space representations, has contributed to their wide success and acceptability. So much so, that a trained DNN is often invoked as a subroutine in bigger architectures, e.g., in the context of reinforcement learning (RL), to compute complicated functions—an area in which DNNs excel.

We focus on a class of performance evaluation problems that demand the derivation of steady-state distributions over geo-coordinates. Such problems arise in mobility modeling using continuum approximation of mobile nodes, seen as points on a disk on a 2-dimensional terrain, and influenced by particular movement-related constraints. The disk area in the systems under study imposes spatial constraints (i.e., a circular boundary) and symmetries (i.e., rotations around the center are, often, inconsequential). Related mobility modeling results have been reported in various previous articles [6, 7, 8]. Some of the models, like the one we use as an example in this paper, capture the “attraction” of mobile nodes at certain locations in space, e.g., because of the presence of chargers etc. Therefore, our approach can be applied to the broader literature of mobility models where the available space is circularly confined but, within it, there are “bumps” in the density owing to various reasons. Note that such density models can be close approximations of non-uniform mobility in cities, where the extent of the city is circumscribed. In summary, we are looking at density models that capture non-trivial shapes of spatial distributions over a disk.

Our methodology employs Mixture Density Networks (MDNs) [9], starting by considering typical Gaussian mixture models [10], due to their universality at approximating densities [11] and because of their extensive use in many practical settings. Yet, once again, it becomes evident that off-the-shelf ML strategies, are not as effective as fine-tuned ones. While the statement might appear obvious, it scratches the surface of what is the additional degree of fine-tuning necessary to produce good results that balance modeling effort for performance (seen as accurate reconstruction of distributions in the MSE or KL-divergence sense). In terms of modeling effort, we will demonstrate that abandoning Gaussian mixtures can be advantageous. At the same time, we need to consider that the training data are few because they are also the result of long computation. This computation involves either simulations or numerical approximations involving complicated formulas.

Existing mobility models range from simple random waypoint (RWP) to fairly elaborate models where nodes are attracted to specific regions/points because there exist spatial features. Closed-form analytical descriptions of the density rarely exist but only for a few cases, one such being the noted example of the RWP density [12]. An example of non-trivial density which lacks a closed-form solution and is computationally demanding even for its numerical approximation, is the one derived from the work in [8] and shown in the example of Figure 1. Zooming into the features of the distribution reveals subtleties like the density “peak” drops in a non-symmetric manner to the left and to the right of the peak, while at some distance around this peak there exists a “dip” in the density. The notable background “dome” is the RWP residual component of the distribution.

For an ML model to be convincing, it is necessary that it is also interpretable. The examples shown here illustrate what is a criterion for such interpretability. For example, a model that captures the RWP dome in one component distribution, and the off-center peak in another component can be claimed to be more interpretable than one that does not provide such separation.

The contributions of the current paper are all related to learning spatial distributions by means of mixtures of Möbius distributions and have the characteristics of:

• •

  capturing distributions defined on the disk, as required by many existing mobility models,

• •

  capturing distributions that include “features” such as local attraction points/peaks,

• •

  requiring a small number of training set data to produce a good model, and

• •

  resulting in component distributions that are more interpretable than Gaussians.

In Section 2 we provide a specific problem definition from the point of view of approximating a particular variety of density functions derived from mobility. The point is that a wider set of similar distributions can be handled using similar models. Section 3 introduces the family of Möbius distributions, in a mixture setting analogous to mixture-of-Gaussians, providing motivating examples of how (better) Möbius performs compared to Gaussians. The summary of learning results is given in Section 4 where we also comment on the fragility and worst-case outcomes of the Möbius mixtures compared to Gaussians. Section 5 links our work with a number of current trends in learning models comparing and contrasting them. Finally, Section 6 summarizes the paper and points to generalizations like non-circular boundaries, more general distributions, etc.

When engineering systems, such as communication systems, it is convenient to have closed-form formulations for relevant probability density functions and/or their related statistical properties. For example, the closed-form relation between the average delay and arrival rates at a queueing system, as per Little’s law, is often employed in solving bigger network design problems [13]. Here we are interested in aiding the solutions of spatial optimization problems, i.e., placement of resources, where the placement impacts a probability density function. Such optimization problems arise e.g., in the study of node density of mobile nodes moving around under certain mobility rules. Part of the assumptions when studying such spatial densities is how far is the extent of the space over which the distribution should be described. Motivated by the fact that many applications involve the placement of resources in cities and other population centers, an often applied assumption is that the area is circumscribed by a circle. Hence, densities over a disk are commonplace [14].

A well-studied baseline mobility model is the RWP model, which results in a closed-form probability density function of nodes across space [12]. RWP assumes that nodes continuously move around along straight paths from waypoint to next waypoint. Yet, when RWP is extended to capture more subtle movement rules, e.g., movement predicated on the state of the mobile nodes, and non-uniformity of space, extracting a closed-form solution is challenging. For example, RWP alone does not have a concept of points of attraction. Attraction points could be inflicting a change in the mobility behavior depending on the state and proximity of the mobile nodes. An example of this is the case of placing charging stations at specific locations and, as a result, attracting users to go to them to charge their mobile devices. One expects the density of mobile nodes to be higher around chargers. The less the overall available energy, the more the likelihood they will be attracted to the charger and the density will “bulge” around the charger which will influence their trip between successive waypoints. Deriving a function to describe the density based on pertinent parameters, e.g., location of the charger, average energy capacity of nodes, etc., is necessary in such an application.

In previous work [8], we have derived the node density distribution for a network of energy-constrained mobile nodes. We are repeating here some of the key formulas behind the analytical expression of the node density to illustrate their intertwined nature. Note that we have thus far found no closed-form solution for the distribution. Additionally, it is computationally expensive to produce a numerical approximation to the solution, in particular due to the multiple integrals involved. The purpose of the learning described in the next section is to learn an efficient function that evaluates to this complicated density using a small training set since each training point requires substantial numerical processing or lengthy simulations.

The density is a complicated one, owing to the fact that it models chargers that “attract” nodes whose energy is depleted. The trajectory of the mobile nodes is altered, from continuing to their destinations, to diverting to a charger before proceeding to their destinations. The decision to divert is only taken if the mobile node is a certain, relatively small, distance from the charger—equivalently seen as the “attraction” range of a charger. More formally, suppose that $c=(r,\theta)$ (our first two parameters) denotes the charger’s location in polar coordinates. That is, $r$ and $\theta$ are radial and angular coordinates respectively, where the radius of the disk (centered at (0,0)) is normalized to a unit (hence $r<1$). The (energy-depleted) mobile node will divert for recharging, provided that its current location, denoted by $z$, is within the (closed) disk, $\overline{\mathcal{D}}(c,d^{\prime})=\left\{z\mid|z-c|\leq d^{\prime}\right\}$, i.e., if it is within $d^{\prime}>0$ distance from the charger ($d^{\prime}$ being the third parameter of our model). The limitation of attracting a depleted node only if it is within an attraction range models the reluctance of a node to divert significantly from their path to the destination for the sake of recharging (which they can do at their destination point anyway).

Besides adequate proximity to the charger, the other premise of the diversion for recharging, i.e., that the node is depleted, is made tantamount to the event that the node has traveled a distance $d>0$ ($d$ is our fourth parameter) after departure from its last waypoint. This projection of energy consumption processes onto distance is a modeling convenience allowing us to lump together the effect of many factors (e.g., rated power, moving speed) into a typical “energy budget” expressed as a distance threshold $d$, in the same units as those for determining the attraction ($d^{\prime}$) to the charger. An implicit assumption is that the energy for locomotion is separate from the energy for computation and communication—a scenario which fits the case of human movement with portable devices (e.g., smartphones, wearables) being carried around.

Even in the presence of a single charger, the formula we derived for the probability density is non–trivial and has so far evaded a closed form expression [8], leading to needs for numerical evaluation. Note that we have been able to derive closed form expressions for the 1-dimensional case [15], but its components and constraints are significantly less challenging than the 2-dimensional case. The following formula of node density for the 2-dimensional case applies when $r<1-(d+d^{\prime})$, i.e., when the charger is placed relatively close to the center to allow nodes that start moving from the boundary to have a chance to be attracted to the charger when their energy is depleted. The high-level structure of the spatial density $f_{Z(t)}(z)$ satisfies the following equation:

where $\alpha$ denotes the expected distance (including diversions to the charger) between any two waypoints, $f_{Z(t)}(z)$ is the (stationary) probability density of the node’s location $Z(t)$ over time $t$, and $f_{\mathrm{RWP}}(z)$ is the density assuming only RWP mobility (i.e., with no diversion to the charger).

Note that the movement area and thus the support of $f_{Z(t)}(z)$ is a unit disk $\Omega=\{z\in\mathbb{R}^{2}\mid\left|z\right|\leq 1\}$. The density within the attraction range of the charger (i.e., for $z\in\overline{\mathcal{D}}(c,d^{\prime})$) can be captured, excluding the basic RWP component $f_{\mathrm{RWP}}(z)$ (and the coefficient $(\pi^{2}\alpha)^{-1}$), by an infinity at the charger (i.e., for $z\in\{c\}$) and a subtraction of two non-trivial functions $g_{1}(z)$ and $h_{1}(z)$; a similar superimposition onto $f_{\mathrm{RWP}}(z)$, with the subtraction $g_{0}(z)-h_{0}(z)$, is also found for formulating the nodal density beyond the attraction range (for $z\in\Omega\setminus\overline{\mathcal{D}}(c,d^{\prime})$). The respective forms of $g_{1}(z)$, $h_{1}(z)$, $g_{0}(z)$, and $h_{0}(z)$ are given by

where

and where $l_{2}$, $l_{3}$, $l_{4}$, $l_{12}$ ($\hat{l}_{12}$), $l_{13}$ ($\hat{l}_{13}$), $l_{14}$, $l_{16}$, and $l_{17}$ are all lengths of line intervals between the node’s location (or, its current destination) and a certain intersection point along the node’s moving direction, with or without recharging at the charger. For instance, $l_{2}$ is the distance from the node’s location $z$ straight to the boundary of the (disk-shaped) movement area along its moving direction (at angle $\theta^{\prime}_{0}$), in the case of no need (or no opportunity) for recharging, while $l_{16}$ is the distance from $z$ to its original destination $w_{1}$ assuming that the node just departed from the charger, after recharging, for $w_{1}$. Detailed illustration of the relevant line segments and other variables can be found in [8]. Nevertheless, the density function given by Equation 1 is too complicated to derive a closed-form expression for it. In fact, it might be computationally more convenient to perform simulations than attempt a numerical solution altogether! Yet, both simulation and numerical evaluation are undesirable when the density function is just a subroutine repetitively called from a larger optimization engine to solve a spatial placement problem.

The question we are trying to answer is, what learning technique could we employ to derive, from a small training set, an accurate node density distribution across the domain of interest (a disk in our case)? The parameters of this function will be the $d$, $d^{\prime}$, $r$ and $\theta$ as presented above. If indeed such a form is learned as a mixture of component distributions, we could also, depending on the learning technique used, tease out its constituent components, to gain further insights to the density distribution [8], especially if the components come from elementary well-understood functions, e.g., from Gaussians. To this end, we next focus on putting the research question in the context of learning mixtures of distributions.

The general framework of density learning based on mixtures of Gaussian (or Möbius for that matter) distributions can be found in Figure 2, which is inspired by the mixture density network (MDN) first introduced in [9]. The basic idea of MDNs is to approximate an unknown probability distribution conditioned on the input, through a weighted sum of $K\in\mathbb{Z}_{\geq 1}$ distributions that are more manageable.

Note that Gaussian distributions are the most commonly assumed kernels in MDNs, since it has been proved that any density function in question can be approximated as accurately as possible by a Gaussian mixture model that is properly parameterized [16]. Despite the high flexibility and generality of Gaussian distributions, however, it remains an open question as to whether and how new distributions can be used to be better tailored (e.g., at a lower cost) to the needs of case-specific density approximations. In this paper, we answer the question of matching the finitely-supported density observed in (charging aware) mobile nodes by focusing and further elaborating on the (beta type III) Möbius distribution, which is introduced to substitute for the Gaussians.

As seen in Figure 2, the input to the learning models is the 4-tuple characterizing an instance, i.e., $(d,d^{\prime},r,\theta)$, corresponding to the energy availability ($d$) of the node and the radius ($d^{\prime}$) and location ($r$, $\theta$) of the charger (we start with the single-charger case). Through a fully connected deep neural network (containing two hidden layers of size $8K$ in the example), a number of $6K$ parameters, for Gaussian distributions (see Section 3.1), or $5K$ parameters, for Möbius distributions (see Section 3.2), are produced at the output layer.

Regardless of whether we train for Gaussian or Möbius mixtures, the training loss function (Equation (10)) is evaluated over a set of density values sampled from the continuous derived distribution, at coordinate points defined by the set $\mathcal{S}$. While we are free to define set $\mathcal{S}$ arbitrarily, our work is sampling uniformly on the x- and y- axis, along a $100\times 100$ grid.¹¹1That is, $\left|\mathcal{S}\right|=10^{4}$, following the same setting as in [8]. Only points within the circular support are of relevance to the density while points outside of it have zero density. For the loss function, we define

where $\hat{f}(z)$ denotes the density mixture obtained based on the learned parameters while $f(z)$ is the ground truth acquired either by simulation or by numerical evaluation of the density in Equation 1. It is worth noting that we are not using the mean squared error (MSE) nor the Kullback–Leibler (KL) divergence for the loss measure, despite their prevalence in the literature. Instead, we consider the mean “quartic” error, i.e., the error of $\left(\hat{f}(z)-f(z)\right)$ is raised to the power of 4, so as to enhance the ability to capture the “outlier”, which in this case is the peak of density around the charger’s location, as shown in Figures 1 and 2. MSE and KL divergence are used only as metrics of relative performance comparison in the experiments (see Section 4).

The results shown here are the outcome of ten separate training sessions using a training set of only five examples. A summary of the training results across training sessions is shown in Table 1 where the columns indicate the use of a different number, $K$, of distributions mixed. The row Params. lists the number of determining parameters needed by each model, which is $5K$ for Gaussian (mean and variance for x- and y-axis, as well as covariance) and $4K$ for Möbius ($a$, $\beta$, $\gamma$, and $\mu$), and the numbers in parentheses is the number inclusive of the additional weight parameters when $K>1$. The results are shown both in terms of MSE and KL-divergence, evaluated over the set $\mathcal{S}$, accompanied by 95% confidence intervals. Given that the training data set is the same, the benefit of using Möbius is obvious under a KL-divergence metric; a distinction in terms of MSE is less obvious. The improved KL score suggests a better matching of the “shape” of the distribution using Möbius.

The same results are plotted in Figure 5 albeit the errorbars indicate the range from minimum to maximum performance across the ten independent training runs. In terms of variability of training, the KL-divergence again points to a better albeit more variable outcome for Möbius at small $K$. Yet, for $K>4$ the KL-divergence range of results is almost always non-intersecting with the range of results from the corresponding Gaussian. Again, the MSE results do not allow for a more crisp separation between Gaussian and Möbius.

A demonstration that the training has been effective can be witnessed in the prediction results shown in Figure 6 which presents the average MSE and KL-divergence for the distribution prediction for 25 different parameter settings not present in the training set. The results appear similar, and occasionally better compared to the training results in terms of KL-divergence for Möbius (for $K>4$) over Gaussian mixtures. The good observed generalization performance suggests that the models do not overfit to the, small as it is, training data set.

The inability of the Gaussian mixture to mimic the density distribution can be explained by considering the shape of the constituent distributions. As an example, the Gaussian distributions ($K=3$) plotted in the two first rows of Figure 7 are superposed to produce the result shown in Figure 3 for a point belonging to the training set. While the peak at the charger’s location is captured, that is the limit of how interpretable the Gaussian results can be. The remaining background “dome” of the distribution is poorly reconstructed by the two remaining Gaussians both of which are haphazard in their placement. While the MSE of this mixture may indeed be small, its shape (and its KL-divergence) leaves a lot to be desired, lacking the symmetries we would expect to see in a better quality reconstruction.

The Möbius distributions ($K=3$), on the other hand, for the exact same training set input, corresponding to the bottom two rows of Figure 7 combined to produce Figure 4, are quite capable of capturing the requisite symmetries and be very interpretable. One Möbius distribution corresponds to the charger peak feature, and another one corresponds to the main RWP dome feature, with additional distributions contributing to enhance either one of the two features. Albeit it is also evident that at $K=3$ is the cusp of sufficiently good Möbius-based reconstruction as the error bars in Figure 5 attest. To this end we provide examples of the worst Möbius reconstructions for $K=3$ collected during our runs for points in the training set.

Figure 8 is for $(0.2,0.2,0.6,\pi)$ exhibiting an MSE of 0.0177 and KL-divergence of 0.1228. Figure 9, for the same parameters shows an MSE of 0.0276 and KL-divergence of 0.6863. Figure 8 captures an instance of a failure to capture the peak at the charger while able to capture a significant part of the mass on the main dome (Figure 8). Figure 9 represents the opposite behavior, where the density around the charger is captured (Figure 9) but not the main dome. As we comment next, the training phase for Möbius can be fragile owing to the interplay of the parameters and sub-expressions that define it.

The problem under consideration falls under the general category of inverse problems, i.e., problems where the observed data are assumed to come from a function whose parameters need to be derived. Solving inverse problems using learning algorithms has met with some success in various domains in the natural sciences [26, 27, 28]. The shape of the density function seen in [8] and [15] suggest a rather well-behaving, continuous and smooth function with the notable exception of the singular point at the exact location of the charger(s). Nevertheless, it is difficult to find a tractable density function that fits the shape while being analytic, and we resort to mixture models, although the choice of a suitable kernel distribution is still non-trivial. Note that our attention is given to inverse modeling, meaning the derivation of a closed form is considered as valuable, not only for accelerating the simulation process but also for facilitating further mathematical analyses of network performance (e.g., wireless capacity [15]). Hence, we disregard the class of sampling-based generative models (e.g., [29]) that can only numerically output the density value given the location.

To our understanding, such spatial density learning using finite-supported, non-Gaussian mixtures for improving mobility modeling and simulation is still lacking in the networking community. There are related studies on Gaussian-based mixture density learning for urban transportation, whose interests or assumptions are nevertheless at variance with our work. For instance, in order to alleviate traffic congestion in urban areas, the GPS traces of over 1,000 taxis in a city are analyzed in [18], where a Gaussian mixture model is adopted to model the taxis’ location distribution for hotspot evaluation. The mobility considered there is constrained, in the sense that the traces are limited to graph-like road networks. This is different from our assumption of a bounded circular movement area in the Euclidean plane. In [30], the patterns of cellphone users, who have transitional states of location driven by temporal periodicity and social interplay, are modeled by mixing pairs of Gaussians (accounting for periodicity), together with a power-law distribution (for social dynamics). In their example, the spatial model assumes a square-shaped support for mobility, but it is unclear how the potential density leakage is addressed. Also, the number of component distributions is supposedly fixed in [30]. In [31], Gaussian kernels with adaptive bandwidths (parameterizing the correlation matrices) are introduced for estimating spatial density of mobile users at different scales. The bandwidths are adapted by the kernel density estimation method to account for density discrepancy between highly (e.g., roads) and sparsely (e.g., deserts) populated areas. Like [18, 30], the focus is mainly on accommodating the heterogeneous dependence of density on locations, with no concern or solution for fitting specifically confined movement areas.

The application of non-Gaussian mixtures, with finite or semi-finite supports, appears to be more common in learning non-geographic density distributions. Such applications have occurred in various domains like astrophysics [32], economics [33], and power engineering [34]. For instance, beta distributions (with a support of $[0,1]$) are mixed in [32] to infer cosmological parameters of one-dimensional CDM (cold dark matter) models. The betas are also employed in [34] for an improved MDN such that the univariate wind power output of wind farms can be forecast with no issue of density leakage. In [33], an MDN using the multivariate Gamma distribution (having a support of $[0,+\infty)$ for each marginal variable) as its kernel is proposed to model the distribution of insurance claim amounts. The methods mentioned above are certainly not applicable to our scenario of mobility modeling for mobile networks.

The Möbius-based mixture learning model we have proposed is inspired by the work of [21], where a new class of disk-supported distributions (including the beta type III Möbius distribution) is proposed to model the bivariate density of wind speed and wind direction over Marion Island. Their work focuses on deriving fitted distribution forms for accurate wind modeling, making no mention of density mixture for potential fitness improvement.

Lastly, it is worth mentioning another type of mixture density learning problem that is often considered in autonomous driving and robot navigation, e.g., [35, 36, 37], which appear to be related but are in essence distinct from our scenario. There, the basic idea is to use mostly Gaussian and occasionally non-Gaussian (e.g., Laplacian [36]) mixture models to predict near-future trajectories of vehicles or robots in the locality (e.g., at an intersection) given a relatively short past history (e.g., within seconds), such that the moving agents can navigate efficiently and safely without collision. As can be seen, the main concern of such studies is prescriptive route planning or maneuver control at relatively small scales in both spatial and temporal aspects, or equivalently, involving one or at most a few pairs of random “waypoints” in the whole learning process. In contrast, the (disk-shaped) movement area we have assumed can be as broad as a city, and the span of a nodal trajectory simulated corresponds to the node moving for extended periods of time such that a steady state can be reached. Therefore, applying short-term prediction methods is at odds with our steady-state based problem formulation.

We have demonstrated the usefulness of Möbius distributions as component distributions for mixed density networks, as it pertains to capturing distributions of mobile nodes confined on disks. The results are promising in the sense that they produce small, interpretable models that retain the symmetries present in such distributions. Our objective is to use the learned model as a subroutine in an optimization task, e.g., a charger location assignment problem. The methodology we followed exposes also several shortcomings and areas for future improvement.

As mentioned in Section 4.1, the stability of general purpose optimizers like Adam for MDNs of Möbius distributions is relatively poor. We will pursue alternatives, e.g., based on particular forms of regularization [25]. Another extension will be towards distributions with other (non-disk) forms of support, e.g., rectangular or square, as they are also often seen in mobile network performance studies.

It is also worth pointing out that extensions can be applied also when some of the density-influencing parameters are stochastic, rather than constants. As an example, the distance $d$ which acts as a proxy of energy availability upon departure from a waypoint, can be replaced by distributions with a mean $d$, and there may be benefit to also include higher moments of such distributions.