Uncertainty-aware demand management of water distribution networks in deregulated energy markets

TL;DR

This paper introduces an open-source, GPU-accelerated stochastic model predictive control approach for managing water distribution networks under demand and electricity price uncertainties in deregulated energy markets, ensuring economic and safe operation.

Contribution

It develops a novel, parallelized scenario-based control method that effectively handles uncertainty in water demand and electricity prices, improving operational reliability and cost-efficiency.

Findings

GPU acceleration significantly reduces computation time.

The approach effectively manages demand and price uncertainties.

Ensures reliable and economical water supply operation.

Abstract

We present an open-source solution for the operational control of drinking water distribution networks which accounts for the inherent uncertainty in water demand and electricity prices in the day-ahead market of a volatile deregulated economy. As increasingly more energy markets adopt this trading scheme, the operation of drinking water networks requires uncertainty-aware control approaches that mitigate the effect of volatility and result in an economic and safe operation of the network that meets the consumers' need for uninterrupted water supply. We propose the use of scenario-based stochastic model predictive control: an advanced control methodology which comes at a considerable computation cost which is overcome by harnessing the parallelization capabilities of graphics processing units (GPUs) and using a massively parallelizable algorithm based on the accelerated proximal gradient method.

Figures

Tables

Table 1. Table 1 : Description of the major class and the associated JSON files in RapidNet .

Class name	Description	JSON file
SmpcController	Runs the accelerated dual proximal gradient algorithm can computes a control action to be applied to the water network. Computations are carried out on GPU and the output (flow set-points) can be stored in a JSON file.	controlOutput.json
DwnNetwork	Encapsulates all information related to the topology, dynamics and constraints of the water network.	network.json
ScenarioTree	Scenario-tree representation of the uncertainty in electricity prices and water demands.	scenarioTree.json
Forecaster	An abstract forecaster which predicts the upcoming electricity prices and water demands using some predictive model (implemented by subclassing Forecaster) or reads the forecasts from a JSON file (so that the user can use forecasts from third-party software).	forecaster.json
Engine	Provides essential functionality to SmpcController and manages the GPU-side memory.
SmpcConfiguration	Contains configuration parameters that are relevant for SmpcController (tuning parameters, solver tolerance, maximum number of iterations).	controllerconfig.json

Table 2. Table 2 : Key methods in RapidNet .

Class	Method	Description
SmpcController	controlAction	Computes control action using the accelerated proximal gradient algorithm to solve the scenario-based stochastic optimal control problem.
Forecaster	predictDemand	Returns water demand forecasts.
	predictPrice	Returns energy price forecasts.
Engine	factorStep	Precomputes certain quantities that facilitate and accelerate the computation of the dual gradient in the algorithm.

Equations

$$\frac{\text{d}{V}^j(t)}{\text{d}t}=\sum _{i\in {\mathcal{I}}_j}{q}^i(t)-\sum _{i\in {\mathcal{J}}_j}{q}^i(t).$$

$$V_{\min }^j\le V^j(t)\le V_{\max }^j.$$

$$\sum _{i\in {\mathcal{K}}_s}{q}^i(t)=\sum _{i\in {\mathcal{C}}_s}{q}^i(t),$$

$$q^i(t)=u^i(t),i\in \mathcal{P}.$$

$$q^i(t)=u^i(t),i\in \mathcal{V}.$$

$$0\le q^i(t)\le q_{\max }^i,i\in \mathcal{P}.$$

$$q^i(t)=d^i(t),i\in \mathcal{D}.$$

$$x_{k+1}$$

$$0$$

$$x_{\min }\le$$

$$u_{\min }\le$$

$$\left[\begin{array}{c}{d}_{k+j}\\ {\alpha }_{k+j}\end{array}\right]=\left[\begin{array}{c}{\hat{d}}_{k+j\mid k}\\ {\hat{\alpha }}_{k+j\mid k}\end{array}\right]+{ϵ}_{j\mid k},$$

$${\text{boldmathd^}}_k=({\hat{d}}_{k+1\mid k},{\hat{\alpha }}_{k+1\mid k},\dots ,{\hat{d}}_{k+H_p\mid k},{\hat{\alpha }}_{k+H_p\mid k}),$$

$$\left[\begin{array}{c}{d}_{k+j}^i\\ {\alpha }_{k+j}^i\end{array}\right]=\left[\begin{array}{c}{\hat{d}}_{k+j\mid k}\\ {\hat{\alpha }}_{k+j\mid k}\end{array}\right]+{ϵ}_{j\mid k}^i,$$

$$E{u}_{k+j\mid k}^i+E_d{d}_{k+j\mid k}^i=0,$$

$$x_{k+j+1\mid k}^l=A{x}_{k+j\mid k}^i+B{u}_{k+j\mid k}^l+G_d{d}_{k+j\mid k}^l,$$

$${\ell }^w(u_k,k)=W_{\alpha }({\alpha }_0+{\alpha }_k{)}^{\prime }{u}_k,$$

$${\ell }^{\Delta }(\Delta u_k)=\Delta u_k^{\prime }{W}_u\Delta u_k,$$

$$\ell (u_k,u_{k-1},k)={\ell }^w(u_k,k)+{\ell }^{\Delta }(\Delta u_k).$$

$${\ell }^S(x_k)=W_s\Vert \max \{0,x_s-x_k\}\Vert ,$$

$${\ell }^x(x_k)$$

$$\displaystyle\qquad\ +\|\max\{0,x_{k}-x_{\max}\}\|\big{)},$$

$$\underset{\text{boldmathx}}{\mathrm{minimise}}\text{E}\left[V(\text{boldmathx},p,q,k)\right]+V_s(\text{boldmathx}),$$

$$x_{k\mid k}=p,u_{k-1\mid k}=q,$$

$$x_{k+j+1\mid k}^l=A{x}_{k+j\mid k}^i+B{u}_{k+j\mid k}^l+G_d{d}_{k+j\mid k}^l,$$

$$E{u}_{k+j\mid k}^i+E_d{d}_{k+j\mid k}^i({ϵ}_j)=0,$$

