plays an important role in best fits. Furthermore, there is a big difference
between $M_{1}$ and $M_{2}$ at the best-fit point without $(g-2)_{\mu}$.

The most visible impact of the LHC 13-TeV constraints has been on the masses
of the strongly-interacting sparticles: see the left panels of Figs. 8 and 9
and compare the solid and dashed curves. On the other hand, the impact of the
LHC constraints on electroweak inos has been less marked: see Fig. 11. As was
to be expected, the importance of the $(g-2)_{\mu}$ constraint is seen in the
likelihood functions for charged slepton masses and electroweak inos: compare
the blue and green curves in Figs. 10 and 11. The composition of the LSP
$\tilde{\chi}^{0}_{1}$ is also different in the cases with and without
$(g-2)_{\mu}$: as seen in Fig. 12 and Table 6, a $\tilde{B}$ LSP is preferred
when $(g-2)_{\mu}$ is included, whereas a $\tilde{H}$ LSP is preferred when
$(g-2)_{\mu}$ is dropped. Moreover, the inclusion of the $(g-2)_{\mu}$
constraint also has significant indirect implications for the squark masses,
as also seen in Figs. 8 and 9. This analysis reinforces the importance of
clarifying the interpretation of the difference between the experimental
measurement and the SM calculation of $(g-2)_{\mu}$. We therefore welcome the
advent of the Fermilab $(g-2)_{\mu}$ experiment [111] and continued efforts to
refine the SM calculation.

We have also analyzed in this paper the importances of different mechanisms
for bringing the relic LSP density into the range favoured by Planck 2015 and
other data: see the shadings in Figs. 2, 4, 5, 6, 19, 20 and 21, and the
profile $\chi^{2}$ functions for the DM measures in Fig. 17. As we see there,
important roles are played by chargino coannihilation, slepton coannihilation
and rapid annihilation via direct-channel $H/A$ boson exchange, though other
mechanisms such as stau and squark coannihilation may be important in limited
regions of parameter space 303030Compared to the pMSSM7 analysis in [24], we
find that stop coannihilation is less prominent, and that rapid annihilation
through the $Z$ and the light Higgs boson is of very limited importance. In
these respects the more general realization of the MSSM with four additional
free parameters yields substantially different results.. In the case where the
$(g-2)_{\mu}$ constraint is dropped, there is a preference for a region where
$m_{\tilde{\chi}^{0}_{1}}\sim m_{\tilde{\chi}^{\pm}_{1}}\sim m_{\tilde{q}}\sim
m_{\tilde{g}}$ where multiple coannihilation processes play a role, and the
compressed spectrum reduces the sensitivity of the LHC sparticle searches.

In general, our analysis favours quite small deviations from the SM
predictions for electroweak, flavour and Higgs observables: see Figs. 12 and
14, in particular. We have also analyzed the pMSSM11 predictions for the
$A_{\Delta\Gamma}$ and $\tau(B_{s}\to\mu^{+}\mu^{-})$ observables recently
measured for the first time by the LHCb Collaboration [57]. As seen in Fig.
13, the pMSSM11 predictions for these observables are very similar to those in
the SM, deviating by much less than the current experimental uncertainties.
Accordingly, we do not include $A_{\Delta\Gamma}$ and
$\tau(B_{s}\to\mu^{+}\mu^{-})$ in our global fits.

We find that current LHC searches for long-lived particles do not impact our
scan of the pMSSM11 parameter space. However, the pMSSM11 still offers
significant prospects for the discovery of long-lived particles. When the
$(g-2)_{\mu}$ constraint is imposed, we find that $\Delta\chi^{2}\gtrsim 4$
for $\tau_{\rm NLSP}\gtrsim 10^{-10}$ s. However, when the $(g-2)_{\mu}$
constraint is dropped, values of $\tau_{\rm NLSP}$ as long as $10^{3}$ s (the
limit we impose in order to maintain successful Big Bang nucleosynthesis) are
allowed at the $\Delta\chi^{2}\lesssim 4$ level,

As seen in Figs. 20 and 21, the pMSSM11 offers interesting prospects for the
detection of supersymmetric DM. In both the spin-independent and -dependent
cases, cross sections close to the present experimental upper limits are
favoured at the 68% CL, whether or not $(g-2)_{\mu}$ is included in the set of
constraints. Interestingly, in the case of $\sigma^{\rm SI}_{p}$ with
$(g-2)_{\mu}$ included, there is a lower limit that is not far below the
neutrino ‘floor’ 313131However, we repeat that the uncertainties in the
calculation of $\sigma^{\rm SI}_{p}$ are large, and these remarks apply within
the framework of a calculation of $\sigma^{\rm SI}_{p}$ using SSARD., whereas
$\sigma^{\rm SI}_{p}$ may be much lower when $(g-2)_{\mu}$ is dropped, and low
values of $\sigma^{\rm SD}_{p}$ are allowed in both cases.

We turn finally to the prospects for discovering sparticles in future runs of
the LHC, or with a future linear $e^{+}e^{-}$ collider. As seen in Fig. 21,
whether or not $(g-2)_{\mu}$ is included in the global fit, the third-
generation squarks may well be within reach of future LHC runs, and the first-
and second-generation squarks and the gluino may also be accessible if the
$(g-2)_{\mu}$ constraint is dropped. If it is included, on the other hand,
there are also good prospects for discovering electroweakly-interacting
sparticles at an $e^{+}e^{-}$ collider, in particular the
${\tilde{e}},{\tilde{\mu}},\tilde{\chi}^{0}_{1},\tilde{\chi}^{0}_{2}$ and
$\tilde{\chi}^{\pm}_{1}$.

It is often said that the night is darkest just before dawn, and the same may
be true for supersymmetry.

## Acknowledgements

We thank Gino Isidori for useful discussions. The work of E.B. and G.W. is
supported in part by the Collaborative Research Center SFB676 of the DFG,
“Particles, Strings and the early Universe”. The work of K.S. is partially
supported by the National Science Centre, Poland, under research grants
DEC-2014/15/B/ST2/02157, DEC-2015/18/M/ST2/00054 and DEC-2015/19/D/ST2/03136.
K.S. thanks the TU Munich for hospitality during the final stages of this work
and has been partially supported by the DFG cluster of excellence EXC 153
“Origin and Structure of the Universe”, by the Collaborative Research Center
SFB1258. The work of M.B., I.S.F. and D.M.S. is supported by the European
Research Council via Grant BSMFLEET 639068. The work of J.C.C. is supported by
CNPq (Brazil). The work of M.J.D. is supported in part by the Australia
Research Council. The work of J.E. is supported in part by STFC (UK) via the
research grant ST/L000326/1 and in part via the Estonian Research Council via
a Mobilitas Pluss grant, and the work of H.F. is also supported in part by
STFC (UK). The work of S.H. is supported in part by the MEINCOP Spain under
contract FPA2016-78022-P, in part by the Spanish Agencia Estatal de
Investigación (AEI) and the EU Fondo Europeo de Desarrollo Regional (FEDER)
through the project FPA2016-78645-P, in part by the AEI through the grant IFT
Centro de Excelencia Severo Ochoa SEV-2016-0597, and by the Spanish MICINN
Consolider-Ingenio 2010 Program under Grant MultiDark CSD2009-00064. The work
of M.L. and I.S.F. is supported by XuntaGal. The work of K.A.O. is supported
in part by DOE grant de-sc0011842 at the University of Minnesota. The work of
G.W. is also supported in part by the European Commission through the
“HiggsTools” Initial Training Network PITN-GA-2012-316704. During part of this
work we used the middleware suite udocker [112] to deploy MasterCode on
clusters, developed by the EC H2020 project INDIGO-Datacloud (RIA 653549). We
are particularly grateful to Jorge Gomes for his kind support. We also thank
DESY and especially the DESY IT department for making us available the
computational resources of the BIRD/NAF2 cluster, which have been used
intensively to carry out this work.

