Adaptive stochastic management of the storage function for a large open reservoir using an artificial intelligence method

J. Hydrol. Hydromech., 67, 2019, 4, 314–321 DOI: 10.2478/johh-2019-0021 314 Adaptive stochastic management of the storage function for a large open reservoir using an artificial intelligence method Tomas Kozel*, Milos Stary Brno University of Technology, Faculty of Civil Engineering, Institute of Landscape Water Management, Veveří 331/95, Brno, Czech Republic. *Corresponding author. E-mail: [email protected]tbr.cz Abstract: The design and evaluation of algorithms for adaptive stochastic control of reservoir function of the water reservoir using artificial intelligence methods (learning fuzzy model and neural networks) are described in this article. This procedure was tested on an artificial reservoir. Reservoir parameters have been designed to cause critical disturbances during the control process, and therefore the influences of control algorithms can be demonstrated in the course of controlled outflow of water from the reservoir. The results of the stochastic adaptive models were compared. Further, stochastic model results were compared with a resultant course of management obtained using the method of classical optimisation (differential evolution), which used stochastic forecast data from real series (100% forecast). Finally, the results of the dispatcher graph and adaptive stochastic control were compared. Achieved results of adaptive stochastic management provide inspiration for continuing research in the field. Keywords: Stochastic; Artificial intelligence; Storage function; Optimisation. INTRODUCTION Water shortage problems have begun to appear in the whole of the Czech Republic (Crhová et al., 2019). Water shortage manifests from the prolonging and deepening of dry seasons. This reality leads to growing tension between capacity of water resources and water user demand. Many problems with water demand occurred in 2015, 2016 and 2018. Long-term average flow Qa is decreasing due to lower values of flow in river networks and prolonging of the dry season, because the low flow rates will not be able to sufficiently dilute the pollutants entering into them. Assumed future values could decrease to 0.8 Qa or even lower values (Kašpárek, 2005). This decrease will not only impact water supply, but it could even influence quality of water. Some rivers can be transformed to drains in the worst scenario, because the low flow rates will not be able to sufficiently dilute the pollutants entering into them. If new large water reservoirs were built, it would lead to an improvement of the current situation. But construction of new large water reservoirs is complicated nowadays. Therefore, optimised current management of reservoirs is required. Controlling the outflow of water from the reservoir will be understood in the following text as a strategic management of storage function of the reservoir using monthly time steps. Reservoir control is the management of an isolated reservoir with a single inflow of water and one regulated flow. At the level of average monthly flows, water flows in streams can be considered as random (stochastic) processes (sequences) for which their future values cannot be accurately determined (Hirsh, 1979; Svanidze, 1961). Probability of their future occurrence can be estimated. Therefore, it is a great simplification to predict their values deterministically. However, a deterministic prediction of a random process can be burdened with a major error, and its use for system control can be misleading. This problem should be approached stochastically, at least to quantify the range of their possible occurrence with a certain probability distribution. A range of possible occurrences increases with increasing prediction length. It follows from the previous paragraph that, when using only one value (deterministic control), there is a significant simplification of the problem, loss of management accuracy or mistaken assessment of the situation (values of real inflow to reservoir may differ significantly from what was assumed). On the other hand, stochastic control allows us to work with a certain scatter of values of controlled outflows (with a given probability distribution). Proper risk assessment and grasping of the options offered by the approach is able to significantly reduce failure risk of the managed reservoir’s supply function. A range of controlled outflows will therefore provide us with a choice of managed outflow based on the probability of overtaking. Therefore, it is desirable to shift from a deterministic control to a stochastic one. METHODS Managing the storage function of the reservoir in an adaptive way (one of the methods of artificial intelligence) allows the problem of controlling the storage function of the reservoir, with a consideration of stochastic flow forecasting, to be very well described. In the transition from deterministic to stochastic control of outflows in adaptive management, use of models based on the Monte Carlo method is offered. The principle of the Monte Carlo method is applied in constructing forecasts that are extrapolations of historically measured flow lines to which a random component is repeatedly added. In this way, a fan of possible future combinations of water inflows into the reservoir is created. Use of optimisation algorithms to control outflow of water from a reservoir brings the problem of computer performance limitations (high machine time requirements for calculation). For a reasonable evaluation of repeated random states, at least 300 repetitions of the calculation are required. Stochastic adaptive control will be understood in the following text as a control of storage function of reservoir, in which future controlled outflows are calculated from actual value of storage volume and the spectrum of random generated inflows of given length (number of members of the series of forecasted average flows). Therefore, for each forecast the Adaptive stochastic management of the storage function for a large open reservoir using an artificial intelligence method 315 Fig. 1. Schema of main algorithm. optimal course of controlled outflow is repeatedly sought. Controlled outflow values are then processed into an exceedance curve, and a controlled outflow value with a specified probability of exceedance is used for control. Simplified schema of the main process is in Figure 1. If the reservoir fill volume values are considered based on actual measured values and inflow forecasts based on the measured data, control algorithms described in the article can be used for operational control of the reservoir's storage function in real time. Here, the Differential Evolution method DE as an optimization method for controlling was chosen (Price et al., 2005; Storn and Price, 1997), and the DE method search area of acceptable solutions Ω, with a number of dimensions equal to number of forecasted members in forecasted inflow line, was used. The DE method quantifies a series of outflows of water from the reservoir for specified initial volume of water in the reservoir and predicted water inflows into the reservoir during the solving period. The optimisation criterion is the sum of differences of second powers between target (improved) average monthly water outflow from the reservoir Op and a series of controlled average monthly outflows of water from reservoir O which is minimized: (1) where OP is value of the target outflow (mean monthly improved outflow), Oj is value of the calculated controlled average monthly outflow, N is number of months, j is number of months and π is value of the critical function. The choice of this criterion implicitly introduces long and shallow disturbances in the malfunctioning fault periods (which is desirable) in an effort to avoid short deep failures that are problematic from the point of view of reservoir management. The disadvantage of this control method using the DE optimisation method is the large machine demands for calculation. It is possible to use Artificial Intelligence (UI) methods, which are capable of replacing the DE method with a certain loss of accuracy. From Artificial Intelligence Methods, a learning fuzzy model of the Mamdani type (Sugeno, 1977; Tagaki and Sugeno, 1985) was chosen and a neural network containing a three-layer perceptron network (Caudill and Butler, 1992) was selected. UI methods are very well described (Donald et al., 1994). Model schemes are shown in Figures 2 and 3. The learned fuzzy model uses the step-by-step aggregation of inputs described (Janál and Starý, 2009; Janál and Starý, 2012), which allows easy creation of a matrix of rules for fuzzy models with more than two inputs. For this purpose, the DE method was used in advance to construct a matrix of target behaviour patterns of the controlled reservoir (training matrix). It then served to train the learning system, i.e. to find a generalised relationship between the status of the reservoir and the used controlled outflow. The state of the reservoir is described by its immediate filling and the predicted monthly inflow vector. When using the basic equation of the reservoir in differential shape, the relation between the average monthly inflows Qτ and the average monthly outflows Oτ and the volume of the water in the tank at the beginning of the time step Vτ-1 and the volume at the end of the time step Vτ is given in the time step τ with duration Δt by relationship (2) 1 VV QO t ττ ττ − − −= Δ (2) For the time step τ = 1, V0 is the initial condition. The members of the series O τ for τ = 1, 2, ..., N can acquire infinitely many values that depend on the filling reservoir and how to control outflow of water from the reservoir. () 2 1 , N Pj j OO MIN =   π= − →    Tomas Kozel, Milos Stary 316 Fig. 2. Schema of learned fuzzy model SA for forecast with length N. Fig. 3. Schema of trained NS model for forecast with length N. For the control algorithm alone, Q series (boundary conditions) is replaced by a series of forecast Qp (inflow forecasting vector) for each solution of task in relation (2) at each step and adaptive recurrent control is performed for all time steps τ. The initial volume of water in reservoir V0 is replaced by volume obtained in the real-time operation. For simulation of control, V0 for repeated calculation is replaced by V1 of the previous calculation. If the value of τ > 1 is the initial condition, Vτ-1 for each additional time step τ is calculated according to Equation (2) in which the predicted value is used instead of real inflow value, which is considered to be a real value for calculation, and the value of the outflow is calculated from the previous time step. Before starting the next calculation step, it is necessary to calculate actual value of the storage volume according to Equation (2). The outflow of the water from the reservoir, which is controlled to Op (target outflow) value, can take values from the interval (0, Op >. If the capacity of the reservoir volume is unable to absorb excess water, the outflow of water from reservoir may be higher than Op. If the controlled outflow values discharged from the reservoir are lower than Op, a fault occurs. The aim of the algorithm is to perform a control in which the value of the π criterion is minimised according to Equation (1). Forecasts are generated by the hybrid zonal model described (Kozel and Starý, 2016). In the following text, the forecasting model is only briefly indicated. The forecasting model is a combination of the linear autoregression model described, for example, by Brockwell and Davis (1991) and the zonal model described by Marton et al., 2015. Average monthly flows of historical series are sorted from the smallest to the largest by month with the last measured flow and divided into a predetermined number of zones. Average monthly flows of the real flow line in the given zone and their subsequent course (the length is determined by length of the flow forecast) form the working area of flows (zones). Other historical flows of the real flow line are not used for forecasts. The modified zonal model differs by applying a linear autoregressive model to the selected zone. Historically, measured flows are converted to the level Z (standard normal distribution), so the first data is deprived of asymmetry by the Box-Cox equation (level Y, assumption of the normalised distribution) (Box and Cox, 1964) and subsequently transformed to level Z by standard transformation relationships between normal and normal standard divisions. The zone is determined for the working month according to the latest measured flow. The correlation matrix, which is the basic input for Yule-Walker equations (Yule, 1927; Walker, 1931), is calculated only from historically measured flows occurring in the assigned zone. For the rest of the data, the model does not have access to the correlation matrix. From the assembled correlation matrix, using the Yule-Walker equations, regression coefficients and Qp predictions are calculated for a given period. When applying a forecasting model, simplistic assumptions have been introduced so that average monthly inflow to the reservoir is considered a random process. During construction of artificial flow generators, the uncertainty of flow measurement in the measurement station profile was neglected. To simplify the task, the uncertainties of input data as well as water losses associated with the reservoir operation were neglected. All of these algorithms have been programmed in the Matlab software environment. Adaptive stochastic management of the storage function for a large open reservoir using an artificial intelligence method 317 APPLICATION Fig. 4. Location of measured profile Bílovice nad Svitavou. Stochastic control was applied to a fictitious water tank designed for this purpose in the measuring profile of Bílovice nad Svitavou. Figure 4 shows the position of the profile. The profile was chosen due to availability of data at the workplace and a long series of real average monthly flows that is very little influenced by management of the large water reservoir located at the top of the river basin. Reservoir parameters and target outflow are designed to cause sufficiently long and deep failures during its management, allowing the methods used to control its storage function to demonstrate their effectiveness. The storage volume was set at 51 811 000 m3 and the target outflow from the reservoir Op at a constant value of 4.25 m3/s. An 89-year long series of average monthly flows (from 1921 to 2009) was used for construction and subsequent validation of the forecasting model of water inflow into the reservoir. The series was made by measuring the profile of Bílovice nad Svitavou, located on the Svitavy river. The series was divided into two parts. The first 75 years were used to calibrate forecasting models (modules) and were used for model validation for the last 14 years. Flows in each month have different probability distributions, so their transformation has been transformed into a single distribution (normalised normal distribution). There were 12 transformational relationships. The transformation itself was described in the previous text. For the construction of the target matrix TM, the period 1981–1995 was chosen and the years 1996–2009 were used as the validation period. The TM construction period was selected in view of increased occurrence of drought periods, which are problematic from the point of view of the storage function of reservoir management. If an entire real series (outside the validation period) is selected for the TM construct, the TM would contain a large amount of data which is not problematic from the point of view of storage function, and would obscure the relationship between inputs and outputs. The TM was created from results of the DE method which used real-time series segment (100% forecast). For each forecasting length (number of forward months to be driven) 1 to 12, TM was compiled. For the training of models based on NS neural networks, the backward propagation method was used and fuzzy C-means methods (Bezdec, 1981) were used to learn the fuzzy model for a pre-assembled matrix of rules. After calibration of both models based on the UI methods (fuzzy model SA, a model containing the NS neural network) validation was made, where progressive slides from the real flow line (100% forecast) were used instead of average monthly flow forecasts. In this way, the maximum achievable effects of management using UI methods were determined. The results for the SA model are shown in Figure 5 (prediction length is 6 months). There is a certain loss of accuracy over the DE method shown in Figure 5 which consists of three graphs. The first contains the average monthly inflows for the validation period. The other two charts show the course of controlled volumes (middle chart) and controlled outflow for the SA model and DE method for the same period. After learning the SA and NS control models, they were applied in the validation period using stochastic predictions of water inflows into the reservoir. As inputs to the models, the volume of water in the reservoir at the start of the solution has always been used for each predicted period and a set of corresponding forecasted vectors (created by the hybrid zonal model). For each predicted period, 1000 prediction vectors were always created. After the model provided control for all predictions, an empirical overflow line was constructed for the first controlled outflow of water from the reservoir by Čegodajev formula (Keylock, 2012). The overflow line was further smoothed out by polynomial using the built-in function of program matlab pchip (Kahaner et al., 1988). Then, the value of the controlled outflow that was applied for the first step was deducted for the chosen probability of the exceedance. The above procedure was applied for given probability of overrun at all time steps of the validation period of successive validation. In the first phase of validation, a probability of exceedance of outflow P was maintained at all times. In order to test sensitivity of the achieved effects arising out of control, depending on probability of exceeding, the value gradually changed. The elected P values were 99.9, 99, 95, 90, ..., 5, 1. The EA model was the first to be applied. The results show that the model is not capable of managing a malfunction and generates a very short but deep fault. That is why correction was introduced in the application of controlled outflows. The correction consisted in averaging calculated controlled flow value with its previous realisations (idea taken from the conjugate gradient method). After the correction was introduced, the results improved significantly. Tomas Kozel, Milos Stary 318 Fig. 5. Results of chosen models (for 100% forecast). Fig. 6. Results of chosen models. Additionally, the SA and NS models were applied. Both models again use described correction. After the correction was introduced, the results of both UI models resembled the results provided by the EA model. Selected results of models are plotted in Figure 6; the horizontal axis shows time in months and the vertical axis shows the average monthly discharge control water in the tank in m3/s. The following figures always show the best course for the model achieved for the chosen probability of exceeding P. It can be seen from Figure 6 that all the used control models give a similar course of controlled outflows. Therefore, other criteria were introduced to assess success. Figure 6 shows the results of the individual models and also draws the course obtained by the EA model, which used the real-time series segment as forecast (100% forecast). These results (course) will be referred to as an ideal course in the following text. The criterion E was the sum of the second powers of differences between controlled outflow provided by stochastic control using sto- Adaptive stochastic management of the storage function for a large open reservoir using an artificial intelligence method 319 chastic inflow forecasts and the value of the target outflow. This is the application of the relationship (1) over the entire validation period. The second criterion was the sum of the undelivered water flow rate Er [m3/s]. The resulting values of 1 and 2 are shown in Table 1. In Table 1 the values for Criteria 1 and 2 (E and Er) are given for the individual methods used for control. Table 1 shows that for the criterion 1, better results were achieved by the SA model than the ideal course (EA model with 100% predictions). According to criterion 2, the best result was the ideal management. During the validation of the model, the total number of forecasts used and their influence on control were also tested. If the total forecasts were less than 400, there were significant differences between different control patterns for generated forecasts. For a total number of 500 or more forecasts, individual courses of controlled outflows for repeated generated predictions varied only slightly. Furthermore, in the second phase of validation, the changes of values P (P95 and P60 quantile) during the management process were tested. Selecting P leads to change of strength management (reduction of controlled outflow at a time when enough water is in the reservoir is unnecessary). If the volume of water in the reservoir drops below 0.6 from storage volume, the value P changes from P60 to P95 and vice versa. The combination above illustrates the possibility of using a fan of possible values that provide stochastic control. The resulting process was then compared to the results of the first validation phase. The comparison between the results of the SA model is shown in Figure 7. Furthermore, a comparison was made between the results of the ideal course and the results provided by the stochastic control of the SA model using the combination P and the stochastic prediction. The results are shown in Figure 7. Finally, a comparison of the method of adaptive stochastic control (model SA with correction) with the method of dispatcher graphs, which is Table 1. Values of criteria for chosen results. Model/criterium E [m3/s]2 Er [m3/s] SA P95 23.8 39.8 NS P85 25.6 37.9 EA P95 28.6 41.6 DG 50.5 46.5 EA 100% forecast 24.6 28.6 SA 100% forecast 24.4 40.1 SA combination P95+P60 21.6 33.7 shown in Figure 8, was performed. The dispatcher charts were designed as zonal (five zones) according to the method mentioned (Broža, 1981). DISCUSSION AND CONCLUSION Stochastic adaptive control is an appropriate method for controlling the storage function of a large open reservoir. In particular, its contribution to suppressing the influence of uncertainties in development of future trends of water inflow into the reservoir can be expected. It is to be expected that future climate change and subsequent flow changes will not work well with the commonly used dispatching graph (DG) methods. Future changes occurring in flow lines are not contained in existing historical flow lines from which DG are constructed. On the other hand, the stochastic adaptive control described is able to capture high diversity and variability of future inflows. The strong aspect of the processed algorithms is primarily generalization of input/output relationships contained in the matrix of patterns that implicitly carry artificial intelligence methods (learning fuzzy models, neural networks). Adaptive stochastic control is able to capture future climate change if the adaptive principle is applied to the target behaviour matrix and construction of forecast generation. Fig. 7. Results of chosen models. Tomas Kozel, Milos Stary 320 Fig. 8. Results of chosen models. The initially poor EA results were attributed to the optimization that used the entire volume of water in the reservoir, including forecasted inflow values. This property led to lower values of P for emptying storage volume of the reservoir in longer dry periods. Applying the described correction has made it possible to overcome this deficiency. The correction used has proven to be very effective and desirable. The results of the SA model were even better than the ideal course for criterion E. A better result can be attributed to generalizing capabilities of the learned fuzzy model. From the results obtained (see Table 1 and Figure 6), after introduction of the correction, SA and NS models are able to replace the EA optimization model, with a significant shortening of calculation time. The calculation time EA needed was about 7 hours for one P value. Model SA needed only 15 minutes for the same calculation and the NS model needed about 16 minutes. During stochastic control, the suitability of using parallel computations in clusters (EA model), which are able to significantly speed up often very time-consuming operations, has been proven. When using a cluster (6 PC, AMD Phenom X4 9550, 4 cores), the time required to calculate, from 7 hours to 25 minutes, has been shortened to 45 minutes. The use of the cluster opens up possibilities for using methods that would otherwise be inappropriate for their time-consuming performance. The cluster price may be considerably lower than the cost of a powerful computing centre, as clusters can be built from older PCs. The benefit of stochastic adaptive control over use of DG is excellent in the multi-annual management of water outflow from the reservoir where use of DG is problematic. The method of adaptive stochastic control is suitable for multi-year management and, on the basis of the results obtained, it is able to successfully manage the outflow of the water in the reservoir in short and long drought seasons with suitably adjusted strength of control. In the work, a suitable probability of exceeding outflow (the appropriate quantity) was searched for each model. The results showed that the models provided very good results for P85 to P95. Higher probabilities of exceedance represent higher strength of control used. If higher probabilities of exceedance are used for choosing the controlled outflow, the corresponding value of controlled outflow is given from exceeding curve (Higher value of P corresponding with lower values of controlled outflow). It cannot be forgotten that the results depend on the total number and length of the forecast. For most models, the best results were achieved for a forecasting period of six months and a value of 3 for the correction. At the value of three for the correction, the calculated controlled outflow is averaged with the two previously controlled water outflows that were used to control two previous time steps. When replacing predictions with real inflows from the real flow line (when testing models), a longer series of average monthly flows is of course a benefit. Using stochastic predictions, the extension of a range of predicted inflow lines is beneficial only to a certain length, then increasing forecast length acts counterproductively as the forecast becomes largely inaccurate. Additionally, a suitable number of forecasts vectors was tested. Tests have shown that it is appropriate to use at least 500 forecast vectors so that the course of controlled outflows does not vary significantly with repeated generated forecasts. A higher number of forecasts generated did not lead to a significant increase in the effects of management. The method described above is applicable to any reservoir with a storage function. Some reservoirs in the world must work with a very high state of stress between inflows and controlled outflows. In general, the situation in the Czech Republic is different and the operation of the reservoirs is safe. The storage volume of reservoirs was designed for very high security. In operating anomalies or large changes in climate development and increasing demand for water supply, we can expect an increase and a higher frequency of disturbances and the methods described may be used. In conclusion, the results provided by the adaptive stochastic control of the storage function of the reservoir were sufficiently positive to justify further examination. Adaptive stochastic management of the storage function for a large open reservoir using an artificial intelligence method 321 Acknowledgements. The article was supported by SEDECO - Sediments, ecosystem services and interrelation with floods and droughts in the AT-CZ border region, ATCZ28. REFERENCES Bezdec, J.C., 1981. Pattern Recognition with Fuzzy Objective Function Algorithms. Plenum Press, New York. Box, G.E.P., Cox, D.R., 1964. An analysis of transformations. Journal of the Royal Statistical Society, Series B, 26, 211-252. Brockwell, P.J., Davis, R.A., 1991. Time Series: Theory and Methods. 2nd edition. Springer. Broža, V., 1981. Methodological instructions for water management solutions of reservoirs (Metodické návody k vodohospodářským řešením nádrží). ČVUT v Praze, Praha. Caudill, M., Butler, C., 1992. Understanding Neural Networks: Computer Explorations. Vols. 1 and 2. MIT Press, Cambridge, MA. Crhová, L., Čekal, R., Černá, L., Kimlová, M., Krejčová, K., Šádková, E., Štěpánková, B., Vrabec, M., 2019. Annual report on the hydrometeorological situation in the Czech Republic 2018. Czech Hydrometeorological Institute, Prague, Czech Republic. http://portal.chmi.cz/files/portal/docs/hydro/sucho/Zpravy/ROK_ 2018.pdf Donald, M., Spiegelhalter, D.J., Taylor, C.C., 1994. Machine Learning, Neural and Statistical Classification. NJ, USA. ISBN 0-13-106360-X. Hirsh, R.M., 1979. Synthetic hydrology and water supply reliability. Water Resources Research, 15, 6, 1603–1615. Janál, P., Starý, M., 2009. Fuzzy model for predicting the degree of vulnerability basin flooding from torrential rains. Journal of Hydrology and Hydromechanics, 57, 3, 145–153. (In Czech.) Janál, P., Starý, M., 2012. Fuzzy model used for the prediction of a state of emergency for a river basin in the case of a flash flood - PART 2. Journal of Hydrology and Hydromechanic, 60, 3, 162– 173. Kahaner, D., Moler, C., Stephen, N., 1988. Numerical Methods and Software. Prentice Hall, Upper Saddle River, NJ. Kašpárek, L., 2005. Estimation of the volume of reservoir required to compensate for the drop in runoff due to climate change. VÚV, Prague. (In Czech.) Keylock, C.J., 2012. A resampling method for generating synthetic hydrological time series with preservation of cross-correlative structure and higher-order, Properties. Water Resour. Res., 48, 12. Kozel, T., Starý, M., 2016. Zone stochastic forecasting model for management of large open water reservoir with storage function. In: SGEM Conference Proceedings. STEF92 Technology Ltd., Sofia, Bulgaria. Marton, D., Menšík, P., Starý, M., 2015. Using predictive model for strategic control of multi-reservoir system storage capacity. Procedia Engineering, 119, 994–1002. Price, K., Storn, R., Lampinen, J., 2005. Differential Evolution: A Practical Approach to Global Optimization. Springer-Verlag Berlin. Storn, R., Price, K., 1997. Differential evolution - a simple and efficient heuristic for global optimization ov er continuous spaces. Journal of Global Optimization, 11, 341–359. Sugeno, M., 1977. Fuzzy Measures and Fuzzy Integrals. In: Gupta M.M., Saridis, G.N., Ganies, B.R. (Eds.): Fuzzy Automata and Decision Processes. North-Holand, New York, pp. 89–102. Svanidze, G.G., 1961. Mathematical Modelling of Hydrological Series and Some Problems of Long-Term River Runoff Control, AN Gruz. SSSR, 1961, Vol. 14. p. 189–216. Tagaki, H., Sugeno, M., 1985. Fuzzy identification of systems and its applications to modelling and control. IEEE Transactions on Systems, Man and Cybernetics, SMC-15, 1, 116–132. Yule, G.U., 1927. On a method of investigating periodicities in disturbed series, with special reference to Wolfer's sunspot numbers. Philosophical Transactions of the Royal Society of London, Series A, 226, 267–298. Walker, G., 1931. On periodicity in series of related terms. Proceedings of the Royal Society of London, Series A, 131, 518– 532. Received 9 November 2018 Accepted 9 September 2019