Particle Swarm Optimisation

Particle Swarm Optimisation in Health Economic Modelling

Particle swarm optimisation searches for good parameter values or decisions by moving a group of trial solutions, called particles, through the space of possible values. Each particle steers by the best result it has found itself and by the best result found by the particles it communicates with, so the group gathers around promising regions. Health economists use it when a disease model has to be calibrated to observed data, or a simulation searched for the best allocation of vaccines or services, and the model can score a candidate but gives no gradient to follow. It is one branch of numerical optimisation and returns the best point found in the time allowed, which need not be the true optimum.

Origins in a simulated bird flock

Kennedy and Eberhart introduced the method in 1995 for optimising continuous nonlinear functions. It began as a simulation of a flock searching for food, in which each bird remembered the best position it had visited and also knew the best position any member had found. Velocities were nudged towards both points by random amounts. The random pull was later multiplied by 2, giving it a mean of 1, so that particles would overfly the target about half the time. A version that discarded the old velocity at every step proved ineffective at finding global optima, and the authors noted that the method needs only primitive mathematical operators.

The velocity and position updates

Bratton, in a doctoral thesis that defines a standard version, describes each particle by its current position, the best position it has found and its velocity, all originally drawn at random. The inertia weight form of the update is shown below.

$$v_{id}^{t+1} = w \, v_{id}^{t} + c_1 r_1 \left( p_{id} - x_{id}^{t} \right) + c_2 r_2 \left( g_{d} - x_{id}^{t} \right), \qquad x_{id}^{t+1} = x_{id}^{t} + v_{id}^{t+1}$$

where $x_{id}^{t}$ is the position of particle $i$ in dimension $d$ at step $t$, $v_{id}^{t}$ is its velocity, $p_{id}$ is the best position the particle has found, $g_{d}$ is the best position found by any particle in its neighbourhood, $w$ is the inertia weight, $c_1$ and $c_2$ are the weights on the personal and social pulls, and $r_1$ and $r_2$ are random numbers between 0 and 1 drawn afresh for each particle, dimension and step. In a calibration each dimension is one unknown model parameter.

The first term keeps a particle moving, the second pulls it towards its own best point and the third towards the best point its neighbours know. The original equations had no inertia weight and a constant of 2 on each pull. Bratton explains that velocities then had to be capped, because random combinations could otherwise drive them towards infinity and carry particles far outside the relevant region. Shi and Eberhart added the inertia weight in 1998 to control this growth.

Clerc and Kennedy formally published a constriction factor in 2002, derived from an analysis of particle trajectories and swarm stability.

$$\chi = \frac{2}{\left| 2 - \varphi - \sqrt{\varphi^{2} - 4\varphi} \right|}, \qquad \varphi = c_1 + c_2$$

where $\chi$ is the constriction factor, which multiplies the whole velocity update in place of the inertia weight, and $\varphi$ is the sum of the two pull weights, which must exceed 4 for quick and guaranteed convergence. With $c_1 = c_2 = 2.05$, $\varphi$ is 4.1, the square root is 0.640 and $\chi$ is 2 / 2.740 = 0.730. Bratton gives the standard value as about 0.72984, algebraically the same as an inertia weight of 0.72984 with pull weights of 1.49618.

Neighbourhoods and the balance between spread and speed

In the global best form every particle sees the best point found by the whole swarm, while in the local best ring each sees only its two neighbours. Bratton reports that the ring converges more slowly but is less likely to be trapped at a local optimum on problems with several peaks, while the global form usually does better where there is a single optimum. Kennedy and Eberhart saw the same tension, since a strong social pull sent the flock rushing prematurely towards local minima. In Bratton's benchmarks no swarm size between 20 and 100 particles was clearly better or worse for most problems, and his standard used 50.

Calibrating disease and transmission models

Natural history and transmission parameters that no study measures are set by calibration, searching for values whose outputs match observed data under a goodness-of-fit score. Vahdat and colleagues, calibrating a colorectal cancer microsimulation, list particle swarm optimisation with genetic algorithms and simulated annealing as heuristics used extensively for calibration, adding that their complexity has limited their use in real-world applications.

Michels, Kang and Wang calibrated three parameters of a spatially explicit agent-based model of influenza in Miami against weekly case data from the Florida Department of Health. A standard swarm of 20 particles was run 35 times for each setting and compared with Monte Carlo random sampling of 500, 750 and 1,000 parameter sets. The swarm achieved a mean error 9.4% to 20.5% lower even when it used less than half as many model runs. All particles in a step were evaluated in parallel, so a slow model's runs could be spread across processors.