## References

  * [1] For a compendium of CMS searches for supersymmetry, see https://twiki.cern.ch/twiki/bin/view/ CMSPublic/PhysicsResultsSUS.
  * [2] For a compendium of ATLAS searches for supersymmetry, see https://twiki.cern.ch/twiki/bin/view/ AtlasPublic/SupersymmetryPublicResults.
  * [3] X. Cui et al. [PandaX-II Collaboration], arXiv:1708.06917 [hep-ex], which updates A. Tan et al. [PandaX-II Collaboration], Phys. Rev. Lett. 117 (2016) no. 12, 121303 [arXiv:1607.07400 [hep-ex]].
  * [4] D. S. Akerib et al. [LUX Collaboration], Phys. Rev. Lett. 118 (2017) no.2, 021303 [arXiv:1608.07648 [astro-ph.CO]].
  * [5] C. Amole et al. [PICO Collaboration], arXiv:1702.07666 [astro-ph.CO].
  * [6] E. Aprile et al. [XENON Collaboration], arXiv:1705.06655 [astro-ph.CO].
  * [7] O. Buchmueller et al., Phys. Lett. B 657 (2007) 87 doi:10.1016/j.physletb.2007.09.058 [arXiv:0707.3447 [hep-ph]]; Eur. Phys. J. C 72 (2012) 1878 [arXiv:1110.3568 [hep-ph]]; Eur. Phys. J. C 72 (2012) 2243 [arXiv:1207.7315]; Eur. Phys. J. C 74 (2014) 2809 [arXiv:1312.5233 [hep-ph]].
  * [8] O. Buchmueller et al., Eur. Phys. J. C 74 (2014) 2922 [arXiv:1312.5250 [hep-ph]].
  * [9] O. Buchmueller et al., Eur. Phys. J. C 74 (2014) 12, 3212 [arXiv:1408.4060 [hep-ph]].
  * [10] K. J. de Vries et al., Eur. Phys. J. C 75 (2015) no.9, 422 [arXiv:1504.03260 [hep-ph]].
  * [11] A. M. Sirunyan et al. [CMS Collaboration], arXiv:1705.04650 [hep-ex].
  * [12] A. M. Sirunyan et al. [CMS Collaboration], arXiv:1705.04673 [hep-ex].
  * [13] A. M. Sirunyan et al. [CMS Collaboration], arXiv:1709.05406 [hep-ex].
  * [14] J. Billard, L. Strigari and E. Figueroa-Feliciano, Phys. Rev. D 89 (2014) no.2, 023524 doi:10.1103/PhysRevD.89.023524 [arXiv:1307.5458 [hep-ph]]; P. Cushman, C. Galbiati, D. N. McKinsey, H. Robertson, T. M. P. Tait, D. Bauer, A. Borgland and B. Cabrera et al., Snowmass Working Group Report: WIMP Dark Matter Direct Detection, arXiv:1310.8327 [hep-ex].
  * [15] H. Goldberg, Phys. Rev. Lett. 50 (1983) 1419; J. Ellis, J. Hagelin, D. Nanopoulos, K. Olive and M. Srednicki, Nucl. Phys. B 238 (1984) 453.
  * [16] G. Bertone, D. G. Cerdeno, M. Fornasa, R. Ruiz de Austri, C. Strege and R. Trotta, JCAP 1201 (2012) 015 doi:10.1088/1475-7516/2012/01/015 [arXiv:1107.1715 [hep-ph]]; A. Fowlie, A. Kalinowski, M. Kazana, L. Roszkowski and Y. L. S. Tsai, Phys. Rev. D 85 (2012) 075012 doi:10.1103/PhysRevD.85.075012 [arXiv:1111.6098 [hep-ph]]; C. Strege, G. Bertone, D. G. Cerdeno, M. Fornasa, R. Ruiz de Austri and R. Trotta, JCAP 1203 (2012) 030 doi:10.1088/1475-7516/2012/03/030 [arXiv:1112.4192 [hep-ph]]; L. Roszkowski, E. M. Sessolo and Y. L. S. Tsai, Phys. Rev. D 86 (2012) 095005 doi:10.1103/PhysRevD.86.095005 [arXiv:1202.1503 [hep-ph]]; P. Bechtle et al., JHEP 1206 (2012) 098 doi:10.1007/JHEP06(2012)098 [arXiv:1204.4199 [hep-ph]]; A. Fowlie, M. Kazana, K. Kowalska, S. Munir, L. Roszkowski, E. M. Sessolo, S. Trojanowski and Y. L. S. Tsai, Phys. Rev. D 86 (2012) 075010 doi:10.1103/PhysRevD.86.075010 [arXiv:1206.0264 [hep-ph]]; K. Kowalska, S. Munir, L. Roszkowski, E. M. Sessolo, S. Trojanowski and Y. L. S. Tsai, Phys. Rev. D 87 (2013) 115010 doi:10.1103/PhysRevD.87.115010 [arXiv:1211.1693 [hep-ph]]; C. Strege, G. Bertone, F. Feroz, M. Fornasa, R. Ruiz de Austri and R. Trotta, JCAP 1304 (2013) 013 doi:10.1088/1475-7516/2013/04/013 [arXiv:1212.2636 [hep-ph]]; L. Roszkowski, E. M. Sessolo and A. J. Williams, JHEP 1408 (2014) 067 doi:10.1007/JHEP08(2014)067 [arXiv:1405.4289 [hep-ph]]; P. Bechtle et al., Eur. Phys. J. C 76 (2016) no.2, 96 doi:10.1140/epjc/s10052-015-3864-0 [arXiv:1508.05951 [hep-ph]]; C. Han, K. i. Hikasa, L. Wu, J. M. Yang and Y. Zhang, Phys. Lett. B 769 (2017) 470 doi:10.1016/j.physletb.2017.04.026 [arXiv:1612.02296 [hep-ph]].
  * [17] P. Athron et al. [GAMBIT Collaboration], arXiv:1705.07935 [hep-ph].
  * [18] E. Bagnaschi et al., Eur. Phys. J. C 77 (2017) no.4, 268 [arXiv:1612.05210 [hep-ph]].
  * [19] E. Bagnaschi et al., Eur. Phys. J. C 77 (2017) no.2, 104 [arXiv:1610.10084 [hep-ph]].
  * [20] See, for example, C. F. Berger, J. S. Gainer, J. L. Hewett and T. G. Rizzo, JHEP 0902, 023 (2009) [arXiv:0812.0980 [hep-ph]]; S. S. AbdusSalam, B. C. Allanach, F. Quevedo, F. Feroz and M. Hobson, Phys. Rev. D 81, 095012 (2010) [arXiv:0904.2548 [hep-ph]]; J. A. Conley, J. S. Gainer, J. L. Hewett, M. P. Le and T. G. Rizzo, Eur. Phys. J. C 71, 1697 (2011) [arXiv:1009.2539 [hep-ph]]; J. A. Conley, J. S. Gainer, J. L. Hewett, M. P. Le and T. G. Rizzo, [arXiv:1103.1697 [hep-ph]]; S. S. AbdusSalam, et al., Eur. Phys. J. C 71 (2011) 1835 [arXiv:1109.3859 [hep-ph]]; S. Sekmen, S. Kraml, J. Lykken, F. Moortgat, S. Padhi, L. Pape, M. Pierini and H. B. Prosper et al., JHEP 1202 (2012) 075 [arXiv:1109.5119 [hep-ph]]; A. Arbey, M. Battaglia and F. Mahmoudi, Eur. Phys. J. C 72 (2012) 1847 [arXiv:1110.3726 [hep-ph]]; A. Arbey, M. Battaglia, A. Djouadi and F. Mahmoudi, Phys. Lett. B 720 (2013) 153 [arXiv:1211.4004 [hep-ph]]; M. W. Cahill-Rowley, J. L. Hewett, A. Ismail and T. G. Rizzo, Phys. Rev. D 88 (2013) 3, 035002 [arXiv:1211.1981 [hep-ph]].
  * [21] A. Fowlie, K. Kowalska, L. Roszkowski, E. M. Sessolo and Y. L. S. Tsai, Phys. Rev. D 88 (2013) 055012 doi:10.1103/PhysRevD.88.055012 [arXiv:1306.1567 [hep-ph]]; C. Strege, G. Bertone, G. J. Besjes, S. Caron, R. Ruiz de Austri, A. Strubig and R. Trotta, JHEP 1409 (2014) 081 doi:10.1007/JHEP09(2014)081 [arXiv:1405.0622 [hep-ph]]; M. Cahill-Rowley, J. L. Hewett, A. Ismail and T. G. Rizzo, Phys. Rev. D 91 (2015) no.5, 055002 [arXiv:1407.4130 [hep-ph]]; L. Roszkowski, E. M. Sessolo and A. J. Williams, JHEP 1502 (2015) 014 doi:10.1007/JHEP02(2015)014 [arXiv:1411.5214 [hep-ph]]; G. Bertone, F. Calore, S. Caron, R. Ruiz, J. S. Kim, R. Trotta and C. Weniger, JCAP 1604 (2016) no.04, 037 doi:10.1088/1475-7516/2016/04/037 [arXiv:1507.07008 [hep-ph]].
  * [22] R. K. Barman, B. Bhattacherjee, A. Choudhury, D. Chowdhury, J. Lahiri and S. Ray, arXiv:1608.02573 [hep-ph]; C. Boehm, P. S. B. Dev, A. Mazumdar and E. Pukartas, JHEP 1306 (2013) 113 doi:10.1007/JHEP06(2013)113 [arXiv:1303.5386 [hep-ph]].
  * [23] E. A. Bagnaschi et al., Eur. Phys. J. C 75 (2015) 500 [arXiv:1508.01173 [hep-ph]].
  * [24] P. Athron et al. [GAMBIT Collaboration], arXiv:1705.07917 [hep-ph].
  * [25] J. Ellis, F. Luo and K. A. Olive, JHEP 1509 (2015) 127 [arXiv:1503.07142 [hep-ph]].
  * [26] J. Ellis, J. L. Evans, F. Luo and K. A. Olive, JHEP 1602 (2016) 071 [arXiv:1510.03498 [hep-ph]].
  * [27] [ATLAS and CDF and CMS and D0 Collaborations], arXiv:1403.4427 [hep-ex].
  * [28] K. A. Olive et al. [Particle Data Group], Chin. Phys. C 38 (2014) 090001. doi:10.1088/1674-1137/38/9/090001
  * [29] F. Feroz and M.P. Hobson, Mon. Not. Roy. Astron. Soc. 384 (2008) 449 [arXiv:0704.3704 [astro-ph]]. F. Feroz, M.P. Hobson and M. Bridges, Mon. Not. Roy. Astron. Soc. 398 (2009) 1601-1614 [arXiv:0809.3437 [astro-ph]]. F. Feroz, M.P. Hobson, E. Cameron and A.N. Pettitt, [arXiv:1306.2144 [astro-ph]].
  * [30] For more information and updates, please see http://cern.ch/mastercode/. A survey of LHC Run 1 constraints can be found in K. J. de Vries, Global Supersymmetric Models after LHC Fits of Run 1 (PhD thesis, 2015), available on this website.
  * [31] P. Skands et al., JHEP 0407 (2004) 036 [arXiv:hep-ph/0311123];
  * [32] B. C. Allanach, Comput. Phys. Commun. 143 (2002) 305 [arXiv:hep-ph/0104145].
  * [33] S. Heinemeyer et al., JHEP 0608 (2006) 052 [arXiv:hep-ph/0604147]; S. Heinemeyer, W. Hollik, A. M. Weber and G. Weiglein, JHEP 0804 (2008) 039 [arXiv:0710.2972 [hep-ph]].
  * [34] S. Heinemeyer, W. Hollik, G. Weiglein and L. Zeune, JHEP 1312 (2013) 084 doi:10.1007/JHEP12(2013)084 [arXiv:1311.1663 [hep-ph]].
  * [35] S. Heinemeyer, W. Hollik and G. Weiglein, Comput. Phys. Commun. 124 (2000) 76 [arXiv:hep-ph/9812320]; S. Heinemeyer, W. Hollik and G. Weiglein, Eur. Phys. J. C 9 (1999) 343 [arXiv:hep-ph/9812472]; G. Degrassi, S. Heinemeyer, W. Hollik, P. Slavich and G. Weiglein, Eur. Phys. J. C 28 (2003) 133 [arXiv:hep-ph/0212020]; M. Frank et al., JHEP 0702 (2007) 047 [arXiv:hep-ph/0611326]; T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak and G. Weiglein, Comput. Phys. Commun. 180 (2009) 1426; T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak and G. Weiglein, Phys. Rev. Lett. 112 (2014) 14, 141801 [arXiv:1312.4937 [hep-ph]]; See http://www.feynhiggs.de .
  * [36] P. Athron, J. h. Park, D. St喞kinger and A. Voigt, Comput. Phys. Commun. 190 (2015) 139 [arXiv:1406.2319 [hep-ph]].
  * [37] P. Athron, M. Bach, D. Harries, T. Kwasnitza, J. h. Park, D. Stöckinger, A. Voigt and J. Ziebell, arXiv:1710.03760 [hep-ph]; E. Bagnaschi, G. F. Giudice, P. Slavich and A. Strumia, JHEP 1409 (2014) 092 doi:10.1007/JHEP09(2014)092 [arXiv:1407.4081 [hep-ph]]; E. Bagnaschi, J. Pardo Vega and P. Slavich, Eur. Phys. J. C 77 (2017) no.5, 334 doi:10.1140/epjc/s10052-017-4885-7 [arXiv:1703.08166 [hep-ph]].
  * [38] P. Athron, J. h. Park, T. Steudtner, D. Stöckinger and A. Voigt, JHEP 1701 (2017) 079 doi:10.1007/JHEP01(2017)079 [arXiv:1609.00371 [hep-ph]].
  * [39] H. Bahl and W. Hollik, Eur. Phys. J. C 76 (2016) 499 [arXiv:1608.01880 [hep-ph]]; H. Bahl, S. Heinemeyer, W. Hollik and G. Weiglein, arXiv:1706.00346 [hep-ph].
  * [40] G. Isidori and P. Paradisi, Phys. Lett. B 639 (2006) 499 [arXiv:hep-ph/0605012]; G. Isidori, F. Mescia, P. Paradisi and D. Temes, Phys. Rev. D 75 (2007) 115019 [arXiv:hep-ph/0703035], and references therein.
  * [41] F. Mahmoudi, Comput. Phys. Commun. 178 (2008) 745 [arXiv:0710.2067 [hep-ph]]; Comput. Phys. Commun. 180 (2009) 1579 [arXiv:0808.3144 [hep-ph]]; D. Eriksson, F. Mahmoudi and O. Stal, JHEP 0811 (2008) 035 [arXiv:0808.3551 [hep-ph]].
  * [42] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 185 (2014) 960 [arXiv:1305.0237 [hep-ph]], and references therein.
  * [43] Information about this code is available from K. A. Olive: it contains important contributions from J. Evans, T. Falk, A. Ferstl, G. Ganis, F. Luo, A. Mustafayev, J. McDonald, F. Luo, K. A. Olive, P. Sandick, Y. Santoso, C. Savage, V. Spanos and M. Srednicki.
  * [44] J. R. Ellis, A. Ferstl and K. A. Olive, Phys. Lett. B 481, 304 (2000) doi:10.1016/S0370-2693(00)00459-7 [hep-ph/0001005]; J. R. Ellis, K. A. Olive, Y. Santoso and V. C. Spanos, Phys. Rev. D 71, 095007 (2005) doi:10.1103/PhysRevD.71.095007 [hep-ph/0502001].
  * [45] J. R. Ellis, K. A. Olive and C. Savage, Phys. Rev. D 77 (2008) 065026 [arXiv:0801.3656 [hep-ph]].
  * [46] M. Muhlleitner, A. Djouadi and Y. Mambrini, Comput. Phys. Commun. 168 (2005) 46 [hep-ph/0311167].
  * [47] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, Eur. Phys. J. C 74 (2014) 2, 2711 [arXiv:1305.1933 [hep-ph]]; JHEP 1411 (2014) 039 [arXiv:1403.1582 [hep-ph]].
  * [48] P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein and K. E. Williams, Comput. Phys. Commun. 181 (2010) 138 [arXiv:0811.4169 [hep-ph]], Comput. Phys. Commun. 182 (2011) 2605 [arXiv:1102.1898 [hep-ph]]; P. Bechtle et al., Eur. Phys. J. C 74 (2014) 3, 2693 [arXiv:1311.0055 [hep-ph]]; P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, Eur. Phys. J. C 75 (2015) no.9, 421 [arXiv:1507.06706 [hep-ph]].
  * [49] J. de Blas [HEPFit Collaboration], talk given at the First FCC Physics Workshop, CERN, Jan. 2017, https://indico.cern.ch/event/550509/ contributions/2413788/attachments/ 1396663/2130440/1stFCC_Ph_W_deBlas.pdf.
  * [50] P. Athron et al., Eur. Phys. J. C 76 (2016) no.2, 62 doi:10.1140/epjc/s10052-015-3870-2 [arXiv:1510.08071 [hep-ph]].
  * [51] LEP Electroweak Working Group [ALEPH, CDF, D0, DELPHI, L3, OPAL and SLD Collaborations, LEP Electroweak Working Group, Tevatron Electroweak Working Group and SLD Electroweak and Heavy Flavour Groups], arXiv:1012.2367 [hep-ex].
  * [52] M. Baak [Gfitter Collaboration], PoS EPS -HEP2013 (2013) 203.
  * [53] D. Stockinger, J. Phys. G 34 (2007) R45 [arXiv:hep-ph/0609168]; J. Miller, E. de Rafael and B. Roberts, Rept. Prog. Phys. 70 (2007) 795 [arXiv:hep-ph/0703049]; J. Prades, E. de Rafael and A. Vainshtein, arXiv:0901.0306 [hep-ph]; F. Jegerlehner and A. Nyffeler, Phys. Rept. 477, 1 (2009) [arXiv:0902.3360 [hep-ph]]; M. Davier, A. Hoecker, B. Malaescu, C. Z. Yuan and Z. Zhang, Eur. Phys. J. C 66, 1 (2010) [arXiv:0908.4300 [hep-ph]]. J. Prades, Acta Phys. Polon. Supp. 3, 75 (2010) [arXiv:0909.2546 [hep-ph]]; T. Teubner, K. Hagiwara, R. Liao, A. D. Martin and D. Nomura, arXiv:1001.5401 [hep-ph]; M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 71 (2011) 1515 [arXiv:1010.4180 [hep-ph]].
  * [54] G. Bennett et al. [The Muon g-2 Collaboration], Phys. Rev. Lett. 92 (2004) 161802, [arXiv:hep-ex/0401008]; and Phys. Rev. D 73 (2006) 072003 [arXiv:hep-ex/0602035].
  * [55] V. Khachatryan et al. [CMS and LHCb Collaborations], Nature 522 (2015) 68 [arXiv:1411.4413 [hep-ex]]
  * [56] M. Aaboud et al. [ATLAS Collaboration], Eur. Phys. J. C 76 (2016) no.9, 513 [arXiv:1604.04263 [hep-ex]].
  * [57] R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 118 (2017) no.19, 191801 [arXiv:1703.05747 [hep-ex]].
  * [58] M. Misiak, H. M. Asatrian, R. Boughezal, M. Czakon, T. Ewerth, A. Ferroglia, P. Fiedler and P. GamBino et al., arXiv:1503.01789 [hep-ph].
  * [59] Y. Amhis et al. [Heavy Flavor Averaging Group (HFAG) Collaboration], arXiv:1412.7515 [hep-ex].
  * [60] B. Kronenbitter et al. [Belle Collaboration], Phys. Rev. D 92 (2015) no.5, 051102 [arXiv:1503.05613 [hep-ex]].
  * [61] T. Huber, T. Hurth and E. Lunghi, JHEP 1506 (2015) 176 [arXiv:1503.04849 [hep-ph]].
  * [62] A. J. Buras, P. GamBino, M. Gorbahn, S. Jager and L. Silvestrini, Nucl. Phys. B 592 (2001) 55 [hep-ph/0007313].
  * [63] W. J. Marciano, Phys. Rev. Lett. 93 (2004) 231803 [arXiv:hep-ph/0402299].
  * [64] C. Patrignani et al. [Particle Data Group], Chin. Phys. C 40 (2016) no.10, 100001 and 2017 update.
  * [65] A. J. Buras, D. Buttazzo, J. Girrbach-Noe and R. Knegjens, JHEP 1511 (2015) 033 [arXiv:1503.02693 [hep-ph]].
  * [66] A. V. Artamonov et al. [E949 Collaboration], Phys. Rev. Lett. 101 (2008) 191802 [arXiv:0808.2459 [hep-ex]].
  * [67] CMS Collaboration [CMS Collaboration], CMS-PAS-EXO-16-036.
  * [68] V. Khachatryan et al. [CMS Collaboration], Eur. Phys. J. C 75 (2015) no.7, 325 [arXiv:1502.02522 [hep-ex]].
  * [69] V. Khachatryan et al. [ CMS Collaboration], JHEP 1410 (2014) 160 [arXiv:1408.3316 [hep-ex]].
  * [70] P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak and G. Weiglein, Eur. Phys. J. C 75 (2015) 9, 421 [arXiv:1507.06706 [hep-ph]].
  * [71] The ATLAS collaboration [ATLAS Collaboration], ATLAS-CONF-2016-085.
  * [72] ATLAS Collaboration,
