Retail Location Optimisation via Production-Constrained Spatial Interaction Modelling: A Case Study of Havering, London
Published:
This study develops a production-constrained spatial interaction model to evaluate two candidate supermarket locations in the London Borough of Havering. A Poisson log-linear GLM with exponential distance decay is calibrated on 775 Output Areas, using supermarket floor area as the destination attractiveness proxy. Site A, centrally located, consistently outperforms the peripheral Site B by approximately 11% across both small- and large-format scenarios. Sensitivity analysis incorporating a 74% oil price shock — triggered by the closure of the Strait of Hormuz in early 2026 — confirms Site A’s advantage is robust to rising transport costs, while Site B’s predicted patronage declines steadily.
1. Models and Calibration
1.1 Spatial Interaction Model Families
Spatial interaction models predict the volume of movement between origins and destinations. Wilson (1971) formalised the gravity model as a family sharing the general form:
\[T_{ij} = \text{balancing factors} \times \text{origin mass} \times \text{destination mass} \times \text{deterrence function}\]Unconstrained (Gravity) Model
Building on Zipf (1946)’s $P_1P_2/D$ formulation and Stewart (1948)’s demographic gravitation framework, the unconstrained gravity model takes the form $T_{ij} = k O_i^\alpha D_j^\gamma d_{ij}^\beta$. Wilson (1971) generalised the distance decay term to $f(d_{ij})$, producing a more flexible form:
\[T_{ij} = k O_i^\alpha D_j^\gamma f(d_{ij})\]where $k$ is a global scaling constant; $O_i$ and $D_j$ are the origin and destination masses; $f(d_{ij})$ is a distance decay function, taking two standard forms: the power law function $d_{ij}^{-\beta}$ and the exponential function $e^{-\beta d_{ij}}$; and $\alpha, \gamma$ control sensitivity to origin and destination size. All parameters ($k, \alpha, \gamma$, and the decay parameter $\beta$) are estimated freely, making it the least constrained model variant. It is applicable where neither origin outflows nor destination inflows are observed and must instead be approximated by proxy variables such as population or GDP (Grosche et al., 2007).
Note: In Wilson (1971)’s original formulation, the equation is expressed as $T_{ij} = kW_i^{(1)}W_j^{(2)}f(c_{ij})$, where $W_i^{(1)}$ and $W_j^{(2)}$ denote the mass terms of origin and destination, and $c_{ij}$ represents the generalised cost of travel. This study substitutes $O_i$, $D_j$, and $d_{ij}$ to denote origin mass, destination mass, and distance respectively, aligning with the variables used in the empirical analysis.
Production-Constrained Model
Where origin totals $O_i$ are known, the balancing factor $A_i$ ensures $\sum_j T_{ij} = O_i$ for all $i$ (Wilson, 1971):
\[T_{ij} = A_i O_i D_j^{\gamma} f(d_{ij}), \qquad A_i = \frac{1}{\sum_j D_j^{\gamma} f(d_{ij})}\]Free parameters reduce to $\gamma$ and $\beta$ only; destination inflows become a model output. It is widely used in retail geography, where residential population data provide reliable origin totals and the model allocates flows to competing destinations in proportion to their relative attractiveness (Lakshmanan & Hansen, 1965; Wilson, 1971; Harris & Wilson, 1978).
Attraction-Constrained Model
Where destination totals $D_j$ are known, the balancing factor $B_j$ ensures $\sum_i T_{ij} = D_j$ for all $j$ (Wilson, 1971):
\[T_{ij} = B_j O_i^{\alpha} D_j f(d_{ij}), \qquad B_j = \frac{1}{\sum_i O_i^{\alpha} f(d_{ij})}\]Free parameters are $\alpha$ and $\beta$. The model fixes known destination inflows and allocates them across origins, suited to contexts such as residential location modelling around fixed workplaces (Wilson, 1971; O’Kelly, 2009).
Doubly Constrained Model
When both $O_i$ and $D_j$ are observed, balancing factors are applied at both ends (Wilson, 1971):
\[T_{ij} = A_i O_i B_j D_j f(d_{ij})\] \[A_i = \frac{1}{\sum_j B_j D_j f(d_{ij})}, \qquad B_j = \frac{1}{\sum_i A_i O_i f(d_{ij})}\]Because $A_i$ and $B_j$ are mutually dependent, they are obtained by iterating these two equations until convergence (Wilson, 1971). With $\beta$ as the sole free parameter, this model yields the highest predictive accuracy and is standard for journey-to-work and transport demand forecasting where both origin outflows and destination inflows are known (Wilson, 1971; Flowerdew & Lovett, 1988; O’Kelly, 2009).
1.2 Model Selection, Calibration and Flow Conservation
A production-constrained model is selected because origin totals are known: the ONS Census provides fixed populations $O_i$ for 775 Output Areas (OAs) in the Borough of Havering, while the Geolytix dataset supplies 56 supermarket locations with floor areas $D_j$ as the destination attractiveness proxy, and shortest-path distances $d_{ij}$ are computed on the OSM pedestrian network.
Note: Output Areas (OAs) are the finest-grained census geography in England and Wales.
As shown in Figure 1, observed flows are heavily right-skewed, with most origin–destination pairs recording zero or near-zero trips. This violates the normality and homoscedasticity assumptions of OLS. A Poisson log-linear GLM is therefore adopted: its log-link ensures all predicted flows remain positive, and maximum-likelihood estimation is better suited to the discrete, non-negative nature of trip counts (Flowerdew & Aitkin, 1982).
The exponential function $e^{-\beta d_{ij}}$ imposes a constant proportional reduction per unit distance, favouring short-range interaction, whereas the power-law function $d_{ij}^{-\beta}$ decays more slowly and better accommodates long-range flows (Fotheringham & O’Kelly, 1989; Wilson, 2010; Chen, 2015). Thus, in this borough-level retail modelling, the exponential distance decay function is chosen. The performance also confirms it empirically (Table 1).
Table 1. Performance Comparison of Distance Decay Functions under Production-Constrained Model.
| Distance Decay | R² | RMSE |
|---|---|---|
| Power law ($d_{ij}^{-\beta}$) | 0.1164 | 4.120 |
| Exponential ($e^{-\beta d_{ij}}$) | 0.1172 | 4.119 |
Finally, the production-constrained Poisson regression model is linearised as:
\[\lambda_{ij} = \exp\!\left(\alpha_i + \gamma \ln D_j - \beta d_{ij}\right)\]Note: In this study, $\beta$ is constrained to be strictly positive, ensuring that interaction intensity decreases with distance.
where $\alpha_i$ is an origin fixed effect absorbing $\ln A_i + \ln O_i$. The remaining parameters $\gamma$ and $\beta$ are estimated by maximum likelihood, yielding $\hat{\gamma} = 0.1168$ and $\hat{\beta} = 2.43 \times 10^{-5}$.
Flow conservation is guaranteed structurally by the balancing factor $A_i$, defined as:
\[A_i = \frac{1}{\displaystyle\sum_j D_j^{\gamma} \cdot e^{-\beta d_{ij}}}\]which ensures that $\sum_j \hat{T}_{ij} = O_i$ for all $i$, confirmed empirically in Figure 2.
2. Scenarios
2.1 Assessment of Two Candidate Supermarket Locations
Two candidate sites in the London Borough of Havering are evaluated (Figure 3). Site A (E00011326), located centrally, sits adjacent to the borough’s principal population concentration but faces intense competition from numerous nearby supermarkets. Site B (E00011710), in the southern periphery, occupies a less contested market but draws from a substantially smaller population base, and its surrounding walking network is notably sparser. This demand accessibility versus competitive saturation trade-off cannot be resolved through descriptive analysis alone, but the production-constrained SIM captures it directly, as the balancing factor $A_i$ redistributes flows across competing destinations in proportion to their relative attractiveness and proximity.
The site assessment proceeds under two scenarios, each introducing a single new store into the destination set independently. The calibrated $\hat{\gamma}$ and $\hat{\beta}$ are held fixed from the calibration stage; only $A_i$ is recomputed to accommodate the expanded destination set, with flow conservation verified after each run. The results are shown in Table 2.
Note: Since each candidate site is represented by the centroid of its host OA, the network distance from that OA to the new supermarket is zero by construction. Although the exponential decay function does not produce the mathematical singularity that a power-law function would at zero distance, a zero-distance input still artificially inflates the predicted flow for the host OA. This distance is therefore replaced with a nominal value of 5 m for the host OA only; all other OA-to-site distances remain unchanged.
Under Scenario 1 (both 280 m²), Site A attracts 1,384 predicted flows compared to 1,243 for Site B, a difference of 141 flows. Under Scenario 2 (both 2,800 m²), Site A captures 1,801 flows versus 1,619 for Site B, with the absolute gap widening to 182 flows.
Table 2. Predicted total patronage directed to the new supermarket under each scenario.
| Scenario | Site A (E00011326) | Site B (E00011710) |
|---|---|---|
| Small (280 m²) | 1,384 | 1,243 |
| Large (2,800 m²) | 1,801 | 1,619 |
Site A is recommended across both scenarios. Both scenarios yield an approximately 11% advantage for Site A over Site B. Its mean network distance to all OA centroids is 4,806 m, roughly half that of Site B (9,201 m); under the calibrated distance decay function this difference is amplified, making Site A’s locational proximity the dominant factor that overwhelms the competitive saturation it faces.
Two principal limitations should be noted. Firstly, the calibrated model yields a low goodness-of-fit ($R^2 = 0.1172$), largely attributable to the coarse three-tier classification of supermarket floor area (280, 1,400, and 2,800 m²) used in the Geolytix dataset. This simplification of destination attractiveness propagates into the scenario analysis, as the same calibrated parameters govern the predicted flows. Secondly, although a walk network was chosen to match the fine-grained OA-level resolution and Havering is partially bounded by the Thames and extensive green belt which naturally limit cross-boundary retail interaction, the closed-system assumption remains the binding constraint on model realism.
2.2 Sharp Increase in Transport Cost
The US–Israel military campaign against Iran beginning in late February 2026 prompted Iran to close the Strait of Hormuz to international shipping. With approximately 20% of globally traded oil transiting this waterway, the disruption caused Brent crude to surge from roughly $64 to over $111 per barrel, a rise of approximately 74% (EIA, 2026; Figure 4).
In Wilson (1971)’s entropy-maximising derivation, $\beta$ governs the strength of distance decay as the Lagrange multiplier of the total travel cost constraint, including cost of time and money. Within a single borough, however, travel time differences between OD pairs are relatively limited, and the monetary cost of travel, principally fuel, becomes the dominant source of variation in generalised cost. An oil price shock therefore acts on this dominant component. In Havering, 78.47% of households own at least one car (highest rate in London), making the effect particularly direct. The adjusted $\beta$ is obtained as:
\[\beta_{\text{new}} = \beta_0 \cdot \left(1 + \phi \cdot \frac{\Delta P_{\text{oil}}}{P_{\text{oil},0}}\right)\]where $\phi$ is the fuel share of personal transport operating expenditure and $\Delta P_{\text{oil}} / P_{\text{oil},0}$ is the proportional oil price change. Taking $\phi \approx 0.59$ (ONS, 2025) and $\Delta P_{\text{oil}} / P_{\text{oil},0} \approx 0.74$ gives $\beta_{\text{new}} \approx 1.44\,\beta_0$.
Under the production constraint, balancing factors are recomputed to ensure flow conservation; increasing $\beta$ redistributes flows towards nearer destinations without altering origin totals. A supplementary analysis of 44,950 OD pairs confirms high correlation between walk and drive network distances (Pearson $r = 0.96$, $R^2 = 0.92$), indicating that network choice does not materially affect relative site performance.
At the current Brent price of $111/bbl (\(\beta = 1.44\beta_0\)), Site A attracts 1,810 and 1,390 predicted flows under the large- and small-format scenarios respectively, compared to 1,555 and 1,194 for Site B, consistent with the baseline findings. A broader sensitivity analysis across a range of oil prices (Figure 5) reveals a notable asymmetry: Site A’s flows remain near-flat as transport costs rise, whereas Site B exhibits a steady decline, losing flows more rapidly. These results reinforce the recommendation of Site A, whose performance is not only higher at baseline but also robust to external transport cost shocks.
2.3 15-Minute Walking Catchment
To estimate the population within a 15-minute walk of each site, the following procedure is applied (Figure 6): (a) OAs whose centroids fall within 2,000 metres (straight-line distance) of the site are selected; (b) within each selected OA, one random point per resident is generated assuming uniform spatial distribution; (c) each point is snapped to its nearest node on the OSM walking network; (d) a single-source Dijkstra algorithm from the site node determines whether each point falls within the 1,260-metre walking threshold, corresponding to a 15-minute walk at 1.4 m/s (Bosina & Weidmann, 2017). The process was repeated 30 times to converge random effects; Site A yields a mean accessible population of 5,000 and Site B 7,165.
References
Bosina, E., & Weidmann, U. (2017). Estimating pedestrian speed using aggregated literature data. Physica A: Statistical Mechanics and its Applications, 468, 1–29. https://doi.org/10.1016/j.physa.2016.09.044
Chen, Y. (2015). The distance-decay function of geographical gravity model: Power law or exponential law? Chaos, Solitons & Fractals, 77, 174–189. https://doi.org/10.1016/j.chaos.2015.05.022
Flowerdew, R., & Aitkin, M. (1982). A method of fitting the gravity model based on the Poisson distribution. Journal of Regional Science, 22(2), 191–202. https://doi.org/10.1111/j.1467-9787.1982.tb00744.x
Flowerdew, R., & Lovett, A. (1988). Fitting constrained Poisson regression models to interurban migration flows. Geographical Analysis, 20(4), 297–307. https://doi.org/10.1111/j.1538-4632.1988.tb00184.x
Fotheringham, A., & O’Kelly, M. E. (1989). Spatial interaction models: Formulations and applications. Springer.
Grosche, T., Rothlauf, F., & Heinzl, A. (2007). Gravity models for airline passenger volume estimation. Journal of Air Transport Management, 13(4), 175–183. https://doi.org/10.1016/j.jairtraman.2007.02.001
Harris, B., & Wilson, A. G. (1978). Equilibrium values and dynamics of attractiveness terms in production-constrained spatial-interaction models. Environment and Planning A: Economy and Space, 10(4), 371–388. https://doi.org/10.1068/a100371
Lakshmanan, J. R., & Hansen, W. G. (1965). A retail market potential model. Journal of the American Institute of Planners, 31(2), 134–143. https://doi.org/10.1080/01944366508978155
O’Kelly, M. E. (2009). Spatial interaction models. In R. Kitchin & N. Thrift (Eds.), International encyclopedia of human geography (pp. 365–368). Elsevier. https://doi.org/10.1016/B978-008044910-4.00529-0
Office for National Statistics. (2025). Family spending in the UK: April 2023 to March 2024. https://www.ons.gov.uk/releases/familyspendingintheukapril2023tomarch2024
Stewart, J. Q. (1948). Demographic gravitation: Evidence and applications. Sociometry, 11(1/2), 31–58. https://doi.org/10.2307/2785468
U.S. Energy Information Administration. (2026). Petroleum & other liquids: Spot prices. https://www.eia.gov/dnav/pet/pet_pri_spt_s1_d.htm
Wilson, A. G. (1971). A family of spatial interaction models, and associated developments. Environment and Planning A: Economy and Space, 3(1), 1–32. https://doi.org/10.1068/a030001
Wilson, A. (2010). Entropy in urban and regional modelling: Retrospect and prospect. Geographical Analysis, 42(4), 364–394. https://doi.org/10.1111/j.1538-4632.2010.00799.x
Zipf, G. K. (1946). The P1 P2/D hypothesis: On the intercity movement of persons. American Sociological Review, 11(6), 677–686. https://doi.org/10.2307/2087063