Searching allocation and planning decisions

Yang and colleagues combined a particle swarm with a cross-regional SEIR model to allocate COVID-19 vaccines between regions, minimising total confirmed cases. They noted that the result had a certain black box quality that made stable rules hard to extract, and that economic costs were not considered. An allocation for an economic evaluation would need an objective in health and cost terms and a budget, which makes it a problem of constrained optimisation.

Hirpara and colleagues used a heuristic based on particle swarm optimisation to choose the number and location of trauma centres in Ohio, minimising two weighted kinds of mistriage. On 2012 data, optimising the network with one additional centre cut the objective by over 31.5%, and redistributing the existing 21 centres cut it by 30.4%, and the solutions were sensitive to the weights, as any weighted sum in multi-objective optimisation can be. The ISPOR Optimization Methods Emerging Good Practices Task Force lists swarm optimisation among heuristics that give approximate solutions where an exact approach is unavailable or computationally expensive.

Worked example calibrating two natural history parameters

All numbers in this example are illustrative. A cohort model has an unknown annual probability of disease onset $a$ and an unknown annual probability of death after diagnosis $q$. The targets are a 10-year cumulative incidence of 0.20 and a 5-year survival after diagnosis of 0.60, and the swarm minimises the sum of squared relative deviations from them.

$$f(a, q) = \left( \frac{1 - (1 - a)^{10} - 0.20}{0.20} \right)^{2} + \left( \frac{(1 - q)^{5} - 0.60}{0.60} \right)^{2}$$

where $f$ is the lack-of-fit score, lower being better, the first bracket compares modelled 10-year incidence with its target and the second compares modelled 5-year survival with its target. The settings are $w = 0.7$ and $c_1 = c_2 = 1.5$, rounded values close to the standard ones, with three particles sharing a global best.

1. Evaluate the starting swarm. Velocities start at zero and each particle's best point is its starting point.

For B, incidence is 1 minus 0.97 to the 10th power, or 0.2626, and survival is 0.95 to the 5th power, or 0.7738. The score is (0.0626 / 0.20) squared plus (0.1738 / 0.60) squared, or 0.0979 + 0.0839 = 0.1818, so B holds the global best.

2. Move A and C. Only the social term acts, since velocities are zero and each particle sits at its own best. For A, with social draws of 0.6 for $a$ and 0.4 for $q$, the velocity is 1.5 × 0.6 × 0.020 = 0.018 in $a$ and 1.5 × 0.4 × (minus 0.150) = minus 0.090 in $q$, moving A to (0.028, 0.110). For C, with social draws of 0.8 for $a$ and 0.5 for $q$, it is 1.5 × 0.8 × (minus 0.020) = minus 0.024 and 1.5 × 0.5 × (minus 0.070) = minus 0.0525, moving C to (0.026, 0.0675). B stays put.

3. Re-evaluate. A now scores 0.0606 (incidence 0.2472, survival 0.5584) and C scores 0.0556 (incidence 0.2316, survival 0.7051), so C becomes the global best.

4. Overshoot by momentum. In the second step C sits at both its own best and the global best, so only inertia acts. Its velocity becomes 0.7 × minus 0.024 = minus 0.0168 and 0.7 × minus 0.0525 = minus 0.03675, carrying it to (0.0092, 0.03075), where the score worsens to 0.4932. A and B, pulled towards C with all draws at 0.5, score 0.7118 and 0.0801, so the global best stays at 0.0556.

5. Pull back towards the remembered best. In the third step, with all draws at 0.5, both pulls on C point to (0.026, 0.0675). In $a$ each pull is 0.75 × (0.026 minus 0.0092) = 0.0126 and inertia gives 0.7 × minus 0.0168 = minus 0.01176, so the velocity is 0.0252 less 0.01176, or 0.01344. In $q$ each pull is 0.02756 and inertia gives minus 0.02572, so the velocity is 0.0294. C moves to (0.02264, 0.06015) and scores 0.0499, a new global best.

6. Interpret. The exact solution is $a = 1 - 0.8^{0.1} = 0.0221$ and $q = 1 - 0.6^{0.2} = 0.0971$. After three steps the best particle has $a$ close to this value while $q$ is still well short of it.

Settings that change the answer