http://cds.cern.ch/record/2273866/ files/ATLAS-CONF-2017-050.pdf.

  * [73] P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 594 (2016) A13 [arXiv:1502.01589 [astro-ph.CO]].
  * [74] B. Borasoy and U. G. Meissner, Annals Phys. 254 (1997) 192 [hep-ph/9607432].
  * [75] J. M. Alarcon, L. S. Geng, J. Martin Camalich and J. A. Oller, Phys. Lett. B 730 (2014) 342 [arXiv:1209.2870 [hep-ph]].
  * [76] X. Z. Ling, X. L. Ren and L. S. Geng, arXiv:1710.07164 [hep-ph].
  * [77] M. Hoferichter, J. Ruiz de Elvira, B. Kubis and U. G. Mei?ner, Phys. Rev. Lett. 115 (2015) 092301 [arXiv:1506.04142 [hep-ph]].
  * [78] J. Ruiz de Elvira, M. Hoferichter, B. Kubis and U. G. Meissner, arXiv:1706.01465 [hep-ph].
  * [79] J. M. Alarcon, J. Martin Camalich and J. A. Oller, Phys. Rev. D 85 (2012) 051503 [arXiv:1110.3797 [hep-ph]].
  * [80] R. D. Young and A. W. Thomas, Phys. Rev. D 81 (2010) 014503 doi:10.1103/PhysRevD.81.014503 [arXiv:0901.3310 [hep-lat]]; J. Giedt, A. W. Thomas and R. D. Young, Phys. Rev. Lett. 103 (2009) 201802 doi:10.1103/PhysRevLett.103.201802 [arXiv:0907.4177 [hep-ph]]; H. W. Lin, arXiv:1112.2435 [hep-lat]; L. Alvarez-Ruso, T. Ledwig, J. Martin Camalich and M. J. Vicente-Vacas, Phys. Rev. D 88 (2013) no.5, 054507 doi:10.1103/PhysRevD.88.054507 [arXiv:1304.0483 [hep-ph]].
  * [81] J. R. Ellis, R. A. Flores and S. Ritz, Phys. Lett. B 198 (1987) 393.
  * [82] J. Ellis, K. A. Olive and V. C. Spanos, JCAP 1110 (2011) 024 [arXiv:1106.0768 [hep-ph]].
  * [83] M. G. Aartsen et al. [IceCube Collaboration], Eur. Phys. J. C 77 (2017) no.3, 146 [arXiv:1612.05949 [astro-ph.HE]].
  * [84] J. Ellis, K. A. Olive, C. Savage and V. C. Spanos, Phys. Rev. D 81 (2010) 085004 [arXiv:0912.3137 [hep-ph]].
  * [85] M. Drees, H. Dreiner, D. Schmeier, J. Tattersall and J. S. Kim, Comput. Phys. Commun. 187 (2015) 227 [arXiv:1312.2591 [hep-ph]].
  * [86] C. Bal痙s et al. [GAMBIT Collaboration], Eur. Phys. J. C 77 (2017) no.11, 795 [arXiv:1705.07919 [hep-ph]].
  * [87] E. Conte, B. Fuks and G. Serret, Comput. Phys. Commun. 184 (2013) 222 [arXiv:1206.1599 [hep-ph]].
  * [88] M. Papucci, K. Sakurai, A. Weiler and L. Zeune, Eur. Phys. J. C 74 (2014) no.11, 3163 [arXiv:1402.0492 [hep-ph]].
  * [89] S. Kraml, S. Kulkarni, U. Laa, A. Lessa, W. Magerl, D. Proschofsky-Spindler and W. Waltenberger, Eur. Phys. J. C 74 (2014) 2868 [arXiv:1312.4175 [hep-ph]].
  * [90] F. Ambrogi et al., arXiv:1701.06586 [hep-ph].
  * [91] W. Beenakker, C. Borschensky, M. Krﾃ､mer, A. Kulesza, E. Laenen, S. Marzani and J. Rojo, Eur. Phys. J. C 76 (2016) no.2, 53 [arXiv:1510.00375 [hep-ph]]; W. Beenakker, M. Kramer, T. Plehn, M. Spira and P. M. Zerwas, Nucl. Phys. B 515 (1998) 3 [hep-ph/9710451];
  * [92] E. A. Bagnaschi, M. Papucci, K. Sakurai, A. Weiler and L. Zeune, in preparation.
  * [93] G. Aad et al. [ATLAS and CMS Collaborations], JHEP 1608 (2016) 045
[arXiv:1606.02266 [hep-ex]].

  * [94] CMS Collaboration, CMS-PAS-EXO-16-036.
  * [95] T. Jittoh, K. Kohri, M. Koike, J. Sato, T. Shimomura and M. Yamanaka, Phys. Rev. D 76 (2007) 125023 [arXiv:0704.2914 [hep-ph]]; K. Jedamzik, Phys. Rev. D 77 (2008) 063524 [arXiv:0707.2070 [astro-ph]] and JCAP 0803 (2008) 008 [arXiv:0710.5153 [hep-ph]]; T. Jittoh, K. Kohri, M. Koike, J. Sato, T. Shimomura and M. Yamanaka, Phys. Rev. D 78 (2008) 055007 [arXiv:0805.3389 [hep-ph]]; T. Jittoh, K. Kohri, M. Koike, J. Sato, K. Sugai, M. Yamanaka and K. Yazaki, Phys. Rev. D 84 (2011) 035008 [arXiv:1105.1431 [hep-ph]]; R. H. Cyburt, J. Ellis, B. D. Fields, F. Luo, K. A. Olive and V. C. Spanos, JCAP 1212 (2012) 037 [arXiv:1209.1347 [astro-ph.CO]]; M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu, arXiv:1709.01211 [hep-ph].
  * [96] A. Buckley, Eur. Phys. J. C 75, no. 10, 467 (2015) [arXiv:1305.4194 [hep-ph]].
  * [97] Hunter, J. D., Computing In Science & Engineering, 9.3 (2007), 90-95.
  * [98] W. G. Hollik, JHEP 1608 (2016) 126 doi:10.1007/JHEP08(2016)126 [arXiv:1606.08356 [hep-ph]].
  * [99] J. E. Camargo-Molina, B. O’Leary, W. Porod and F. Staub, Eur. Phys. J. C 73 (2013) no.10, 2588 doi:10.1140/epjc/s10052-013-2588-2 [arXiv:1307.1477 [hep-ph]].
  * [100] S. Profumo and A. Provenza, JCAP 0612 (2006) 019 doi:10.1088/1475-7516/2006/12/019 [hep-ph/0609290].
  * [101] G. Angloher et al. [CRESST Collaboration], Eur. Phys. J. C 76 (2016) no.1, 25 doi:10.1140/epjc/s10052-016-3877-3 [arXiv:1509.01515 [astro-ph.CO]].
  * [102] R. Agnese et al. [SuperCDMS Collaboration], Phys. Rev. Lett. 116 (2016) no.7, 071301 doi:10.1103/PhysRevLett.116.071301 [arXiv:1509.02448 [astro-ph.CO]].
  * [103] L. T. Yang et al. [CDEX Collaboration], arXiv:1710.06650 [hep-ex].
  * [104] B. J. Mount et al., arXiv:1703.09144 [physics.ins-det].
  * [105] E. Aprile et al. [XENON Collaboration], JCAP 1604 (2016) no.04, 027 doi:10.1088/1475-7516/2016/04/027 [arXiv:1512.07501 [physics.ins-det]].
  * [106] K. C. Y. Ng, J. F. Beacom, A. H. G. Peter and C. Rott, arXiv:1703.10280 [astro-ph.HE].
  * [107] K. Choi et al. [Super-Kamiokande Collaboration], Phys. Rev. Lett. 114 (2015) no.14, 141301 doi:10.1103/PhysRevLett.114.141301 [arXiv:1503.04858 [hep-ex]].
  * [108] H. Baer et al., arXiv:1306.6352 [hep-ph].
  * [109] G. Moortgat-Pick et al., Eur. Phys. J. C 75 (2015) no.8, 371 [arXiv:1504.01726 [hep-ph]].
  * [110] M. J. Boland et al. [CLIC and CLICdp Collaborations], arXiv:1608.07537 [physics.acc-ph].
  * [111] J. Grange et al. [Muon g-2 Collaboration], Fermilab Muon $(g-2)_{\mu}$ Technical Design Report, arXiv:1501.06858 [physics.ins-det].
  * [112] J. Gomes et al., arXiv:1711.01758 [cs.SE].


# Optimal Scheduling of Variable Speed Pumps using Mixed Integer Linear
Programming - Towards an Automated Approach