$$u_{\min }\le u_{k+j\mid k}\le u_{\max }^i.$$

$$V(\text{boldmathx},p,q,k)=\sum _{j=0}^{H_p-1}\ell (u_{k+j\mid k},u_{k+j-1\mid k},k+j),$$

$$V_s(\text{boldmathx})=\sum _{j=0}^{H_p}\sum _{i=1}^{{\mu }_{j-1}}{\ell }^x(x_{k+j\mid k}^i)+{\ell }^S(x_{k+j\mid k}^i).$$

$$\mathcal{P}:\underset{\text{boldmathx}\in {\Re }^n}{\mathrm{minimise}}f(\text{boldmathx})+g(H\text{boldmathx}),$$

$$\delta (x\mid C)=\{\begin{array}{ll}0,& \text{ if }x\in C\\ +\infty ,& \text{ otherwise}\end{array}$$

$$f(\text{boldmathx})$$

$$+\sum _{j=0}^{H_p-1}\sum _{i=1}^{{\mu }_j}{p}_j^i({\ell }^w(u_j^i)+{\ell }^{\Delta }(\Delta u_j^i)),$$

$${\Phi }_1(d)=\{u:Eu+E_{d}d=0\},$$

Full text

††journal: Environmental Modelling & Software

Software availability

Name of software: RapidNet.

Hardware required: CUDA-compliant GPU (tested on NVIDIA Tesla 2075).

Software required: CUDA framework v6.0 or higher (including cuBLAS) and rapidjson (http://rapidjson.org/).

Availability: Open-source software, licence: GNU LGPL v3.0. Available online at https://github.com/GPUEngineering/RapidNet.

Program language: CUDA-C++/C++.

First release: 2017.

1 Introduction

1.1 State of the art

The explosive proliferation of interconnected sensing, computing and communication devices has marked the advent of the concept of cyber-physical systems — ensembles of computational and physical components. In drinking water networks this trend has ushered in new control paradigms where the profusion of data, produced by a network of sensors and stored in a database, is used to prescribe informed control actions [1, 2, 3, 4]. Nevertheless, as these data, be they water demand values or electricity prices, cannot be modeled perfectly, the associated uncertainty is shifted to the decision making process.

The high uncertainty in the operation of drinking water networks, as a result of the volatility of future demands as well as energy prices (in a deregulated energy market) is likely to lead to a rather expensive operating mode with poor quality of service (the network may not always be able to provide the necessary amount of water to the consumers). In control engineering practice, this uncertainty is often addressed in a worst-case fashion [5, 6] — if not neglected at all — leading to conservative and suboptimal control policies. It is evident that it is necessary to devise control methods which take into account the probabilistic nature of the underlying uncertainty making use of the wealth of available historical data aiming at a proactive and foresightful control scheme which leads to an improved closed-loop performance. These requirements necessitate the use of stochastic model predictive control: an advanced control methodology where at every time instant we determine a sequence of control laws which minimizes the expected value of a performance index taken with respect to the distribution of the uncertainty [7]. Optimization-based approaches for the operational management of water networks have been studied and are well established in engineering practice [8].

Indeed, scenario-based stochastic model predictive control (SSMPC) has been shown to lead to remarkable decrease in the operating cost and improvement in the quality of service of drinking water networks [9]. In SSMPC, the uncertain disturbances are treated as random variables on a discrete sample space without assuming any parametric form for their distribution [10]. The scenario approach was identified in a recent review as a powerful method for mitigating uncertainty in environmental modeling related to water management [11]. Although this approach offers a realistic control solution as it is entirely data-driven, this comes with considerable computational burden as the resulting optimization problems are of particularly large scale [12, 13]. This has rendered the use of SSMPC prohibitive and has hindered its applicability. Indeed, hitherto there have been used only conventional model predictive control approaches [14, 15], robust worst-case formulations [5, 16, 17] and stochastic formulations where the underlying uncertainty is assumed to be normally identically independently distributed [18, 19]. Note that it has been observed that demand prediction errors are typically follow heavy-tail distributions which cannot be well approximated by normal ones [20].

In this paper, we present a software for the fast and efficient solution of such problems harnessing the immense computational capabilities of graphics processing units (GPU) building up on our previous work [9, 21, 22].

There has been recently a lot of interest in the development of efficient methods for stochastic optimal control problems such as stochastic gradient methods [23], the alternating directions method of multipliers (ADMM) [24] and various decomposition methods which can lead to parallelizable methods [25, 26] (the most popular being the stochastic dual approximate dynamic programming [27], the progressive hedging approach [28] and dynamic programming [29]). There have been proposed parallelizable interior point algorithms for two-stage stochastic optimal control problems such as [30, 31, 32, 33] and an ad hoc interior point solver for multi-stage problems [34]. However, interior point algorithms involve complex steps and are not suitable for an implementation on GPUs which can make the most of the capabilities of the hardware. Additionally, interior point methods cannot accommodate complex non-quadratic terms in the cost function such as soft constraints (distance-to-set functions).

At large, not many software and libraries are available for stochastic optimal control; one of the very few one may find on the web is JSPD, a generic Java stochastic dynamic programming library. QUASAR is a commercial tool for scenario-based stochastic optimization. One of the most popular tools in the toolbox of the water networks engineer is PLIO [35], which implements MPC algorithms. This work covers the yawning gap between engineering practice and the latest developments in control and optimization theory for drinking water networks. These results can also be applied for the control of other infrastructure with similar structure such as power grids [36].

1.2 Contributions and Novelty

Despite the fact that SSMPC problems typically involve millions of decision variables, the associated optimization problems possess a rich structure which can be exploited to devise parallelizable ad hoc methods to solve the problem more than an order of magnitude faster than commercial solvers running on CPU.

The architecture of our implementation comprises three independent modules: (i) the network module, (ii) the energy prices and water demands forecaster and (iii) the control module. The network module provides a dynamical system model which describes the flow of water across the network together with the storage limits of the tanks and the constraints on pumping capacities. The network module defines a safety storage level for each tank — a level which ensures the availability of water in case of high demand and the maintenance of a minimum required pressure. The forecasters produce a scenario tree, that is, a tree of likely future water demands and energy prices, upon which a contingency plan is made by minimizing a cost function which quantifies the operating cost and the quality of service. Such scenario trees are constructed from historical data of energy prices and water demands. The control module computes flow set-points, which are sent to the pumping stations and valves, by solving a scenario-based stochastic model predictive control problem over a finite prediction horizon.

The proposed stochastic model predictive controller leads to measurable benefits for the operation of the water network. It leads to a more economic operation compared to methods which do not take into consideration the stochastic nature of the energy prices and water demands. In this paper, we assess the performance of the controlled network using three key performance indicators: (i) the economic index, (ii) the safety index, which quantifies the extent of violation of the safety storage level requirement and (iii) the computational complexity index with which we assess the computational feasibility of the controller. Simulation results are provided using data from the water network of Barcelona and the energy market of Austria. The advantages of the adopted control methodology are combined with the computational power of GPUs, which enables us to solve problems of very large scale.

1.3 Software

Our implementation is available as an open-source and free software which can be readily tailored to the needs of different water networks modifying the parameters of its modules. The implementation is done entirely in CUDA-C++ and it can be configured either programmatically or using configuration files. The adopted object-oriented programming model is amenable to extensions and users may specify their own predictive models, scenario trees, cost functions, dynamical models and constraints. Our results are accompanied by extensive benchmarks and the software is verified with unit tests.

2 Modeling

2.1 Hydraulic modeling

The water network involves four types of elements: water tanks, active elements (pumps and valves), mixing nodes and demand sectors. We focus on flow-based water networks where all flows are manipulated variables based on the modeling approaches presented in [37, 14, 38, 13, 5].

Water tanks which play a crucial role in demand management as they ensure water supply, obviate the need for continuous pumping and allow water to be pumped when the price of electricity is lower and provide water in cases of unexpected peaks in demand. Tank dynamics are modelled by simple mass balance equations: let $V^{j}(t)$ be the volume of tank $j=1,\ldots,N_{t}$ at time $t$. We denote all flows in the network by $q^{i}$, $i=1,\ldots,N_{f}$. Let $q^{i}(t)$, $i\in\mathcal{I}_{j}$ be the controlled inflowing streams to tank $j$ and $q^{i}(t)$, $i\in\mathcal{J}_{j}$ be the controlled outflowing streams. Then, the mass balance equation becomes (see Fig. 1(a))

[TABLE]

The volume in each tank should never exceed a maximum limit $V^{j}_{\max}$ and it should always be above a hard lower limit $V^{j}_{\min}$, that is,

[TABLE]

In addition, for the sake of service reliability (availability of water when demand rises unexpectedly) and safety, it is required that the level of water remains above a certain level $V^{j}_{\mathrm{safe}}$. This is allowed to be violated occasionally, when the demand happens to be too high.

Mixing nodes are intersections where flows of water are merged or separated. The mass balance equations for mixing nodes give rise to algebraic constraints of the form

[TABLE]

where $\mathcal{K}_{s}$ and $\mathcal{C}_{s}$ are the sets of indices of the incoming and outgoing flows at node $s=1,\ldots,N_{s}$ as shown in Fig. 1(b).

Pumps and valves are used to control the flow of water in the network and transfer it across tanks and to the demand sectors. We treat these as controlled systems — indeed, pumping stations and valves are equipped with local controllers — to which we prescribe flow set-points. The local control systems operate at a sampling rate of about $1\,\mathrm{Hz}$, while the operational management of the network updates its decisions at a much slower rate (e.g., hourly). It is reasonable to assume that the local control system equilibrates fast enough to neglect its dynamics in the context of operational control. That said, the flow determined by each pump $i\in\mathcal{P}$ is equal to its prescribed set-point $u^{i}$. As shown in Fig. 1(c), that is

[TABLE]

Similarly, as shown in Fig. 1(d) the flow through each valve is

[TABLE]

All flows in the network are unidirectional, so we require that $q^{i}\geq 0$ for all $i=1,\ldots,N_{f}$. Each pump $i\in\mathcal{P}$ has a maximum pumping capacity $q^{i}_{\max}$, that is we require that

[TABLE]

Demand sectors are the exit nodes of the water from the network towards the consumers as shown in Fig. 1(e). At each demand sector $i\in\mathcal{D}$, the mass balance yields

[TABLE]

We treat $d^{i}(t)$ as a random process and elaborate on that in Section 2.2.

We may now describe the system dynamics in terms of the state variable $x(t)=(V^{j}(t))_{j=1,\ldots,N_{t}}$, the input variable $u(t)=(u^{j}(t))_{j\in\mathcal{P}\cup\mathcal{V}}$ and the disturbance $d(t)=(d^{i}(t))_{i\in\mathcal{D}}$.

Let $n_{u}=|\mathcal{P}\cup\mathcal{V}|$ be the dimension of $u(t)$, that is, the total number of pumps and valves. Let $n_{d}=|\mathcal{D}|$ be the total number of demand nodes. By discretizing the dynamical equation (1) using the exact discretization method [39, Sec. 4.2.1] and taking into account the algebraic constraints stated above, we may write the system dynamics in the following form of a discrete-time linear time-invariant system

[TABLE]

where $A\in\Re^{N_{t}\times N_{t}}$, $B\in\Re^{N_{t}\times n_{u}}$, $G_{d}\in\Re^{N_{t}\times n_{d}}$, $E\in\Re^{N_{s}\times n_{u}}$ and $E_{d}\in\Re^{N_{s}\times n_{d}}$.

The constraints on $x_{k}$ and $u_{k}$ can be concisely written as

[TABLE]

with $x_{\min},\,x_{\max}\in\Re^{N_{t}}$ and $u_{\min},\,u_{\max}\in\Re^{n_{u}}$ and $\leq$ is meant in the element-wise sense. These constraints encompass (2) and (6).

2.2 Demands and Electricity Prices

Water demand has been the main source of uncertainty for the operation of drinking water networks and a lot of attention has been paid on the development of models for its prediction. Prediction methodologies span from simple linear models [40] to neural networks [41, 42] and support vector machines [43, 44], nonlinear multiple linear regression [45], Holt-Winters-type models [46], as well as more complex neuro-fuzzy models [47, 48]. Increased predictive ability can be obtained using exogenous information such as weather forecasts [49] and calendar data [5].

In the context of a deregulated wholesale energy market, prices on the day-ahead market are volatile and are often decided on the basis of an auction (bid-based market) among energy companies instead of bilateral agreements with an energy provider. In such cases, energy prices may change on a daily or hourly basis [50]. It is then necessary to be able to predict the day-ahead evolution of the prices using past data; several time series analysis methodologies have been developed for that purpose — see [51] and references therein.

In our approach the prediction procedure is decoupled from the control

which allows the use of any forecasting methodology without having to modify the controller parameters or implementation. An independent forecaster provides estimates of future demands and electricity prices along with an estimation of their uncertainty which is discussed in the next section. At time $k$, the predicted demands for the future time $k+j$ are estimated by a model which computes $\hat{d}_{k+j{}\mid{}k}$. Likewise, we denote the predicted electricity prices by $\hat{\alpha}_{k+j{}\mid{}k}$. Let $d_{k+j}$ and $\alpha_{k+j}$ denote the actual, but unknown at time $k$, values of the water demands and prices. Then

[TABLE]

where $\epsilon_{j{}\mid{}k}$ is a random variable which corresponds to the $j$-step-ahead prediction error at time $k$. At time $k$, a forecaster provides finite-horizon estimates of the upcoming water demands and electricity prices

[TABLE]

where $H_{p}$ is a prediction horizon. This information will then be provided to the controller as we shall discuss in Section 3.

2.3 Uncertainty

It is common in stochastic control-oriented modeling to assume that the errors $\epsilon_{j{}\mid{}k}$ are independently distributed [18, 19]. This assumption however neglects the covariance across the times stages — indeed, if at the future time $j=1$ the model has a large prediction error we would rather expect that the prediction error at time $j=2$ is likely to be large too. This motivates the use of scenario trees: discrete representations of the random processes $(\epsilon_{j{}\mid{}k})_{j}$ which capture such multistage covariances [52].

To date, stochastic modeling for drinking water networks in presence of price uncertainty has received little attention — to the best of the authors’ knowledge, [53] is the only relevant reference — and the scenario tree approach has not been used previously.

A scenario tree is a structure such as the one illustrated in Fig. 2. The scenario tree is organised into time instants $j=0,\ldots,H_{p}$ called stages and a number of nodes at each stage denoted by $\epsilon_{j{}\mid{}k}^{i}$ — these are treated as the possible values of $\epsilon_{j{\mid}k}$. At time $j=0$, the prediction error is always equal to [math] assuming that we observe the current water demands and electricity prices. The corresponding unique node is called the root node of the tree. The nodes of the tree at the last stage are known as leaf nodes. The number of nodes at stage $j$ is denoted by $\mu_{j}$. All nodes are identified by a pair $\alpha=(j,i)$, where $j=0,\ldots,H_{p}$ is the stage index and $i$ is the index of the node in that stage. Each non-leaf node defines a (nonempty) set of children $\operatorname{ch}(j,i)$ which are those nodes in next $j+1$ which are linked to $(j,i)$. Conversely, each not except for the root node defines a unique ancestor denoted by $\operatorname{anc}(j,i)$. The probability of visiting a node $(j,i)$ starting from the root node and following the tree structure is denoted by $p_{j}^{i}$.

Note that the joint demand-price modeling of the uncertainty allows us to cast possible demand-price correlations as in cases of uncertain volume-based pricing.

Using the scenario tree structure, equation (10) yields

[TABLE]

where $j=0,\ldots,H_{p}$ and $i=1,\ldots,\mu_{j}$. Here we see that the tree structure of $\epsilon_{j\mid k}^{i}$ induces a corresponding tree structure upon the water demands and electricity prices, namely $d_{k+j}^{i}$ and $\alpha_{k+j}^{i}$. These are the contingent future water demand and electricity price values associated with the prediction error $\epsilon_{j\mid k}^{i}$.

Similarly, equation (8b) gives

[TABLE]

where $(u_{k+j{}\mid{}k}^{i})_{j,i}$ become the decision variables of the stochastic optimal control problem we shall present in the following section. The discrete-time system dynamics (8a) becomes

[TABLE]

where $(j+1,l)\in\operatorname{ch}(j,i)$.

Scenario trees can be constructed from observed sequences of prediction errors which can be easily obtained in practice using methodologies such as [54] or the popular scenario reduction method [55]. There exist several other scenario generation algorithms such as clustering-based algorithms [56, 57] and simulation and optimization-based approaches [58, 59]. It is not necessary, however, to update the scenario tree for $\epsilon_{j{}\mid{}k}$ at every time instant $k$ — it should be updated occasionally to detect changes in the predictive ability of the forecaster or whenever the predictive model is updated.

3 Scenario-based stochastic optimal control

3.1 Control Objectives

The cost for the operation of the water network is quantified in terms of three individual costs which have been proposed in the literature [5, 60, 16, 61]: the economic cost which is related to the treatment cost and electricity required for pumping, the smooth operating cost which penalizes the abrupt operation of pumps and valves and the safety storage cost which penalizes the use of water from the reserves (i.e., allowing the level in the tanks to drop below the safety level).

The economic cost quantifies the production and distribution cost and it is computed by

[TABLE]

where $\alpha_{0}^{\prime}u_{k}$ is the water production cost (treatment and acquisition fees), $\alpha_{k}^{\prime}u_{k}$ is the uncertain pumping cost and $W_{\alpha}$ is a positive scaling factor.

The smooth operation cost is defined as

[TABLE]

where $\Delta u_{k}=u_{k}-u_{k-1}$ and $W_{u}\in\Re^{n_{u}\times n_{u}}$ is a symmetric positive definite weight matrix.

The total stage cost at a time instant $k$ is the summation of the above costs and is given by

[TABLE]

The safety storage cost penalizes the drop of water level in the tanks below a given safety level $x_{s}$. An elevation above this safety level ensures that there will be enough water in unforeseen cases of unexpectedly high demand and also maintains a minimum pressure for the flow of water in the network. This is given by

[TABLE]

where $W_{s}$ is a positive scaling factor.

The state constraints (9a) should be satisfied at all times without however jeopardizing the feasibility of the optimal control problem we have to solve at every time instant. For that, we introduce an additional cost which penalizes the violation of the state constraints as follows

[TABLE]

where $W_{x}$ is a positive weight factor.

The scaling factors $W_{\alpha}$, $W_{u}$, $W_{s}$ and $W_{x}$ are the tuning knobs of stochastic MPC as we shall discuss in the following section.

3.2 Stochastic optimal control: problem formulation

In scenario-based stochastic MPC, at every time instant $k$ we solve a stochastic optimal control problem which consists in determining an optimal contingency plan for the future course of actions in a causal fashion, that is, our future decisions $u_{k+j{}\mid{}k}$ are only allowed to depend on information that will be available to the controller at time $k+j$ [62]. This was tacitly stated in equation (14). This concept is illustrated in Fig. 3.

We formulate the following scenario-based stochastic MPC problem with a prediction horizon $H_{p}$ and decision variables \mathchoice{\text{\boldmath{\displaystyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}=\{u_{k+j\mid k},x_{k+j+1\mid k}\}_{j=0,\ldots,H_{p}-1} where we minimize the expected total cost along the prediction horizon

[TABLE]

In (19a), $\E$ is the expectation operator, $V$ is the cost function given by

[TABLE]

and $V_{s}$ is the total state constraint violation penalty defined as

[TABLE]

Note that we have a time-varying cost as the electricity prices change with time.

4 Numerical Algorithm

4.1 Problem reformulation

Most modern numerical optimization algorithms such as the (accelerated) proximal gradient algorithm, the alternating directions method of multipliers (ADMM) [63], the Pock-Chambolle method [64], Tseng’s forward-backward-forward algorithm [65] and many another require that the optimization problem be first written in a form

[TABLE]

where $f:\Re^{n}\to\barre\dfn\Re\cup\{+\infty\}$ and $g:\Re^{m}\to\barre$ are convex, lower semi-continuous extended-real-valued functions and $H:\Re^{n}\to\Re^{m}$ is a linear operator. Functions $f$ and $g$ are allowed to return the value $+\infty$ to encode constraints; for example, the constraint $x\in C$ is encoded by the indicator function of the set $C$ which is

[TABLE]

The key question is how to split the optimization problem in (19) so that the resulting formulation is amenable to a fast numerical solution with massive parallelization. For reasons that will be elucidated in Section 4.3 we choose $f$ to be the smooth part of the cost (which corresponds to the linear function $\ell^{w}$ and the quadratic function $\ell^{\Delta}$) plus the indicator of the input-disturbance coupling given in (8b) plus the indicator of the system dynamics in (8a), that is $f:\Re^{n}\to\barre$ is defined as

[TABLE]

where $\Delta u_{j}^{i}=u_{j}^{i}-u_{j-1}^{\operatorname{anc}(j,i)}$ and $\Phi_{1}(d)$ is the affine subspace of $\Re^{n_{u}}$

[TABLE]

and $\Phi_{2}(d)$ is the affine subspace of $\Re^{2n_{x}+n_{u}}$ defined by the system dynamics

[TABLE]

Function $g$ is naturally chosen to be the indicator of the set of input constraints plus the total constraint violation penalty function $V_{s}$. Note, however, that the same variable $x_{k+j{}\mid{}k}^{i}$ participates in both functions $\ell^{x}$ and $\ell^{S}$. As we shall explain in Section 4.3, this complicates any computations thereon. For that reason we introduce a linear operator H:\mathchoice{\text{\boldmath{\displaystyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}\mapsto\mathchoice{\text{\boldmath{\displaystyle y}}}{\text{\boldmath{\textstyle y}}}{\text{\boldmath{\textstyle y}}}{\text{\boldmath{\textstyle y}}}\dfn H(\mathchoice{\text{\boldmath{\displaystyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}) which maps $x_{k+j{}\mid{}k}^{i}$ to $(x_{k+j{}\mid{}k}^{i},x_{k+j{}\mid{}k}^{i})$ and $u_{k+j{}\mid{}k}^{i}$ to itself, that is

[TABLE]

Then, we define the function $g:\Re^{m}\to\barre$

[TABLE]

where $\mathcal{U}\dfn\{u\in\Re^{n_{u}}{}\mid{}u_{\min}\leq u\leq u_{\max}\}$.

The scenario-based optimization problem (19) is now in the form (22) with $f$ and $g$ given by (24) and (4.1) respectively and $H$ given by (27).

4.2 Convex conjugates and proximal operators

Before we can proceed with the statement of the numerical algorithm for the solution of problem (22) we need to introduce a few mathematical notions. A function $f:\Re^{n}\to\barre$ is called proper if it is not everywhere equal to $+\infty$. It is called lower semi-continuous if for every $x\in\Re^{n}$, $\liminf_{z\to x}f(z)=f(x)$. The domain of $f$ is the set $\operatorname{dom}f\dfn\{x\in\Re^{n}{}\mid{}f(x)<\infty\}$. It is called $\kappa$-strongly convex, for some $\kappa>0$, if the function $f(x)-\nicefrac{{\kappa}}{{2}}\|x\|^{2}$ is convex.

For a proper, convex, lower semi-continuous function $f$, we define its convex conjugate to be the convex function $f^{*}:\Re^{n}\to\barre$ defined as [66, Ch. 11]

[TABLE]

An important property is that if $f$ is $\kappa$-strongly convex, then $f^{*}$ is differentiable with $\nicefrac{{1}}{{\kappa}}$-Lipschitz gradient [66, Prop. 12.60] and the gradient of $f^{*}$ is

[TABLE]

By means of the convex conjugate we may derive the Fenchel dual optimization problem of (22) which is [67, Sec. 15.3]

[TABLE]

Fenchel duality offers a more powerful framework as compared to the classical Lagrangian duality approach as it allows us to dualise functions (by means of their convex conjugates) rather than merely constraints.

Under certain conditions on $f$ and $g$ (see Section 4.3), the two problems are equivalent: their optimal values are equal and given a dual-optimal point $y^{\star}$ which is a minimizer of (31), the optimal solution of $\mathcal{P}$ is $x^{\star}=\nabla f^{*}(-H^{\prime}y^{\star})$. This is referred to as strong duality. The reason why we formulate the dual optimization problem $\mathcal{D}$ is because it possesses a favorable structure which can be exploited in the development of fast and parallelizable numerical algorithms (See Section 4.3).

Lastly, for a proper, lower semi-continuous, extended-real valued function $g:\Re^{m}\to\barre$ we define its proximal operator with parameter $\gamma>0$ to be a function $\operatorname{prox}_{\gamma g}:\Re^{m}\to\Re^{m}$ defined as

[TABLE]

Proximal mappings act as generalized projections. For example, the proximal mapping of the indicator function $\delta({\,}{\cdot}{\,}{\mid}{\,}C)$ — cf. (23) — of a nonempty, closed, convex set is the projection onto that set, i.e., $\operatorname{prox}_{\gamma\delta({}{\cdot}{}\mid{}D)}(v)=\operatorname{proj}(v{}{\mid}{}C)$.

The proximal operators of many convex functions (such as Euclidean norm, norm-1, quadratic and linear functions, distance-to-set functions) are easy to compute and typically consist in element-wise operations which can be fully parallelized on a GPU. We shall refer to such functions as prox-friendly [68].

So long as $\operatorname{prox}_{\gamma g}$ is easy to compute, so is $\operatorname{prox}_{\gamma g^{*}}$ and it can be obtained from the Moreau decomposition formula which is

[TABLE]

4.3 Accelerated proximal algorithm on the dual optimization problem

The function f(\mathchoice{\text{\boldmath{\displaystyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}})+g(H\mathchoice{\text{\boldmath{\displaystyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}{\text{\boldmath{\textstyle x}}}), with $f$ and $g$ defined by (24) and (4.1) respectively, is proper, convex and piecewise linear-quadratic on its domain, so, following [66, Thm. 11.42] there is strong duality. As discussed above, since $f$ is strongly convex, $f^{*}$ is differentiable with Lipschitz gradient. Function $g$ is written in the form of a separable sum — a sum of prox-friendly functions of different arguments [63]. Function $g$ is indeed prox-friendly. Let $(\operatorname{prox}_{\gamma g}(v))_{j,i,(s)}$ denote the part of the vector $\operatorname{prox}_{\gamma g}(v)\in\Re^{m}$, indexed by $j$, $i$ and $s$, which corresponds to $y_{k+j{}\mid{}k}^{i,(s)}$, for $s=1,2,3$. Then, $(\operatorname{prox}_{\gamma g}(v))_{j,i,(1)}$ and $(\operatorname{prox}_{\gamma g}(v))_{j,i,(2)}$ are computed by virtue of the formula

[TABLE]

where $\operatorname{dist}(\cdot{}\mid{}C)$ is the distance-to-set function and the fact that $\ell^{x}(x)=W_{x}\operatorname{dist}(x\mid[x_{\min},x_{\max}])$ and $\ell^{S}(x)=W_{s}\operatorname{dist}(x\mid[x_{s},+\infty))$. The proximal operator $(\operatorname{prox}_{\gamma g}(v))_{j,i,(3)}$ is simply the projection on $\mathcal{U}$.

Note that although $g$ is prox-friendly, $g$ composed with the linear operator $H$ — as it is in (22) — is not. This is the main reason why we resort to the dual problem (31).

Given the properties of functions $f^{*}$, being differentiable with Lipschitz gradient, and $g^{*}$, being prox-friendly, we may use Nesterov’s accelerated proximal gradient method on the dual problem which produces the sequence

[TABLE]

with $y^{0}=0$, $y^{-1}=0$, $\beta^{0}=0$. In (35a) we perform an extrapolation step for some $\beta_{\nu}>0$ and we perform a dual gradient projection update on the extrapolated vector $w^{\nu}$. The extrapolation parameters $\beta-\nu$ are parametrized as $\beta^{\nu}=\theta_{\nu}(\theta_{\nu-1}^{-1}-1)$ with $\theta_{0}=\theta_{-1}=1$. Any choice of $\theta_{\nu}$ so that $1-\theta_{\nu+1}\leq\nicefrac{{\theta_{\nu+1}^{2}}}{{\theta_{\nu}^{2}}}$. A simple choice is

[TABLE]

for $\nu\geq 1$. Here we choose $\theta_{\nu+1}=\nicefrac{{1}}{{2}}(\sqrt{\theta_{\nu}^{4}+4\theta_{\nu}^{2}}-\theta_{\nu}^{2})$ which satisfies the above requirement with equality.

This algorithm fits into Tseng’s Alternating Minimization framework [69, 70]. It has been found to be suitable for embedded applications as it is relatively simple to implement and it has good convergence properties (the dual variable $y^{\nu}$ converges with rate $O(\nicefrac{{1}}{{\nu^{2}}})$ and an averaged primal iterate converges at $O(\nicefrac{{1}}{{\nu^{2}}})$ as well) [71]. It involves only matrix-vector operations (additions and multiplications) and it is numerically stable. In the following section we discuss how the involved operations can be massively parallelized in a lock-step fashion (performing the exact same operation on different memory positions) and how the algorithm can be implemented on a GPU.

4.4 Implementation

Because of the definition of $f$ in (24), the computation of the gradient of $f^{*}$ as defined in problem (30) boils down to the solution of a scenario-based optimal control problem where the only constraints are the ones defined by the system dynamics. This problem can be solved by dynamic programming leading to a Riccati-type recursion from stage $k=H_{p}$ to stage $k=0$. At each stage operations across all nodes can be fully parallelized. In particular, such a parallelization — assuming that full parallelization is supported by the hardware — equalizes the complexity of the scenario-based Riccati recursion to that of a deterministic one. A detailed exposition of the details of this procedure is available in [22, 9].

All operations involved in the computation of $\operatorname{prox}_{\gamma g}$ are element-wise operations and can be fully parallelized on a GPU; therefore, the computational cost for applying $\operatorname{prox}_{\gamma g}$ — or, what is the same — $\operatorname{prox}_{\gamma^{-1}g^{*}}$ via (33) is negligible.

5 Functionality

In this section we present RapidNet, a CUDA-C++ implementation of the accelerated proximal gradient method for the solution of scenario-based stochastic optimal control problems, particularly tailored to the needs of a drinking water network.

5.1 Software structure

The entities involved in RapidNet and the relationships among each other are illustrated in Fig. 4 which correspond to classes in the CUDA-C++ implementation of RapidNet (see also Table 1).

The DwnNetwork class stores data related to the network topology, physical constraints and network dynamics matrices in equations (8a) and (8b). The user can create instances of DwnNetwork simply by passing a JSON file, network.json, with the network information.

The ScenarioTree class models the scenario tree structure that represents the uncertainty associated with the volatile energy prices and water demands. This class describes the structure of the scenario tree by assigning a unique index to each node and storing the indexes of the children of each non-leaf node, the ancestor of all nodes except for the root node and the values $d_{k+j\mid k}^{i}$ and $\alpha_{k+j\mid k}^{i}$ at each node $(j,i)$. An instance of ScenarioTree can be generated using a JSON file, scenarioTree.json.

An SmpcController executes the accelerated proximal gradient algorithm to solve a stochastic optimal control problem at every time instant and compute the control actions to be applied to the water network. Certain functionality is delegated to Engine — a collection of utility methods — which precomputes certain quantities which are associated with the Riccati-type recursion (details can be found in [9]) for the computation of the dual gradient and manages the associated memory on GPU.

The Engine is, in turn, linked with to a DwnNetwork entity which provides all necessary technical specifications for the network (topology, dynamics, constraints) and a ScenarioTree entity which encodes the probabilistic information associated with the prediction errors. Note that the end-user does not have to create instances of Engine or directly interact with it.

A Forecaster provides to the controller estimates of the upcoming water demands and electricity prices. This is an abstract class which can be subclassed with particular model implementations (e.g., ARIMA, SVM, or any other), or the user can provide custom forecasts using any third-party software which exports its forecasts in a JSON file.

SmpcController requires certain configuration parameters which is provided by the entity SmpcConfiguration. There, the user specifies the desired tolerance, maximum number of iterations and can override other solver-specific properties.

Overall, the flow of information in RapidNet is shown in Fig. 5. The end-user initializes an SmpcController object by providing the network topology, the controller configuration and a scenario tree. During real-time operation, the controller receives the network state $x_{k}$ (which can be provided in a JSON file) and, using a demand/price forecaster, computes a control action which is applied to the system.

In Table 2 we list the main methods of RapidNet. Additional getter methods are available in each class for more advanced use-case scenarios.

5.2 GPU Implementation

GPUs were first developed for video applications and, due to the high demand in high-performance graphics, rapidly evolved to powerful hardware featuring hundreds of computation cores. Nowadays, GPUs are used for more than video processing and they are becoming popular for computational purposes including, but not limited to, environmental modeling [72, 73, 74]. By design they are well-suited for data-parallel lockstep applications where the same type of operation is applied to different memory positions. Instructions are sent to the GPU (from the CPU) in the form of compute kernels. GPUs offer unprecedented parallelization capabilities provided that the program can be parallelized in a lockstep fashion (the same operation is executed on different memory positions).

CUDA is a parallel computing framework and application programming interface for NVIDIA GPUs used for general-purpose computing. Part of the CUDA framework is cuBLAS, a parallel counterpart of the popular linear algebra library BLAS.

In RapidNet, at every time instant $k$, the method controlAction in SmpcController returns the control action that is to be applied to the water network (pump and valve set points). All computations involved in this method are either summations or matrix-vector multiplications which can be parallelized across the nodes of a stage. These multiplications are implemented using the function cublasSgemmBatched of cuBLAS and vector additions are performed using the cublasSaxpy of cuBLAS. Apart from the standard cuBLAS methods, we have defined custom kernels for the summation over the set of nodes and to evaluate projections with respect to the box constraints and the proximal operator of the distance function from the set (that is, to compute the proximal operator of $g$ as discussed in Section 4.3.

5.3 Software verification

The validation of the software is done through unit testing. A unit represents the smallest functional part of the software and unit testing involves verification of its functionality through predefined inputs and expected outputs. It assists in the debugging and maintenance of the code and facilitates the integration of the various units reliably. In lack of a standardized testing framework for CUDA applications, we developed our in-house testing framework. Moreover, using cudaMemCheck we have thoroughly tested for GPU-side memory leaks.

6 RapidNet in action: Simulation results

In this section, we present the application of RapidNet for the management of the drinking water network of the city of Barcelona, whose schematic diagram is shown in Fig. 6, using the demand data provided in [75, 9]. The network counts a total of $63$ tanks, $114$ controlled flows (by means of $75$ pumps and $39$ valves), $88$ demand sectors and $17$ mixing nodes.

6.1 Forecasting of water demands

Upcoming water demands are predicted using a radial-basis-function support vector machines (SVM) model from the literature with good predictive ability [5]. The model predicts the water demand using past demand data together with calendar data (day of the week). Validation information is provided in [5]. In Fig. 7 we show an instance of a prediction using this SVM model.

6.2 Forecasting of energy prices

Among other European countries such as Denmark and Sweden, Austria implements a deregulated energy market which induces a volatile stream of electricity prices. Time series of energy prices in Austria have become public by EXAA (Energy Exchange Austria; a central European energy exchange) and are available online at http://www.exaa.at/de/marktdaten/historische-daten (Accessed on ).

In lack of electricity price data from Spain or relevant data from the water network of Barcelona, we used price data from Austria as an indicative dataset.

Out of $8784$ hourly price data which are available for the year 2016, we excluded the last $2000$ data points to be used for testing and using the rest of the data we built an ARIMA$(24,1,4)$ model. ARIMA models have been previously used for the short-term prediction of electricity prices in the day-ahead market [51]. Using a Monte-Carlo method, a set of $10^{4}$ independent scenarios were generated and, subsequently, these were reduced into a scenario tree using the method described in [55]. An instance of a prediction using the trained ARIMA model along with the associated scenario tree is shown in Fig. 8.

The residuals of the model where found to be uncorrelated at the confidence level of $99.9\%$. Indeed, the residuals pass the Ljung-Box Q-test of uncorrelatedness with p-value equal to $1.0000$. The model was selected using the Akaike information criterion (AIC) with value $1.523$.

It is natural to expect that the upcoming energy prices can be only predicted up to moderate accuracy as they do not follow a regular pattern and are influenced by many market-related parameters. Despite our limited predictive capacity, we shall show in Section 6.3 that by taking into account this volatility using stochastic model predictive control we do mitigate the effect of the price uncertainty leading to a more economic operation of the water network.

6.3 Closed-loop simulations

In this section we present closed-loop simulation results on the water network where the sampling time is equal to $1\,\mathrm{hr}$ and the prediction horizon of the SSMPC was fixed to $H_{p}=24$. The weight parameters used to tune the SSMPC are $W_{\alpha}=10^{6}$, $W_{u}=1.3\cdot 10^{7}\cdot I$, $W_{s}=10^{5}$ and $W_{x}=10^{8}$. Note that all units used in this section are SI units (flows in $\nicefrac{\mathrm{m^{3}}}{\mathrm{s}}$ and volumes in $\mathrm{m^{3}}$).

The system was simulated for a period of $H_{s}=168$ time instants, which corresponds to one week of operation. In order to assess the performance of the closed-loop operation, we use three KPIs, in particular: (i) the economic index which is defined as

[TABLE]

which provides an estimation of the average hourly cost of operation of the water network, (ii) the safety index, which is defined as

[TABLE]

which quantifies the total weekly violation of the safety constraint (the “soft” requirement that $x_{k}\geq x_{s}$), and, last, (iii) the complexity index which is the worst-case time required to solve the SSMPC optimization problem up to the specified accuracy, which is defined as

[TABLE]

where $\tau_{k}$ is the computation time required by the solver for solving the SSMPC problem at time instant $k$. The maximum runtime is of higher importance that the average runtime in applications to verify that a decision can be made within the available time period. All three indices are defined so that low values are preferred.

In order to justify the need for and underline the importance of taking into account the uncertainty associated with electricity prices, we performed two sets of simulations where in one case we disregarded that volatility (the SSMPC used only the nominal predictions of the electricity prices). We may observe in Fig. 9 that the hourly operation of the network becomes more expensive by approximately $+5\%$, therefore, the proposed stochastic control approach can lead to significant economic savings for the network operator. For example, for the case of $631$ scenarios, the operation cost when the price uncertainty is disregarded is € 4947/hour, whereas with the proposed approach which takes into account the price uncertainty it drops to € 4748/hour — a saving of € 199/hour, that is, a cost reduction of $4\%$.

We may also see that as more scenarios are considered in the SSMPC formulation, the average cost of operation plummets at around $194$ scenarios, where, however, the safety index is still high. In order to obtain a safer operation we need to pay the “cost of safe operation”: Indeed, at $631$ scenarios, the operation of the network is more expensive compared to the case of $194$ scenarios, but $\mathrm{KPI}_{S}$ is at a minimum. This reflects a trade-off between the economic and safe operation of the network.

At $631$ scenarios, the SSMPC optimization problem involves as many as $2,306,133$ primal variables and $3,126,960$ dual variables. RapidNet solves it with $\mathrm{KPI}_{\tau}=82.7\,\mathrm{s}$ as compared to $\mathrm{KPI}_{\tau}=719\,\mathrm{s}$ for the popular commercial solver Gurobi v7.0 at the same level of accuracy (tolerance $5\cdot 10^{-2}$). The solvers cplex and mosek were also tested, but Gurobi outperformed them consistently. As we may see in Fig. 10, RapidNet exhibits a lower complexity index by approximately an order of magnitude.

7 Conclusions and Future Work

We presented the integrated software solution RapidNet for the control of drinking water networks which accounts for uncertainty both in predicted water demand and in predicted electricity prices in the day-ahead energy market.

RapidNet is a highly inter-operable software as it can be interfaced using JSON files which follow a standard API which is detailed in the software documentation. It can be combined with any forecaster as the controller does not need to know the mechanism with which the forecasts are produced. RapidNet features a parser implemented in MATLAB which allows the conversion of EPANET .inp files to JSON files. It is an open-source and free software which is distributed under the terms of the GNU LGPL v3.0 licence and can be downloaded at https://github.com/GPUEngineering/RapidNet.

In this paper we have shown via simulations that when a water network is operated in the context of a volatile energy market, considerable savings can be obtained by using a predictor of the upcoming energy prices and taking into account the associated price-related uncertainty in the formulation of the scenario tree. Alongside, we advocate that RapidNet can transform SSMPC from a powerful (but too complex) theoretical development to control engineering practice and enabling the solution of very large SSMPC scale problems and their seamless integration into the control system of the water network.

This work focuses on flow-based water distribution networks such as the network of the city of Barcelona. Pressure-based networks lead to optimization problems with nonconvex constraints. These can be approximated by solving a constraints satisfaction problem (CSP) as discussed in [60]. The problem can be then transformed into a convex SSMPC problem which can be solved by RapidNet. The operator splitting concept has been proven to apply to sums of nonconvex functions which allows the solution of nonconvex optimal control problems with input constraints [76]. This allows the application of recent developments in nonconvex optimization such as PANOC [77] for the solution of nonlinear MPC problems.

Future developments in RapidNet will involve the implementation of new faster parallelizable algorithms which make use of quasi-Newton directions based on our recent theoretical work [78] and the exploitation of multiple-GPU architectures. We shall introduce a library of available water network entities from the literature as well as reported water demand and energy price models for different water networks and energy markets.