Swarm size, the inertia weight or constriction factor, the pull weights, the neighbourhood and the handling of particles outside the permitted ranges all change where the swarm settles and how many model runs it takes. Michels and colleagues stopped after a fixed number of generations, the simplest rule, while other rules stop when particles crowd together or improvement stalls. The second ISPOR task force report states that settings such as the convergence level or the maximum number of iterations need to be justified.

Microsimulation and agent-based scores carry Monte Carlo error, so two runs at one point give different scores. A particle can record a lucky draw as its best point, so a remembered best merits re-scoring with more replications before it is trusted.

What the result means in an economic evaluation

A particle swarm returns one best-scoring point with no measure of its uncertainty. Bratton reports work showing that suitable settings guarantee convergence, though not necessarily to an optimal point. If one calibrated set is fixed in the model, its parameter uncertainty is left out of any later probabilistic sensitivity analysis unless several good sets are carried forward. Where several sets fit the targets about equally well, the swarm settles on one without showing that the others exist.

Particle swarms beside other search methods

Kennedy and Eberhart placed the method between genetic algorithms and evolutionary programming and likened the pull towards the two best points to crossover. In a particle swarm, though, candidates move through the space instead of being selected and bred, as in a genetic algorithm. Simulated annealing moves a single candidate and accepts some worse moves, the Nelder-Mead method is a local search, and Bayesian optimisation chooses each new run from a fitted surrogate, which suits very slow models. When the model supplies derivatives, gradient descent is the natural alternative. Wolpert and Macready's no free lunch theorems show that, averaged across all possible problems, any two search algorithms perform identically, so a swarm's success on one model is no guarantee of success on another.

Limitations and common misreadings

A frequent misreading treats the swarm's answer as the optimum. A swarm that converges early, especially in the global best form, can settle at a local optimum, and Michels and colleagues attributed outliers in their fully connected runs to exactly that. Repeated runs from different random starts are the usual check.

The method was designed for continuous parameters, so integer decisions such as numbers of facilities need an adapted version. Every particle at every step needs a full model run, and the answer is only as good as the objective, its weights and the data behind it.

Sources

• Bratton D. Simple and adaptive particle swarms. PhD thesis, Department of Computing, Goldsmiths College, University of London; 2011. https://research.gold.ac.uk/4752/
• Crown W, Buyukkaramikli N, Thokala P, et al. Constrained optimization methods in health services research: an introduction. Report 1 of the ISPOR Optimization Methods Emerging Good Practices Task Force. Value in Health. 2017;20(3):310-319. https://doi.org/10.1016/j.jval.2017.01.013
• Crown W, Buyukkaramikli N, Sir MY, et al. Application of constrained optimization methods in health services research: report 2 of the ISPOR Optimization Methods Emerging Good Practices Task Force. Value in Health. 2018;21(9):1019-1028. https://doi.org/10.1016/j.jval.2018.05.003
• Hirpara S, Vaishnav M, Parikh PJ, Kong N, Parikh P. Locating trauma centers considering patient safety. Health Care Management Science. 2022;25(2):291-310. https://doi.org/10.1007/s10729-021-09576-y
• Kennedy J, Eberhart R. Particle swarm optimization. Proceedings of ICNN'95, International Conference on Neural Networks. 1995;4:1942-1948. https://doi.org/10.1109/ICNN.1995.488968
• Michels A, Kang JY, Wang S. Particle swarm optimization for calibration in spatially explicit agent-based modeling. Journal of Artificial Societies and Social Simulation. 2022;25(2):8. https://doi.org/10.18564/jasss.4796
• Vahdat V, Alagoz O, Chen JV, Saoud L, Borah BJ, Limburg PJ. Calibration and validation of the Colorectal Cancer and Adenoma Incidence and Mortality (CRC-AIM) microsimulation model using deep neural networks. Medical Decision Making. 2023;43(6):719-736. https://doi.org/10.1177/0272989X231184175
• Wolpert DH, Macready WG. No free lunch theorems for optimization. IEEE Transactions on Evolutionary Computation. 1997;1(1):67-82. https://doi.org/10.1109/4235.585893
• Yang T, Deng W, Liu Y, Deng J. Comparison of health-oriented cross-regional allocation strategies for the COVID-19 vaccine: a mathematical modelling study. Annals of Medicine. 2022;54(1):941-952. https://doi.org/10.1080/07853890.2022.2060522