Tomasz Janus
Department of Mechanical, Aerospace and Civil Engineering
The University of Manchester
Sackville Street, Manchester, M13 9PL, UK
<EMAIL_ADDRESS>
&Bogumil Ulanicki
School of Engineering and Sustainable Development
De Montfort University
The Gateway, Leicester LE1 9BH, UK
<EMAIL_ADDRESS>
&Kegong Diao
School of Engineering and Sustainable Development
De Montfort University
The Gateway, Leicester LE1 9BH, UK
<EMAIL_ADDRESS>
Visiting ResearcherEmeritus Professor of Engineering SystemsSenior Lecturer

###### Abstract

This article describes the methodology for formulating and solving optimal
pump scheduling problems with variable-speed pumps as mixed integer linear
programs using piece-linear approximations of the network components. The
water distribution network (WDN) is simulated with an initial pump schedule
for a defined time horizon, e.g. 24 hours, using a nonlinear algebraic solver.
Next, the network element equations including VSPs are approximated with
linear and piece-linear functions around chosen operating point(s). Finally, a
fully parameterized MILP is formulated in which the objective is the total
pumping cost. The method was programmed in MATLAB/OCTAVE and Python and is
publicly available on GitHub 111https://github.com/tomjanus/milp-scheduling.
The method was used to solve a pump scheduling problem on a a simple two
variable speed pump single-tank network that allows the reader to easily
understand how the methodology works and how it is applied in practice. The
obtained results showed that the formulation is robust and the optimizer is
able to return global optimal result in a reliable manner for a range of
operating points. The work summarized here is a prototype of a framework that
is being implemented in a Python package for automated solution of optimal
pump scheduling problems on EPANET networks using mixed integer linear
programming.

_Keywords_ mixed integer linear programming $\cdot$ energy optimization
$\cdot$ global optimization $\cdot$ variable speed pumps $\cdot$ water
distribution networks

##  1 Introduction

Pump scheduling plays a crucial role in optimizing WDN operation and has been
an important subject of research over the last decades. The purpose of pump
scheduling is to find the sequence of pump ON/OFF statuses and pump speeds, in
case of variable-speed pumps (VSPs), that minimize the pumping cost by
shifting pumping to time periods with lower electricity tariffs and by routing
the flow through the network such that the energy is used most efficiently.
Optimal pump schedules depend on the electricity tariff profile, the demand
profile, and the hydraulic characteristic of the network. Pump scheduling is
an important problem in the operation of WDNs because up to 70% of total
operating costs are attributed to electricity consumption for pumping. Past
research and industrial case studies have shown that optimization of pump
schedules can lead to up to 10-20% reduction in pumping costs, i.e. up to
7%-14% of total operating costs of WDNs.

Nevertheless, finding optimal pump schedules is a difficult dynamic mixed-
integer nonlinear nonconvex optimization problem (Bonvin et al., 2017).
Nonconvexity arises from the nonlinearities imposed by the network equations
that act as constraints. Integer variables used for selecting statuses of
pumps and other active network components such as valves, as well as, in case
of piece-wise approximations and relaxations, selection of active domains,
make the problem combinatorial. Dynamics arise in storage tanks and require
the optimization problem to consider variables and constraints from all time-
steps, which in turn substantially increases the problem size. Consequently,
pump scheduling problems are inherently difficult to solve and, depending on
the size of the model and on the optimization method, may get stuck at local
minima, fail to find a global optimum within the allocated time or at all, or
fail to find a feasible solution altogether.

Optimization methods, including those applied to pump scheduling, can be
broadly classified into two categories: (a) deterministic mathematical
programming, such as mixed integer nonlinear programming, mixed integer linear
programming, etc. and (b) stochastic evolutionary searches, such as e.g.
genetic algorithms, particle swarm optimization (PSO), and ant colony
optimization (ACO). The former use full or partial information about the model
to guide the optimization process whilst the latter explore the objective
space without reliance on the information about the model. Consequently,
mathematical programming tends to be significantly faster at the cost of being
more difficult to formulate, while evolutionary algorithms tend to be slow but
relatively easy to set up and able to work in conjunction with an input-output
model of any mathematical form. Slower convergence speeds however, are
leveraged by some GAs’ massive parallelization capabilities, as demonstrated
in Reed and Hadka (2014). Additionally, many GAs support multiple objectives
while multi-objective extensions within mathematical programming frameworks
are far less common. Last but not least, evolutionary algorithms, including
GAs, are very good at exploring large decision spaces. Unfortunately, they
cannot guarantee that the global optimal solution has been found nor provide
bounds on global optimality (Menke et al., 2016) which is in contrast to
(convex) mathematical programming which can provide such guarantees.

The literature on the subject of pump scheduling is extensive and voluminous.
For a more in-depth study the readers are referred to the most recent review
papers on the topic by Mala-Jetmarova et al. (2017) and Wu et al. (2018).
Here, we shall only mention a handful of aspects of pump scheduling that
pertain to mixed integer programming which has recently become more popular,
most likely in response to the recent improvements in speed, reliability,
scalability and stability of numerical solvers, such as CPLEX (Cplex, IBM
ILOG, 2009), GUROBI (Gurobi Optimization, LLC, 2023), MOSEK (MOSEK, 2023),
SCIP (Bestuzheva et al., 2021)

To increase the speed and robustness of finding globally optimal solution, the
originally non-convex mixed integer nonlinear program (MINLP) formulation of
the pump scheduling problem needs to be numerically simplified. The literature
distinguishes between (a) model simplification/reduction that is applied to
the model before optimization problem formulation and (b) simplification of
the optimization problem using various mathematical techniques from the
operational research (OR) field. Model network simplifications were described
in Anderson and Al-Jamal (1995); Deuerlein (2008); Alzamora et al. (2014) and
in the context of real-time pump scheduling, by Shamir and Salomons (2008).
The most popular methods for problem simplification are via (a) various
relaxations and approximations of constraints, objective(s) and the type of
decision variables in order to turn the original problem it into a convex
nonlinear or linear problem, (b) decomposition into several easier to solve
problems, and (c) relaxation of the optimality criterion see e.g. Gleixner et
al. (2012).

Non-convexity was addressed in Fooladivanda and Taylor (2018); Singh and
Kekatos (2019) who turned the non-convex MINLP problem into a convex MINLP
problem via different second order cone (SOC) relaxations. Additionally,
Fooladivanda and Taylor (2018) incorporated VSPs and pressure reducing valves.
Bonvin et al. (2017) addressed the non-convexity in a less formal and more
heuristic way that required restricted formulations supporting only the
network topologies without loops.

Problem decomposition into short-term and long-term optimization was
investigated in Pulido-Calvo and Gutierrez-Estrada (2011). Similarly, Ulanicki
et al. (2007) proposed time decomposition via solution of a relaxed continuous
problem to find optimum reservoir trajectories, followed by a solution of a
mixed-integer pump scheduling problem that tracks those trajectories.
Lagrangian decomposition and Benders decomposition were successfully applied
to MINLP problems by Ghaddar et al. (2015) and Naoum-Sawaya et al. (2015),
respectively. Alternatively, the problem could also be decomposed spatially by
dividing the network into smaller isolated subsystems.

Various mixed-integer problem reductions via approximations and relaxations
were recently performed by a number of authors. Vieira et al. (2020) used an
iterative approach with a feedback loop from EPANET simulator to iteratively
limit the error from MILP with component relaxations. In a conceptually
similar way, Liu et al. (2020) created a MILP formulation with component
relaxations which are tightened by the solutions of a series of EPANET
simulations. The authors claimed that their formulation with pipe
characteristic relaxations over-performs one with piece-linear approximations
due to reduction of binary variables. Salomons and Housh (2020) tested
different levels of reduction of binary variables in order to align
computational times of MILP pump scheduling schemes for real-time pump
scheduling applications. Bonvin et al. (2021) developed relaxation of non-
convex constraints of the original non-convex MINLP problem using Polyhedral
Outer Approximations (OA) and solved the relaxed convex problem with branch
and bound method for convex MINLPs. Their method supports VSPs. Most recently,
Tasseff et al. (2022) developed tight polyhedral relaxations of the original
MINLP, derived novel cuts using duality theory, added novel optimization-based
bound tightening and cut generation procedures and implemented their method in
an open-source Julia package (Tasseff et al., 2019) The authors addressed two
main deficiencies of current state-of-the-art MILP solvers: the slow
improvement of dual bounds and the difficulty in generating feasible primal
solutions. The authors considered fixed-speed pumps only.

Although many advancements in pump scheduling using mixed integer programming
have been introduced recently, the literature on the subject is still
fragmented with papers addressing some crucial aspects of the methodology
whilst omitting others. Additionally, some treatment of crucial network
components such as VSPs is under-represented. Meanwhile, the findings reported
by many authors suggest that current state-of-the-art MILP solvers are able to
find optimal pump schedules for networks with around $100$ nodes Liu et al.
(2020), i.e. of sufficient complexity to make them practical This suggests
that an automated method for solving pump scheduling problems with mixed
integer programming could be of practical value to the community. In this
paper we communicate our initial findings in prototyping a method for
automatic network conversion into a MILP problem and subsequent solution using
one of the available solvers. We provide a complete mathematical description
of the MILP formulation of a network with VSPs, which is novel. The validity
of the approach and its reliability and robustness is tested on a small
network with two parallel VSPs and one tank. This work is a part of a larger
project on automating pump scheduling with mixed integer linear programming
that is being developed in the `dev-python` branch of the GitHub repository of
MILOPS-WDN - the Mixed Integer Linear Optimal Pump Scheduler (Janus and
Ulanicki, 2023). The prototype source code used in this study is contained in
the `main` branch of the same repository.

##  2 Methodology

Figure 1: Block diagram visualising the methodology for formulating and
solving pump scheduling problems using mixed integer linear programming
approach with linear network element approximations.

The procedure described in this paper follows the methodology shown in Fig. 1.
First, the network is simulated using initial schedules in order to obtain
approximate operating points that are used in the subsequent approximations of
the model components and the pump power consumption characteristic. The
nonlinear equations of pipe and pump characteristics are approximated using
linear and piecewise linear approximations so that the network model can be
put into a linear form required by MILP. Other nonlinear components such as
valves, including check valves and PRVs, as well as other nonlinearities such
as leakage could additionally be included in the model formulation. For
simplicity, these are not considered in this study, but the methodology does
not prohibit their inclusion. Subsequently, the MILP optimization problem of
the standard form shown in Eq. 1 is solved for the desired time horizon of
typically $24$ hours. Finally, the optimal schedule of pump ON/OFF statuses
and pump speeds is input back into the simulator and the final simulation
result is compared against the result of the initial simulation.

$\displaystyle\min_{\mathbf{x}}$ $\displaystyle\mathbf{c}^{T}\,\mathbf{x}$ (1)
s.t. $\displaystyle\mathbf{A}_{ineq}\,\mathbf{x}\leq\mathbf{b}_{ineq}$
$\displaystyle\mathbf{A}_{eq}\,\mathbf{x}=\mathbf{b}_{eq}$
$\displaystyle\mathbf{l}\leq\mathbf{x}\leq\mathbf{u}$ $\displaystyle
x_{i}\in\mathbb{Z},\forall_{i}\in\varUpsilon$

The objective function $\mathbf{c}^{T}\,\mathbf{x}$describes the total pumping
cost over the optimization time-horizon. The network equations are described
with inequality and equality constraints using with two tuples:
$\left(\mathbf{A}_{ineq},\mathbf{b}_{ineq}\right)$ and
$\left(\mathbf{A}_{eq},\mathbf{b}_{eq}\right)$. The lower bounds $\mathbf{l}$
and the upper bounds $\mathbf{u}$ enforce physical limits on the decision
variable vector $\mathbf{x}$, such as e.g. minimum and maximum tank levels,
maximum pump flows, etc. $\varUpsilon$ is a nonempty subset of the set
$\\{1\ldots n\\}$ that specifies the indices of the integer variables, where
$n=|\mathbf{x}|$. Integer variables are used for the selection of pumps and
active segments in piece-linear approximations. The decision variable vector
$\mathbf{x}$ includes all state variables of the network model in all time
moments plus the auxiliary (artificial) variables. The auxiliary variables are
introduced to represent the variables bound to the domains of piece-linear
segments obtained via piece-linear linear approximations. By convention, an
auxiliary (continuous) variable for $x$ has a symbol $xx$ with an index
representing the index of the subdomain, e.g. $xx_{i}$ represents the value of
$x$ if $x$ lies inside segment $i$. By definition $x=\sum_{i=1}^{m}xx_{i}$
where $m$ is the number of segments. Auxiliary variables are accompanied by
binary selection integer variables. By definition, these are denoted by
capital letters, e.g. $XX$, and satisfy equation $\sum_{i=1}^{m}XX_{i}=1$,
which means that only one component can be active at a time.

###  2.1 Network equations

Water network is described with standard network equations describing: (1)
headlosses in pipes, (2) flow continuity equations in nodes, incl. tanks, (3)
pump characteristics describing pumping head and power consumption vs. flow,
and (4) head vs. volume relationships in tanks. More information about
modelling of WDNs can be found in Strafaci et al. (2007). The equations are
briefly listed below as their understanding is necessary to follow the model
approximation steps.

####  2.1.1 Pipe headloss equations

The nodal pressures in the network are modelled as follows:

$\mathbf{R}\left|\mathbf{q}(k)\right|\mathbf{q}(k)+\mathbf{\Lambda}_{c}^{T}\mathbf{h}_{c}(k)+\mathbf{\Lambda}_{f}^{T}\,\mathbf{h}_{f}(k)=\mathbf{0}$
(2)

where $\mathbf{\Lambda}_{c}$ and $\mathbf{\Lambda}_{f}$ are node-element
incidence matrices for the connection nodes and the fixed nodes, respectively.
$\mathbf{R}$ are pipe resistances, $\mathbf{q}$ is the vector of element
flows, $\mathbf{h}_{c}$ are the calculated heads and $\mathbf{h}_{f}$ are the
fixed heads, e.g. heads in tanks and reservoirs. The above equation can be
split into multiple equations, each representing a single pipe $j$ for
$j\in\langle 1,\ldots,|\mathbb{P}|\rangle$ where $\mathbb{P}$ denotes the set
of pipes.

$\underbrace{R_{j}\left|q_{j}(k)\right|q_{j}(k)}_{\textrm{pipe
characteristic}}+\mathbf{\Lambda}_{c,j}^{T}\mathbf{h}_{c}(k)+\mathbf{\Lambda}_{f,j}^{T}\,\mathbf{h}_{f}(k)=0$
(3)

where $\mathbf{\Lambda}_{c,j}^{T}\mathbf{h}_{c}(k)=h_{d,j}(k)$ and
$\mathbf{\Lambda}_{f,j}^{T}\,\mathbf{h}_{f}(k)=h_{o,j}(k)$, i.e. the
downstream and the upstream head of pipe $j$, respectively.

####  2.1.2 Mass balance in nodes

Mass balance in nodes for each time step $k\in{1\ldots K}$ is calculated as

$\mathbf{\Lambda}_{c}\,\mathbf{q}(k)-\mathbf{d}(k)=\mathbf{0}$ (4)

where $\mathbf{d}$ is the matrix of nodal demands.

####  2.1.3 Pump power consumption

Pump power consumption is modelled with the relationship described in Ulanicki
et al. (2008) representing power demand of a group of $n$ identical pumps,
each operating at speed $s$.

$P(q,n,s)=n\,s^{3}\,P\left(\frac{q}{n\,s}\right)$ (5)

where

$P\left(\frac{q}{n\,s}\right)=a_{3}\left(\frac{q}{n\,s}\right)^{3}+a_{2}\left(\frac{q}{n\,s}\right)^{2}+a_{1}\left(\frac{q}{n\,s}\right)+a_{0}$
(6)

The coefficients $a_{3}$, $a_{2}$, $a_{1}$, $a_{0}$ are unique for each
individual pump model. Smaller individual VSPs could alternatively be modelled
with the scaling of Sárbu and Borza (1998) as advised in Simpson and Marchi
(2013), although this decision is left to the user.

####  2.1.4 Pump hydraulics

Pump hydraulics are formulated with the pump characteristic model $H=H(q,n,s)$
from Ulanicki et al. (2008) that describes the relationship between the head
gain $H$ and the pump group flow $q$ for a group of $n$ identical pumps
operating at speed $s$.

$\frac{H}{n^{2}\,s^{2}}=A\,\left(\frac{q}{n\,s}\right)^{2}+B\left(\frac{q}{n\,s}\right)+C$
(7)

which translates into

$H=A\,q^{2}+B\,q\,n\,s+C\,n^{2}\,s^{2}$ (8)

Eqs. 5, 6, and 8 are implemented in the `dev-battery` branch of the EPANET
repository (Janus and Ulanicki, 2020).

####  2.1.5 Tank model

Tanks are modelled with a backward finite-difference equation describing the
change in tank head $h_{t}$ between time steps $k$ and $k-1$ as a function of
the net flow $q_{t}(k)$ in/out of the tank in time steps $k=1\ldots K$, where
$K$ is the optimization time horizon.

$h_{t}(k)-h_{t}(k-1)-\frac{1}{A_{t}}\,q_{t}(k)=0\quad\forall k=2\ldots K$ (9)

with initial condition

$h_{t}(1)-h_{t,init}=0$ (10)

For cylindrical tanks, the tank’s surface area $A_{t}=\textrm{const}$. For
other geometries, $A_{t}=A_{t}(h_{t})$ and needs to be included in the model.
If this relationship is not linear, it needs to be approximated with one of
the approximation methods - see below.

###  2.2 Linear and piece-linear approximations of the network components

When approximating nonlinear functions with linear functions, the choice of
the linearization technique is often left to the user. The possible choices
are: (1) tangent line approximation using first two terms of Taylor’s
expansion, (2) linearization by substitution via introduction of additional
variables and transformations, (3) piecewise linearization, (4) convex hull
approximation. The choice of the method should be based on the specific
characteristics of the nonlinear function and the context of the problem. Each
method has its own limitations and applicability and, in the context of pump
scheduling, will affect the accuracy of the solution and the complexity of the
MILP formulation. In the following sections, the choice of linearization
techniques was made by the Authors but the methodology is not limited to those
choices and the readers are encouraged to try different techniques in order to
fine-tune their problem formulations.

####  2.2.1 Linear approximation of the pump power consumption model

Pump power consumption is linearized by finding a tangent line to the power
consumption model of the group of $n$ parallel pumps given in Eq. 5 and Eq. 6
at the linearization point $(q_{0},s_{0})$. Since each pump is linearized
individually, the linearized equations are derived for $n=1$. Linearization
using a tangent at the nominal point is chosen over other linearization
methods due to the fact that power consumption curves tend to be quite flat.
Therefore, linearizing the curve around the nominal operating point should
offer sufficient approximation accuracy whilst keeping the complexity at
minimum, e.g. compared to piece-wise linearization. Linearization of Eq. 5 and
Eq. 6 for $n=1$ around the selected operating point $(q_{0},s_{0})$ yields

$P(q,s)=P(q_{0},s_{0})+3\,a_{3}\,q_{0}^{2}\,\delta
q\,+2\,a_{2}\,q_{0}\,s_{0}\,\delta q+a_{2}\,q_{0}^{2}\,\delta
s\,+a_{1}\,s_{0}^{2}\,\delta q+2\,a_{1}\,q_{0}\,s_{0}\,\delta
s\,+3\,a_{0}\,s_{0}^{2}\,\delta s$ (11)

After grouping similar terms with $\delta s$ and $\delta q$

$P(q,s)=P(q_{0},s_{0})+\left(3\,a_{3}\,q_{0}^{2}+2\,a_{2}\,q_{0}\,s_{0}+a_{1}\,s_{0}^{2}\right)\delta
q+\left(a_{2}\,q_{0}^{2}+2\,a_{1}\,q_{0}\,s_{0}+3\,a_{0}\,s_{0}^{2}\right)\delta
s$ (12)

where $\delta q=q-q_{0}$ and $\delta s=s-s_{0}$. Using short notation
$P(q,s)=P$ and $P(q_{0},s_{0})=P_{0}$ we obtain the following equation for the
difference between the power consumption at point $(q,s)$ and the power
consumption at the linearization point $(q_{0},s_{0})$.

$P^{\prime}=P-P_{0}=\left(3\,a_{3}\,q_{0}^{2}+2\,a_{2}\,q_{0}\,s_{0}+a_{1}\,s_{0}^{2}\right)\,(q-q_{0})+\left(a_{2}\,q_{0}^{2}+2\,a_{1}\,q_{0}\,s_{0}+3\,a_{0}\,s_{0}^{2}\right)\,(s-s_{0})$
(13)

Ultimately, the linearized power consumption of a pump at any time step $k$ is
calculated with Eq. 14.

$P(q(k),s(k))=P(q_{0},s_{0})+P^{\prime}(q(k),s(k))$ (14)

Throughout this paper $(q_{0},s_{0})=(q_{n},s_{n})$, i.e. the power is
linearized around the operating point at the nominal speed $s_{n}=1.0$ and
flow $q_{n}$ equal to the flow at the nominal speed for which the pump
efficiency is at its maximum. After collecting all terms, Eq. 13 can be
simplified to Eq. 15. The expressions for coefficients $m^{q}_{j}$,
$m^{s}_{j}$ and $c_{j}$ can be found in the Appendix in Section 5.1.

$P_{j}(k)=m^{q}_{j}\,q_{j}(k)+m^{s}_{j}\,s_{j}(k)+c_{j}$ (15)

Eq, 15 holds only if the pump status is ON. This is enforced in MILP by
expressing the equality as a double sided inequality with the ‘big U’ trick as
expressed in Eq. 16. $U_{power}$ is a large number in the order of magnitude
but larger in value than the largest power consumption calculated by the model
and $n_{j}(k)$ is the status of $j$-th pump at timestep $k$.

$(n_{j}(k)-1)\,U_{power}\leq
m^{q}_{j}\,q_{j}(k)+m^{s}_{j}\,s_{j}(k)+c_{j}-P_{j}(k)\leq(1-n_{j}(k))\,U_{power}$
(16)

If $n_{j}(k)=1$, Eq. 16 is reduced to equality as both sides of the inequality
are zero. If $n_{j}(k)=0$, $(n_{j}(k)-1)\,U_{power}=-U_{power}$ and
$(1-n_{j}(k))\,U_{power}=+U_{power}$. Consequently, the equality in Eq. 15 is
not enforced. Note that Eq. 15 yields non-zero power consumption for $q=0$ and
$s=0$ as an unwanted byproduct of linearization. To ensure that power is null
when the pump is OFF, $P_{j}(k)$ is made to obey the following two-sided
inequality that forces $P_{j}(k)=0$ for $n_{j}(k)=0$.

$0\leq P_{j}(k)\leq n_{j}(k)\,U_{power}$ (17)

####  2.2.2 Piece-linear approximation of the pipe model

Nonlinear pipe characteristic can approximated by $n_{s}$ piece-linear
segments by defining $n_{s}-1$ breakpoints and exploiting the symmetry of the
characteristic around the origin $(0,0)$. A generic pipe characteristic and
its piece-linear form for $n_{s}=3$ segments :
$[\tilde{q}_{-2},\tilde{q}_{-1}]$, $[\tilde{q}_{-1},\tilde{q}_{1}]$ and
$[\tilde{q}_{1},\tilde{q}_{2}]$, is shown in Figure 2.

Figure 2: A generic form of a quadratic pipe characteristic model describing
headloss $\Delta h$ across the pipe vs flow $q$ and its piece-wise
linearization with $3$ linear segments.

Each segment $i$ of pipe $j$ is described with a linear equation $\Delta
h_{j}=m_{j,i}\,q_{j}+c_{j,i}$ where

$m_{j,i}=\frac{\Delta h_{j,i}-\Delta
h_{j,i-1}}{\tilde{q}_{j,i}-\tilde{q}_{j,i-1}}$ (18)

and

$c_{j,i}=\frac{\Delta h_{j,i-1}\,\tilde{q}_{j,i}-\Delta
h_{j,i}\,\tilde{q}_{j,i-1}}{\tilde{q}_{j,i}-\tilde{q}_{j,i-1}}$ (19)

The breakpoint $(\tilde{q}_{j,1},\Delta h_{j,1})$ is fixed for every pipe $j$
and taken from the simulator, e.g. from the model state at 12.00 o’clock. The
breakpoint $(\tilde{q}_{j,2},\Delta h_{j,2})$ is chosen arbitrarily for every
pipe $j$ and should be made large enough to cover the entire range of observed
pipe flows. The points $(\tilde{q}_{j,-1},\Delta h_{j,-1})$ and
$(\tilde{q}_{j,-2},\Delta h_{j,-2})$ do not need to be calculated but instead,
can be derived by exploiting the symmetry of the pipe’s characteristic around
point $(0,0)$. In order to represent the three linear segments within one pipe
model, two new types of auxiliary variables are required: a continuous
variable $ww_{j,i}$ for the flow in pipe $j$ and segment $i$, and a binary
variable $BB_{j,i}$ for selecting the active segment of the piece-wise
linearized pipe characteristic.

$ww_{j,i}=\left\\{\begin{array}[]{ll}q_{j},&\mbox{if
$q_{j}\in[\tilde{q}_{j,i-1},\tilde{q}_{j,i}]$}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$BB_{j,i}=\left\\{\begin{array}[]{ll}1,&\mbox{if
$q_{j}\in[\tilde{q}_{j,i-1},\tilde{q}_{j,i}]$}.\\\
0,&\mbox{otherwise}.\end{array}\right.$

The piece-linear pipe model is represented by the following 4 relationships.

$BB_{j,i}(k)\,\tilde{q}_{j,i-1}\leq ww_{j,i}(k)\leq
BB_{j,i}(k)\,\tilde{q}_{j,i}$ (20)
$q_{j}(k)-\sum_{i=1}^{n_{s,pipe}}ww_{j,i}(k)=0$ (21)
$\sum_{i=1}^{n_{s,pipe}}\,BB_{j,i}(k)=1$ (22)
$h_{o}^{j}(k)-h_{d}^{j}(k)-\sum_{i=1}^{n_{s,pipe}}\left(m_{j,i}\,ww_{j,i}(k)+c_{j,i}\,BB_{j,i}(k)\right)=0$
(23)

where $h_{o}^{j}(k)-h_{d}^{j}(k)=\Delta
h_{j}(k)=\Lambda_{c,j}^{T}h_{c}(k)+\Lambda_{f,j}^{T}\,h_{f}(k)$

Eq. 20 ‘binary linearizes’ the flow with respect to segment selection. If the
binary segment selection variable $BB_{j,i}$ is zero, $ww_{j,i}$ is forced to
become zero. Otherwise, $ww_{j,i}$ is bound between $\tilde{q}_{j,i-1}$ and
$\tilde{q}_{j,i}$. Eq. 21 relates the pipe flow variable $q$ to the auxiliary
flow variable $ww$. Eq. 22 ascertains that only one segment in the piece-
linear equation is active at a time. Finally, Eq. 23 describes the linear pipe
headloss.

####  2.2.3 Piece-linear approximation of the pump characteristic - VSPs

Pump characteristic is approximated with piece-linear surfaces. A piece linear
approximation of a surface can be constructed in different ways with different
number of linear segments such as in the case of piece-linear curve
linearization, but also with different segment geometries. In this paper, a
piece-linear approximation of the pump characteristic is a top surface of a
polyhedron whose sides are defined by vertices
$p_{1}=(s_{min},q_{min},H_{1})$, $p_{2}=(s_{max},q_{min},H_{2})$,
$p_{3}=(s_{max},q_{int}^{s_{max}},H_{3})$,
$p_{4}=(s_{min},q_{int}^{s_{min}},H_{4})$, and the nominal point
$p_{n}=(s_{n},q_{n},H_{n})$ \- see Fig. 3. The nominal point is derived for
the nominal speed $s_{n}=1.0$ and the maximum efficiency flow
$q_{n}=q^{\eta_{max}}$. $s_{min}$ and $s_{max}$ are the minimum and the
maximum allowed pump speeds, respectively. $q_{min}=0$ is the minimum pump
flow and $q_{int}^{s_{min}}$ and $q_{int}^{s_{max}}$ are the intercept flows,
i.e. the flows for which the pump head $H=0$ at the minimum and at the maximum
pump speed, respectively. Consequently, $H_{3}=0$ and $H_{4}=0$.

Figure 3: A generic pump characteristic $H=H(q,s)$ with its piece-linear
approximation (left) and projection of the piece-linear approximation onto the
$q-s$ plane (right).

The top surface of the polyhedron is defined by the equations of four
Euclidean planes $A_{i}=\\{(q,s,H)\,|\,(q,s)\in\Delta_{i},H\in\mathbb{R}\\}$.
For each plane $A_{i}$, the point $(q,s)$ must lie within the triangular
domain $\Delta_{i}$ defined by the plane’s projection onto the $q-s$ space -
see Fig 3 (right). The equations of each of the four planes $A_{i}$ are
derived from the three vertices of the polyhedron that are contained within
it. The domain for each plane is defined by three inequality constraints
derived from the equations of the lines containing the three sides of its
projection - see Section 5.2 of the Appendix. Direction of each inequality is
indicated in Fig 3 with a small blue arrow. The projections of points $p_{1}$,
$p_{2}$, $p_{3}$, $p_{4}$, and $p_{n}$ onto the $q-s$ space are denoted as
$p^{\prime}_{1}$, $p^{\prime}_{2}$, $p^{\prime}_{3}$, $p^{\prime}_{4}$, and
$p^{\prime}_{n}$, respectively.

Introduction of piece-linear pump characteristic requires three auxiliary
variables: a binary variable $AA_{j,i}(k)$, and two continuous variables
$ss_{j,i}(k)$, $qq_{j,i}(k)$ for each of the four domains, each pump and each
time-step. $AA_{j,i}(k)$ defines whether the current operating point
$(s_{j}(k),q_{j}(k))$ of pump $j$ lies within the $i$-th segment of the
linearized pump characteristic: 1 for YES and 0 for NO. $ss_{j,i}(k)$ and
$qq_{j,i}(k)$ are the speed and the flow of pump $j$ in each domain $i$ of the
piece-linear pump characteristic approximation, respectively. The
approximation is defined as follows.

$s_{j}(k)-\sum_{i=1}^{n_{s,pump}}ss_{j,i}(k)=0$ (24)
$q_{j}(k)-\sum_{i=1}^{n_{s,pump}}qq_{j,i}(k)=0$ (25)

Eq. 24 and Eq. 25 link the auxiliary segment speeds and segment flows to the
original pump speeds and pump flows, respectively. Only one segment is allowed
to be active if the pump is ON, i.e. $n_{j}=1$ Otherwise, if the pump is OFF,
i.e. $n_{j}=0$, no segments are allowed to be active.

$\sum_{i=1}^{n_{s,pump}}AA_{j,i}(k)-n_{j}(k)=0$ (26)

The binary segment selection variable $AA_{j,i}(k)$ is used to ‘binary
linearize’ the pump speed and the pump flow with respect to segment selection.

$AA_{j,i}(k)\,s_{j,min}\leq ss_{j,i}(k)\leq AA_{j,i}(k)\,s_{j,max}$ (27)
$0\leq qq_{j,i}(k)\leq AA_{j,i}(k)\,q_{j,max}$ (28)

If $AA_{j,i}=0$, then $ss_{j,i}$ and $qq_{j,i}$ are forced to be zero.
Otherwise, $ss_{j,i}$ and $qq_{j,i}$ are bound with box constraints between
$s_{j,min}$ and $s_{j,max}$, and $0$ and $q_{j,max}=q_{int,max}$, respectively
- see Fig 3 (right). In summary, pump speeds $s_{j}(k)$ and pump flows
$q_{j}(k)$ are forced to be equal to the sums of the auxiliary variables
$ss_{j,i}(k)$ and $qq_{j,i}(k)$ where, at most one auxiliary segment variable
is active at any time-step $k$. The linearized pump characteristic is
represented with the following formula using the ‘big-U’ trick.

$-U^{\prime}\leq\Delta
h_{j}(k)-\sum_{i=1}^{i=4}\left(dd_{j,i}\,ss_{j,i}(k)+ee_{j,i}\,qq_{j,i}(k)+ff_{j,i}\,AA_{j,i}(k))\right)\leq
U^{\prime}$ (29)

where $U^{\prime}=(1-n_{j}(k))U_{pump}$ and the coefficients $dd_{j,i}$,
$ee_{j,i}$, $ff_{j,i}$ describe the equation of plane $A_{i}$. Derivation of
plane equations is described in Appendix in Section 5.2. Each triangular
domain $\Delta_{i}$ of segment $A_{i}$ is defined with three inequality
constraints:

$\begin{bmatrix}m_{qq}^{(1)}&m_{ss}^{(1)}&c^{(1)}\\\
m_{qq}^{(2)}&m_{ss}^{(2)}&c^{(2)}\\\ m_{qq}^{(3)}&m_{ss}^{(3)}&c^{(3)}\\\
\end{bmatrix}\begin{bmatrix}qq_{j,i}\\\ ss_{j,i}\\\ 1\\\
\end{bmatrix}\leq\begin{bmatrix}0\\\ 0\\\ 0\\\ \end{bmatrix}$ (30)

###  2.3 MILP formulation

MILP formulation of the pump scheduling problem is composed of (1) objective
function representing the total pumping cost, (2) a set of equality
constraints representing the originally linear and linearized network
component equations and auxiliary relationships, (3) a set of inequality
constraints representing binary linearized network component equations and
additional constraints such as e.g. symmetry-breaking constraints, (4) a set
of lower bound (LB) and upper bound (UB), aka. box constraints on selected
decision variables, and (5) a vector of indices of binary decision variables.

####  2.3.1 Objective

The objective is the total cost of pumping over time horizon $K$, given the
energy tariff $T(k)$ and the linearized pumping cost model for each pump
$P_{j}(k)$. $\Delta t$ is the time-step - usually 1h.

$TC=\sum_{k=1}^{K}\sum_{j\in E_{pump}}P_{j}(k)\;T(k)\;\Delta t$ (31)

Note that it is also common for the objective function to include a term
penalizing pump switching, e.g. (Lansey and Awumah, 1994). Since this term
includes a sum of absolute (or squared) differences between consecutive pump
statuses, it will need to be linearized via introduction of additional
variables and constraints (Shanno and Weil, 1971). In order to not
overcomplicate the current study, inclusion pump switching cost in the
objective is left for later.

####  2.3.2 Additional constraints

##### Symmetry Breaking:

Symmetries arise in MILP problems when the same feasible solution can be
represented in more than one way. These symmetries can lead to redundant
computations as they increase the search space and require the branch & bound
algorithms to explore and compare multiple equivalent branches, slowing down
the optimization process. In pump scheduling, symmetries will arise if
parallel pumps within one pumping station have the same characteristic. These
symmetries are removed by introducing an additional set of inequality
constraints which enforces the priority of pumps, as described in Eq. 32.
Consequently, the lower priority pumps can be switched ON iff the higher
priority pumps are also switched ON, thus preventing the optimizer from
needlessly exploring equivalent solutions with different permutations of
ON/OFF statuses among equal pump units.

$-n_{j+1}(k)+n_{j}(k)\leq 0\quad\forall
j\in\\{1\ldots(n_{pumps}-1)\\}\quad\textrm{for every pump group with equal
pumps}$ (32)

Adding this constraint reduces the search space for each pumping station with
$n_{pumps}$ from $2^{n_{pumps}}$ to $n_{pumps}+1$ (Gleixner et al., 2012).

##### Enforcing tank levels:

To prevent the optimizer from emptying the reservoirs as it tries to reduce
the total pumping cost, the tank level difference between the final time $N$
and the initial time $1$ is bound to a small threshold $\delta_{h_{t}}$ (Menke
et al., 2016).

$h_{t,j}(N)-h_{t,j}(1)\leq\delta_{h_{t,j}}\quad\forall j=1,\,\ldots\,n_{t}$
(33)

The summary of equality and inequality constraints required for the
formulation of our pump scheduling problem are listed in Table 1 and Table 2,
respectively.

Table 1: Equality constraints | Name | Equation(s) | No. of constraints
---|---|---|---
1 | Mass balance in nodes | 4 | $n_{n}\times K$
2 | Head-flow relationship in tanks |  9 \+ 10 | $n_{t}\times K$
3 | Pipe segment flows | 21 | $n_{p}\times K$
4 | Pipe segment selection variables | 22 | $n_{p}\times K$
5 | Linearized pipe headlosses | 23 | $n_{p}\times K$
6 | Pump segment speeds | 24 | $n_{pump}\times K$
7 | Pump segment flows | 25 | $n_{pump}\times K$
8 | Pump segment selection variables | 26 | $n_{pump}\times K$
Table 2: Inequality constraints | Name | Equation | No. of constraints
---|---|---|---
1 | Pump power binary linearization with respect to pump status | 16 | $2\times n_{pump}\times K$
2 | ‘Zero power’ enforcement for switched OFF pumps | 17 | $2\times n_{pump}\times K$
3 | Binary linearization of pipe flow with respect to pipe segment selection | 20 | $3\times n_{p}\times K$
4 | Binary linearization of VSP speed with respect to pump segment selection | 27 | $4\times n_{pump}\times K$
5 | Binary linearization of VSP flow with respect to pump segment selection | 28 | $4\times n_{pump}\times K$
6 | Binary linearized VSP characteristic | 29 | $2\times n_{pump}\times K$
7 |  VSP characteristic domain definitions | 30 | $12\times n_{pump}\times K$
8 | Symmetry breaking in pump groups with equal pumps | 32 | $n_{groups}\,\left(n_{eq.pumps}^{group}-1\right)\times K$
9 | Enforcing final tank level | 33 | $n_{t}$

####  2.3.3 Lower and upper bounds on decision variables

Most of the decision variables in vector $\mathbf{x}$ are rather tightly
constrained by the inequality and equality constraints. The exceptions are:
(a) heads in tanks, which need to be additionally constrained with box-
constraints such that the levels do not violate restrictions imposed by the
tanks’ minimum ($h^{j}_{t,min}$) and maximum ($h^{j}_{t,max}$) levels and (b)
integer variables that we restrict to take only binary values. Tank level
constraints for $j\in\\{1,\ldots,n_{t}\\}$ tanks, are listed below

$h^{j}_{t,min}(k)\leq h_{t}^{j}(k)\leq h^{j}_{t,max}(k)$ (34)

The binary variable constraints can be expressed as follows

$0\leq x_{i}\leq 1\quad\forall_{i}\in\varUpsilon$ (35)

where $\varUpsilon$ is the set of indices of integer decision variables - see
Eq. 1

##  3 Case study

Our method was tested on a model of a simple system illustrated in Fig. 4. The
network is composed of one fixed-head reservoir, one variable-head tank, two
equal parallel VSPs, 4 pipes and one fixed demand node. Its purpose is to show
the correctness of our method on a simple enough network for which the results
are easy to interpret and visualise.

Figure 4: Schematic of a simple network with one tank, one demand point and a
single group of two equal VSPs in parallel.

The network was used in two separate analyses. In the first analysis, the pump
schedules were optimized for a 24 hour time horizon for a single default set
of inputs and parameters: tank elevation $z_{t}=230$ m, final tank level
$x_{t}^{end}$ = 2.50 m = initial tank level $x_{t}^{init}$, average demand
$\bar{d}$ = 42.7 L/s, tank diameter $D_{t}$ = 15.00 m. The aim of the analysis
was to study the behaviour of the MILP solver on our problem formulation and
to verify the correctness of the obtained results. In the second analysis, a
batch of pump schedule optimization was performed for 81 combinations of
network parameters and inputs. The ranges of parameters were as follows: tank
elevations $z_{t}\in[225,230,235]$ m, tank level differences
$x_{t}^{end}-x_{t}^{init}\in[-0.5\,\mathrm{m},0.0\,\mathrm{m},+0.5\,\mathrm{m}]$,
average demands $\bar{d}\in[34.16,42.70,47.00]$ L/s, tank diameters
$D_{t}\in[12.75,15.00,17.25]$ m. The goal was to test the reliability and the
robustness of the method under a range of operating points and to measure the
calculation times required by CPLEX solver to find optimal solutions. In both
studies, the MILP solver was set to terminate upon achieving the MIP gap of
0.05.

###  3.1 Results

Results of the initial analysis are shown in Figures 5, 7, 6 and 8. Results of
the initial simulation are shown on the left, the outputs of the MILP solver
are shown in the middle, whereas the results of the final simulation are shown
on the right. As demonstrated in Fig. 5, the MILP scheduler found an
alternative pump schedule to the initial one. The new schedule reduces the
pumps’ energy consumption in high tariff periods - see Fig. 6. Consequently,
the operating cost per day of the network reduced from the initial cost of
70.2 GBP to 64.7 GBP (from final simulation), i.e. a 7.8% reduction. Fig. 5
(middle) illustrates that switching Pump 1 always precedes switching Pump 2
and the speed of inactive pumps is always set to zero. It is a desired
behaviour enforced by the symmetry breaking constraint and binary
linearization of the pump speed, respectively. The schedule produced by the
MILP solver is translated into a schedule supported by the simulator, which
treats multiple pumps as a group of equal pumps operating at equal speeds, not
as separate individual units - see Eqs. 5 and 7.

Figure 5: VSP pump schedules. Figure 6: Pumping energy cost and electricity
tariff.

Flows in the selected network elements and heads in the selected network nodes
are significantly altered by the new pump schedule - see the subplots in the
right and in the middle vs. the left in Fig. 7 and Fig. 8, respectively. As
the pumps’ operation switched from constant to one in which the tank’s storage
capacity is utilized in order to reduce pumping during high tariff periods,
the flows in the elements between the pumps and the tank exhibit more
variability. Consequently, the flows are higher, in absolute values, during
the times when the tank is filling and emptying.

It is interesting to observe the discrepancies in the network state (heads and
flows) between MILP that works with linear approximations of network
components and the simulator that uses a complete (nonlinear) network model.
We can notice an offset between the head at the demand node $h_{6}$ returned
by the optimizer (middle) and the simulator (right), and relatively higher
pump outlet heads $h_{3}$ returned by the optimizer compared to the simulator.
These differences stem from the inaccuracies introduced by the linear and
piece-linear approximations The approximations can be tightened via
introduction of a larger number of piece-linear segments or by iterative
adjustment of the locations of the break-points. Both approaches normally come
at the cost of increasing optimization times. The similarity in tank levels is
preserved in MILP to a greater degree than the heads in non-storage nodes.
This is due to the fact that tank levels change in response to changes in
flows which are affected by approximations less than pressures as the former
are forced inputs in demand driven simulations and the latter are the outputs
and therefore, are determined by the model.

Figure 7: Flows through the pump group, tank feed pipe and demand supply node.
Figure 8: Heads at the pump inlet and outlet, node $n_{c}^{(3)}$ (node 4),
demand node and the tank.

Fig. 10 shows the distribution of optimization times out of 80 successful
CPLEX runs that produced optimal solutions within MIP gap of 0.05. All
calculations were performed on a laptop equipped with Intel(R) Core(TM)
i5-6300U CPU @ 2.40GHz processor and 16GB of RAM and IBM(R) ILOG(R) CPLEX(R)
v. 22.1.1.0 MILP solver. The results indicate that the optimization times are
bounded (on this problem size) with 75% of calculations returning optimal
solution within $\pm 20\%$ of the median time of 1.0 seconds. Fig. 10 provides
a measure of accuracy of the optimal solution against the simulation. This
metric is quantified as mean average error (MAE) between the reservoir level
time-series from the simulator and the optimizer. As demonstrated, 96.6% of
MILP outputs return results in which the metric is between 0,1 and 0.3 metres
- a rather robust outcome.

Figure 9: Histogram of optimization times returned by CPLEX.

Figure 10: Histogram of the mean absolute error (MAE) between the simulated
tank level and the tank level returned from MILP.

The results of the 53 successful optimizations out of 54 attempts (the other
27 runs for zero tank level difference are not shown) are visualised in Fig.
12 and Fig. 12. The purpose of these visualisations is to demonstrate that the
results obtained from the MILP solver are physically correct and smooth, which
means that the optimizer is able to find a global solution in a robust way for
a range of network parameters and inputs. As suspected, larger operating costs
are incurred when the final tank level needs to be 0.5m higher than initial
(see Fig. 12), contrary to Fig. 12 where the opposite is true. Operating costs
also increase with demand, due higher pumped volumes and larger headlosses,
and when the tank is positioned at higher elevations and the pumps need to
overcome larger head differences.

Figure 11: Optimized pumping costs for a range of tank diameters, tank
elevations and water demands, at final tank level difference of +0.5m, i.e.
final level 0.5m above the initial level.

Figure 12: Optimized pumping costs for a range of tank diameters, tank
elevations and water demands, at final tank level difference of -0.5m, i.e.
final level 0.5m below the initial level.

##  4 Conclusions and further work

The study demonstrates that mixed integer linear programming with linear and
piece-linear approximations of the objective and of the model components is a
valid method for finding globally optimal pump schedules in networks
containing variable-speed pumps (VSPs). Although the method was tested on a
very small network, the average calculation time of approx. 1 second is a
promising result indicating that the same approach can be adapted to solving
more complex networks. The method proved to be robust and able to arrive at
global optimal solution within similar calculation times for a range of
operating points. As the solvers for mixed integer programming problems, such
as CPLEX, GUROBI, MOSEK, etc. have become faster and now support parallel
execution, it is perhaps a good idea to start reintroducing mixed integer
linear programming techniques, that are known for their stability and
robustness, into WDNs operation studies. A particularly suited application
would be real-time pump optimization. Another application could be a
hierarchical two-level optimization in which the inner pump schedule
optimization loop is solved using mixed integer linear programming whilst
higher level decisions e.g. long-term policies, design options, etc. are
optimized using evolutionary algorithms. The work presented here shows just
one out of many ways of formulating the problem. It is most likely not the
optimal method nor a complete one as it does not support some of the important
network elements such as e.g. PRVs, CVs, or other aspects such as e.g.
handling pressure-dependent demands. Improvement of solution accuracy and
speed can be attempted by experimenting with different component approximation
and relaxation techniques, different approximation accuracy improvements and
relaxation bound tightening, or different piece-linear approximation
representations such as e.g. special ordered sets. These techniques can be
borrowed from the existing literature or developed new. The scalability and
the speed of the method can be improved via decomposition techniques such as
e.g. Lagrangian or Benders decomposition. We are currently developing a free
open-source Python package that summarizes the state-of-the art in pump
scheduling using mixed integer linear programming and can be used for solving
practical problems on EPANET networks using different MILP formulations and
for adding new features and enhancements. The software (under development) is
currently hosted in the `dev-python` branch of the GitHub repository of
_MILOPS-WDN - the Mixed Integer Linear Optimal Pump Scheduler_ (Janus and
Ulanicki, 2023).

## References

  * Bonvin et al. [2017] Gratien Bonvin, Sophie Demassey, Claude Le Pape, Nadia Maïzi, Vincent Mazauric, and Alfredo Samperio.  A convex mathematical program for pump scheduling in a class of branched water networks.  _Applied Energy_ , 185:1702–1711, 2017.  ISSN 0306-2619.  doi:https://doi.org/10.1016/j.apenergy.2015.12.090.  URL https://www.sciencedirect.com/science/article/pii/S0306261915016633.  Clean, Efficient and Affordable Energy for a Sustainable Future.
  * Reed and Hadka [2014] Patrick M. Reed and David Hadka.  Evolving many-objective water management to exploit exascale computing.  _Water Resources Research_ , 50(10):8367–8373, 2014.  doi:https://doi.org/10.1002/2014WR015976.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2014WR015976.
  * Menke et al. [2016] Ruben Menke, Edo Abraham, Panos Parpas, and Ivan Stoianov.  Exploring optimal pump scheduling in water distribution networks with branch and bound methods.  _Water Resources Management_ , 30(14):5333–5349, Nov 2016.  ISSN 1573-1650.  doi:10.1007/s11269-016-1490-8.  URL https://doi.org/10.1007/s11269-016-1490-8.
  * Mala-Jetmarova et al. [2017] Helena Mala-Jetmarova, Nargiz Sultanova, and Dragan Savic.  Lost in optimisation of water distribution systems? a literature review of system operation.  _Environmental Modelling & Software_, 93:209–254, 2017\.  ISSN 1364-8152.  doi:https://doi.org/10.1016/j.envsoft.2017.02.009.  URL https://www.sciencedirect.com/science/article/pii/S1364815216307769.
  * Wu et al. [2018] Wei Wu, Zhiguo Huang, and Shunli Cai.  Optimal pump scheduling in water distribution systems: A review.  _Journal of Water Resources Planning and Management_ , 144(9):04018053, 2018.
  * Cplex, IBM ILOG [2009] Cplex, IBM ILOG.  V12. 1: User’s Manual for CPLEX.  _International Business Machines Corporation_ , 46(53):157, 2009.
  * Gurobi Optimization, LLC [2023] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2023.  URL https://www.gurobi.com.
  * MOSEK [2023] MOSEK.  _MOSEK Optimizer API for Python 10.0.46_ , 2023.  URL https://docs.mosek.com/latest/pythonapi/index.html.
  * Bestuzheva et al. [2021] Ksenia Bestuzheva, Mathieu Besançon, Wei-Kun Chen, Antonia Chmiela, Tim Donkiewicz, Jasper van Doornmalen, Leon Eifler, Oliver Gaul, Gerald Gamrath, Ambros Gleixner, Leona Gottwald, Christoph Graczyk, Katrin Halbig, Alexander Hoen, Christopher Hojny, Rolf van der Hulst, Thorsten Koch, Marco Lübbecke, Stephen J. Maher, Frederic Matter, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Daniel Rehfeldt, Steffan Schlein, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Boro Sofranac, Mark Turner, Stefan Vigerske, Fabian Wegscheider, Philipp Wellner, Dieter Weninger, and Jakob Witzig.  The SCIP Optimization Suite 8.0.  Technical report, Optimization Online, December 2021.  URL http://www.optimization-online.org/DB_HTML/2021/12/8728.html.
  * Anderson and Al-Jamal [1995] Edward J. Anderson and Khairy H. Al-Jamal.  Hydraulic-network simplification.  _Journal of Water Resources Planning and Management_ , 121(3):235–240, 1995.  doi:10.1061/(ASCE)0733-9496(1995)121:3(235).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9496%281995%29121%3A3%28235%29.
  * Deuerlein [2008] Jochen W. Deuerlein.  Decomposition model of a general water supply network graph.  _Journal of Hydraulic Engineering_ , 134(6):822–832, 2008.  doi:10.1061/(ASCE)0733-9429(2008)134:6(822).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9429%282008%29134%3A6%28822%29.
  * Alzamora et al. [2014] F. Martinez Alzamora, B. Ulanicki, and E. Salomons.  Fast and practical method for model reduction of large-scale water-distribution networks.  _Journal of Water Resources Planning and Management_ , 140(4):444–456, 2014.  doi:10.1061/(ASCE)WR.1943-5452.0000333.  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29WR.1943-5452.0000333.
  * Shamir and Salomons [2008] Uri Shamir and Elad Salomons.  Optimal real-time operation of urban water distribution systems using reduced models.  _Journal of Water Resources Planning and Management_ , 134(2):181–185, 2008.  doi:10.1061/(ASCE)0733-9496(2008)134:2(181).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9496%282008%29134%3A2%28181%29.
  * Gleixner et al. [2012] Ambros M. Gleixner, Harald Held, Wei Huang, and Stefan Vigerske.  Towards globally optimal operation of water supply networks.  _Numerical Algebra, Control and Optimization_ , 2(4):695–711, 2012.  ISSN 2155-3289.  doi:10.3934/naco.2012.2.695.  URL /article/id/80157e62-60a5-499b-a676-285b901f9a33.
  * Fooladivanda and Taylor [2018] Dariush Fooladivanda and Joshua A. Taylor.  Energy-optimal pump scheduling and water flow.  _IEEE Transactions on Control of Network Systems_ , 5(3):1016–1026, 2018.  doi:10.1109/TCNS.2017.2670501.
  * Singh and Kekatos [2019] Manish K. Singh and Vassilis Kekatos.  Optimal scheduling of water distribution systems, 2019.
  * Pulido-Calvo and Gutierrez-Estrada [2011] Inmaculada Pulido-Calvo and Juan Carlos Gutierrez-Estrada.  Selection and operation of pumping stations of water distribution systems.  _Environmental Research Journal_ , 5(Issue 3):1–20, 2011.  ISSN 1935-3049.  Nova Science Publishers, Inc.
  * Ulanicki et al. [2007] B. Ulanicki, J. Kahler, and H. See.  Dynamic optimization approach for solving an optimal scheduling problem in water distribution systems.  _Journal of Water Resources Planning and Management_ , 133(1):23–32, 2007.  doi:10.1061/(ASCE)0733-9496(2007)133:1(23).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9496%282007%29133%3A1%2823%29.
  * Ghaddar et al. [2015] Bissan Ghaddar, Joe Naoum-Sawaya, Akihiro Kishimoto, Nicole Taheri, and Bradley Eck.  A lagrangian decomposition approach for the pump scheduling problem in water networks.  _European Journal of Operational Research_ , 241(2):490–501, 2015.  ISSN 0377-2217.  doi:https://doi.org/10.1016/j.ejor.2014.08.033.  URL https://www.sciencedirect.com/science/article/pii/S0377221714007103.
  * Naoum-Sawaya et al. [2015] Joe Naoum-Sawaya, Bissan Ghaddar, Ernesto Arandia, and Bradley Eck.  Simulation-optimization approaches for water pump scheduling and pipe replacement problems.  _European Journal of Operational Research_ , 246(1):293–306, 2015.  ISSN 0377-2217.  doi:https://doi.org/10.1016/j.ejor.2015.04.028.  URL https://www.sciencedirect.com/science/article/pii/S0377221715003215.
  * Vieira et al. [2020] Bruno S. Vieira, Sérgio F. Mayerle, Lucila M.S. Campos, and Leandro C. Coelho.  Optimizing drinking water distribution system operations.  _European Journal of Operational Research_ , 280(3):1035–1050, 2020.  ISSN 0377-2217.  doi:https://doi.org/10.1016/j.ejor.2019.07.060.  URL https://www.sciencedirect.com/science/article/pii/S0377221719306393.
  * Liu et al. [2020] Yang Liu, Clayton Barrows, Jordan Macknick, and Meagan Mauter.  Optimization framework to assess the demand response capacity of a water distribution system.  _Journal of Water Resources Planning and Management_ , 146(8):04020063, 2020.  doi:10.1061/(ASCE)WR.1943-5452.0001258.  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29WR.1943-5452.0001258.
  * Salomons and Housh [2020] Elad Salomons and Mashor Housh.  A practical optimization scheme for real-time operation of water distribution systems.  _Journal of Water Resources Planning and Management_ , 146(4):04020016, 2020.  doi:10.1061/(ASCE)WR.1943-5452.0001188.  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29WR.1943-5452.0001188.
  * Bonvin et al. [2021] Gratien Bonvin, Sophie Demassey, and Andrea Lodi.  Pump scheduling in drinking water distribution networks with an lp/nlp-based branch and bound.  _Optimization and Engineering_ , 22(3):1275–1313, Sep 2021.  ISSN 1573-2924.  doi:10.1007/s11081-020-09575-y.  URL https://doi.org/10.1007/s11081-020-09575-y.
  * Tasseff et al. [2022] Byron Tasseff, Russell Bent, Carleton Coffrin, Clayton Barrows, Devon Sigler, Jonathan Stickel, Ahmed S. Zamzam, Yang Liu, and Pascal Van Hentenryck.  Polyhedral relaxations for optimal pump scheduling of potable water distribution networks, 2022.
  * Tasseff et al. [2019] Byron Tasseff, Russell Bent, and Carleton Coffrin.  WaterModels.jl.  https://github.com/lanl-ansi/WaterModels.jl, 2019.  Accessed: 2023-07-27.
  * Janus and Ulanicki [2023] Tomasz Janus and Bogumil Ulanicki.  MILOPS-WDN: Mixed Integer Linear Optimal Pump Scheduler for Water Distribution Networks, 2023.  URL https://github.com/tomjanus/milp-scheduling.
  * Strafaci et al. [2007] A. Strafaci, D.V. Chase, T.M. Walski, S. Beckwith, W. Grayman, D. Savic, Inc Haestad Methods, and E. Koelle.  _Advanced Water Distribution Modeling and Management_.  Haestad methods water resources modeling collection. Bentley Institute Press, 2007.  ISBN 9781934493014.  URL https://books.google.co.ma/books?id=N9k-PgAACAAJ.
  * Ulanicki et al. [2008] Bogumil Ulanicki, Jens Kahler, and Brian Coulbeck.  Modeling the efficiency and power characteristics of a pump group.  _Journal of Water Resources Planning and Management_ , 134(1):88–93, 2008.  doi:10.1061/(ASCE)0733-9496(2008)134:1(88).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9496%282008%29134%3A1%2888%29.
  * Sárbu and Borza [1998] Ioan Sárbu and Ioan Borza.  ENERGETIC OPTIMIZATION OF WATER PUMPING IN DISTRIBUTION SYSTEMS.  _Periodica Polytechnica Mechanical Engineering_ , 42(2):141–152, 1998.  URL https://pp.bme.hu/me/article/view/5462.
  * Simpson and Marchi [2013] Angus R. Simpson and Angela Marchi.  Evaluating the approximation of the affinity laws and improving the efficiency estimate for variable speed pumps.  _Journal of Hydraulic Engineering_ , 139(12):1314–1317, 2013.  doi:10.1061/(ASCE)HY.1943-7900.0000776.  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29HY.1943-7900.0000776.
  * Janus and Ulanicki [2020] Tomasz Janus and Bogumil Ulanicki.  Calculating power consumption of a number of parallel identical pumps in EPANET.  https://github.com/OpenWaterAnalytics/EPANET/pull/599, 2020.  Accessed: 2023-07-21.
  * Lansey and Awumah [1994] Kevin E. Lansey and Kofi Awumah.  Optimal pump operations considering pump switches.  _Journal of Water Resources Planning and Management_ , 120(1):17–35, 1994.  doi:10.1061/(ASCE)0733-9496(1994)120:1(17).  URL https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9496%281994%29120%3A1%2817%29.
  * Shanno and Weil [1971] David F. Shanno and Roman L. Weil.  Technical note—“linear” programming with absolute-value functionals.  _Oper. Res._ , 19(1):120–124, feb 1971.  ISSN 0030-364X.  doi:10.1287/opre.19.1.120.  URL https://doi.org/10.1287/opre.19.1.120.

.

2D

    two-dimensional
ACO

    ant colony optimization
ANN

    artificial neural network
CV

    check valve
EA

    evolutionary algorithm
EPS

    extended period simulation
FSP

    fixed-speed pump
GA

    genetic algorithm
GPU

    graphical processing unit
HPC

    high performance computing
LB

    lower bound
LP

    linear programming
LPG

    Lagrangian Relaxation and Primal Greedy
MILP

    mixed integer linear program
MINLP

    mixed integer nonlinear program
NRV

    non-return valve
MAE

    mean average error
MOEA

    multiobjective evolutionary algorithm
NLP

    nonlinear programming
OR

    operational research
OWF

    optimal water flow
WDN

    water distribution network
PSO

    particle swarm optimization
PSP

    pump scheduling problem
PRV

    pressure reducing valve
SOC

    second order cone
SOS

    special ordered set
SOS2

    special ordered set type II
UB

    upper bound
VSP

    variable-speed pump

##  5 Appendix

###  5.1 Expanded expressions for $m_{j}^{q}$, $m_{k}^{s}$ and $c_{j}$ in the
linearized pump power equation 15

$m_{j}^{q}=3\,a_{3,j}\,q_{0}^{2}+2\,a_{2,j}\,s_{0}\,q_{0}+a_{1,j}\,s_{0}^{2}$
$m_{j}^{s}=a_{2,j}\,q_{0}^{2}+2\,a_{1,j}\,q_{0}\,s_{0}+3\,a_{0,j}\,s_{0}^{2}$
$c_{j}=-2\left(a_{0,j}\,s_{0}^{3}+a_{1,j}\,q_{0}\,s_{0}^{2}+a_{2,j}\,q_{0}^{2}\,s_{0}+a_{3,j}\,q_{0}^{3}\right)=-2\,P_{j}\,(q_{0},s_{0})$

###  5.2 Derivation of the linearized pump characteristic

The linearized characteristic i.e. the equations of the four planes containing
$A_{1}$, $A_{2}$, $A_{3}$ and $A_{4}$, respectively, are derived from the
vector normal to each plane. The normal vector $\vec{n}$ is calculated as a
vector product of two vectors formed by the sides of each (triangular) plane.
From definition of orthogonality, the dot product of $\vec{n}$ and
$[s,q,H]-p_{n}$ has to be equal to zero if $[s,n,H]$ lies on the same plane as
$p_{n}$. Let’s follow the procedure for the segment $A_{1}$.

$\vec{n_{A_{1}}}=\vec{L_{1}}\times\vec{L_{2}}=\begin{vmatrix}\vec{i}&\vec{j}&\vec{k}\\\
s_{p_{n}}-s_{p_{1}}&q_{p_{n}}-q_{p_{1}}&H_{p_{n}}-H_{p_{1}}\\\
s_{p_{n}}-s_{p_{2}}&q_{p_{n}}-q_{p_{2}}&H_{p_{n}}-H_{p_{2}}\end{vmatrix}$ (36)
$\vec{n_{A_{1}}}\cdot\left([s,q,H]-[s_{n},q_{n},H_{n}]\right)=0$ (37)

Solving the above equation yields the linear plane equation
$dd_{j,i}\,s_{j,i}+ee_{j,i}\,q_{j,i}+ff_{j,i}\,H_{j,i}=gg_{j,i}$ for each of
the four planes, in which $dd_{j,i}$, $ee_{j,i}$, $ff_{j,i}$ and $gg_{j,i}$
are the calculated coefficients. The procedure is repeated for each of the
four segments of the linearized pump characteristic and for each pump $j$.

  *[WDN]: water distribution network
  *[VSPs]: variable-speed pump
  *[MILP]: mixed integer linear program
  *[WDNs]: water distribution network
  *[PSO]: particle swarm optimization
  *[ACO]: ant colony optimization
  *[GAs]: genetic algorithm
  *[MINLP]: mixed integer nonlinear program
  *[OR]: operational research
  *[SOC]: second order cone
  *[MINLPs]: mixed integer nonlinear program
  *[PRVs]: pressure reducing valve
  *[LB]: lower bound
  *[UB]: upper bound
  *[VSP]: variable-speed pump
  *[MAE]: mean average error
  *[CVs]: check valve