Local Water Storage Control for the Developing World

TL;DR

This paper introduces a stochastic control framework for managing water storage in Indian cities with intermittent supply, balancing costs, health risks, and water wastage, using a novel Binary Dynamic Search algorithm.

Contribution

It presents a new stochastic control model and a Binary Dynamic Search algorithm for optimizing water storage and disinfection timing in developing-world urban settings.

Findings

Effective control policy for water storage and reset timing.

Reduction in health risks from water stagnation.

Cost-efficient water management strategy.

Abstract

Most cities in India do not have water distribution networks that provide water throughout the entire day. As a result, it is common for homes and apartment buildings to utilize water storage systems that are filled during a small window of time in the day when the water distribution network is active. However, these water storage systems do not have disinfection capabilities, and so long durations of storage (i.e., as few as four days) of the same water leads to substantial increases in the amount of bacteria and viruses in that water. This paper considers the stochastic control problem of deciding how much water to store each day in the system, as well as deciding when to completely empty the water system, in order to tradeoff: the financial costs of the water, the health costs implicit in long durations of storing the same water, the potential for a shortfall in the quantity of stored versus demanded water, and water wastage from emptying the system. To solve this problem, we develop a new Binary Dynamic Search (BiDS) algorithm that is able to use binary search in one dimension to compute the value function of stochastic optimal control problems with controlled resets to a single state and with constraints on the maximum time span in between resets of the system.

Equations

$$\displaystyle\begin{aligned} x_{n+1}&=\begin{cases}(x_{n}+u_{n}-d_{n})^{+},&\text{if }f_{n}=0\\ (u_{n}-d_{n})^{+},&\text{if }f_{n}=1\end{cases}\\ &=\big{(}(1-f_{n})\cdot x_{n}+u_{n}-d_{n}\big{)}^{+}\end{aligned}$$

$$\displaystyle\begin{aligned} t_{n+1}&=\begin{cases}t_{n}+1,&\text{if }f_{n}=0\\ 1,&\text{if }f_{n}=1\end{cases}\\ &=\big{(}t_{n}\cdot(1-f_{n})+1\big{)}\end{aligned}$$

$$\min$$

$$\displaystyle\hbox to0.0pt{$\hskip 42.67912ptq(t_{n})\cdot(\xi_{n}+u_{n}-d_{n})^{+}+c_{f}f_{n}x_{n}\Big{)}\Big{]}$\hss}$$

$$x_{n+1}=({\xi }_n+u_n-d_n{)}^{+},$$

$$t_{n+1}=({\tau }_n+1),$$

$${\xi }_n=(1-f_n)\cdot x_n,$$

$${\tau }_n=(1-f_n)\cdot t_n,$$

$$t_n\le k,$$

$$f_n\in \{0,1\},u_n\in {\mathbb{R}}_{+},$$

$$f(d)=\sum _{k=1}^{\infty }\frac{(\lambda \cdot \mu {)}^k{d}^{k-1}{e}^{(-\lambda -\mu d)}}{k!(k-1)!}1[d>0]+e^{-\lambda }\delta (d)$$

$$f(d)=e^{(-\lambda -\mu d)}{I}_1(2\sqrt{\lambda \mu d})\sqrt{\frac{\lambda \mu }{d}}1[d>0]+e^{-\lambda }\delta (d)$$

$$x_{n+1}$$

$$t_{n+1}$$

$${\xi }_n$$

$${\tau }_n$$

$$\min$$

$$(\text{ref}eqn:dyngen),t_n\le k,u_n\in \mathcal{U},$$

$$J(\zeta ,0)$$

$$J(x,t)$$

$$\displaystyle\hskip 14.22636pt\min_{u\in\mathcal{U}}\mathbb{E}\Big{[}g(x,t,u,w)+\gamma J(h(x,t,u,w),t+1)\Big{]}\Big{\}}$$

$$J(x,k)$$

$$\min$$

$$\displaystyle\hbox to0.0pt{$\qquad\qquad\qquad\qquad\qquad\qquad\qquad s(x_{n},t_{n},w_{n})\cdot r_{n}\Big{)}\Big{]}$\hss}$$

$$(\text{ref}eqn:dyngen),t_n\le k,u_n\in \mathcal{U},$$

$$\min$$

$$\displaystyle\hbox to0.0pt{$\quad\textstyle\sum_{n=1}^{\infty}\gamma^{n}\Big{(}g(\xi_{n},\tau_{n},u_{n},w_{n})+s(x_{n},t_{n},w_{n})\cdot r_{n}\Big{)}\Big{]}$\hss}$$

$$(\text{ref}eqn:dyngen),t_n\le k,u_n\in \mathcal{U},$$

$$\min$$

$$\displaystyle\hbox to0.0pt{$\textstyle\sum_{n=1}^{\infty}\gamma^{n}\Big{(}g(\xi_{n},\tau_{n},u_{n},w_{n})+s(x_{n},t_{n},w_{n})\cdot r_{n}\Big{)}\Big{]}$\hss}$$

$$(\text{ref}eqn:dyngen),t_n\le k,u_n\in \mathcal{U},$$

$${\rm Y}(v)$$

$$V(x,t,v)$$

$$\displaystyle\min_{u\in\mathcal{U}}\mathbb{E}\Big{[}g(x,t,u,w)+\gamma V(h(x,t,u,w),t+1,v)\Big{]}\Big{\}}$$

$$V(x,k,v)$$

$$\displaystyle V(x,t,J_{0})=\min\Big{\{}J_{0}+\mathbb{E}\big{(}s(x,t,w)\big{)},$$

$$\displaystyle\qquad\min_{u\in\mathcal{U}}\mathbb{E}\Big{[}g(x,t,u,w)+\gamma V(h(x,t,u,w),t+1,J_{0})\Big{]}\Big{\}}.$$

$$\nu ={\rm Y}(\nu )$$

$$V(x,t,\nu )$$

$$\displaystyle\min_{u\in\mathcal{U}}\mathbb{E}\Big{[}g(x,t,u,w)+\gamma V(h(x,t,u,w),t+1,\nu)\Big{]}\Big{\}}$$

$$V(x,k,\nu )$$

$$\overline{J}=\textstyle\Big{(}\frac{1}{1-\gamma}\Big{)}\min_{u\in\mathcal{U}}\mathbb{E}\Big{[}g(\zeta,0,u,w_{0})+\\ \gamma\cdot s(h(\zeta,0,u,w_{0}),1,w_{1})\Big{]},$$

$$J_{0}\leq\min_{u}\mathbb{E}\Big{[}g(\zeta,0,u,w_{0})+\\ \gamma(J_{0}+s(h(\zeta,0,u,w_{0}),1,w_{1})\Big{]}.$$

$$\textstyle\max_{x}\big{|}V(x,k,v^{*})-J(x,k)\big{|}=\big{|}v^{*}-J_{0}\big{|}\leq\epsilon$$

Full text

I INTRODUCTION

Safe drinking water is not widely available in much of the developing world. For instance, the majority of cities in India do not have continuous water supply within the existing distribution infrastructure, and on average citizens get water for a 3 to 4 hour period during the day [1, 2]. As a result, most households use communal and personal water storage containers to ensure adequate water supply throughout the day; these water storage containers are filled every 1 or 2 days during periods of water availability [2, 3]. However, water stored in these containers often becomes contaminated [4, 5, 6]: Longer durations of storing the same water leads to increased amounts of bacterial and viral pathogens.

In this paper, we formulate the problem of deciding how much water to store each day in a local water storage container as a stochastic optimal control problem. The challenge is how to tradeoff the financial costs of purchasing the water, the degradation in water quality (and associated implicit health costs) from long storage durations, the potential for a shortfall in the quantity of stored versus demanded water, and water wastage. The inherent uncertainty in daily water demand makes our choice of stochastic optimal control natural. Furthermore, our goal is to construct a control policy that is easily implementable via a lookup table that could be distributed through paper pamphlets or internet websites.

In addition to studying the stochastic optimal control problem associated with local water storage container, our second contribution in this paper is to develop a new algorithm for solving stochastic optimal control problems with the specific structure of controlled resets to a single state and with constraints on the maximum time in between system resets. Unlike value and policy iteration, which require finding a fixed point in an infinite dimensional functional space [7], we develop a Binary Dynamic Search (BiDS) algorithm that converts the control problem into finding a fixed point in a vector-space using binary search. Though value and policy iteration converge exponentially, the convergence rate can be practically slow when the discount factor is close to 1 [7, 8]. Our BiDS algorithm can compute the value functions of problems with controlled resets to a single state using substantially less computation because it uses binary search.

I-A Energy Storage Control

Though local water storage is nominally similar to control of (mainly electrical) energy storage [9, 10, 11, 12, 13, 14, 15, 16], there are several differences in the underlying physics and domain settings that make the two problems very distinct. The biggest difference is that stored energy dissipates (or decays) in quantity over time, whereas local water storage tanks are typically closed and hence have minimal water evaporation [4, 5, 6]. Another important contrast is that the quality of stored water decays [4, 5, 6], whereas the quality of stored energy is directly related to the stored quantity. Because the quality of stored water decays, our formulation considers the possibility of flushing all water from the storage container; in contrast, energy storage does not have a corresponding control action to empty the storage.

I-B Inventory and Supply Chain Management

Our problem of local water storage control is also similar to inventory and supply management problems [17, 18]. More specifically, it is related to the setting of managing perishable inventory [19, 20, 21, 22, 23, 24], which has been mainly explored in the context of healthcare inventory management. In these perishable inventory models, it is assumed that expired units of inventory can be disposed of individually, which allows salvage and disposal costs to be incorporated into the inventory’s holding cost; hence, policies of strategic inventory removal can be employed [21, 22, 23, 24]. However, in our setting, since fresh water may mix with still water, individuals cannot dispose of single units of their stock strategically. Instead, individuals must make a decision of either keeping or disposing of their entire stock of water to reduce the risk of drinking contaminated water.

I-C Numerical Solution of Stochastic Optimal Control

The workhorse for numerical solution of stochastic optimal control is value and policy iteration [7], which require finding a fixed point in an infinite dimensional space. For high-dimensional state spaces, exact calculation for these approaches is numerically difficult (i.e., the so-called curse of dimensionality). And so a number of approximate dynamic programming approaches [7, 25, 26, 27, 28, 29] have been developed, where the main idea is that approximate dynamic programming performs computation with an inexact but tractable representation of the value function or policy.

However, the convergence rate of (exact or approximate) value and policy iteration is governed by the discount factor, and so convergence is slow for discount factors close to 1 [7, 8]. This motivates development of our BiDS algorithm, which uses the special structure of stochastic control problems with controlled resets to a single state and with a constraint on the amount of time in between resets. Our BiDS algorithm finds a fixed point in a vector-space using binary search. Since the state space is small for our local water storage control problem, our numerical results solve the exact dynamic program; however, in principle our BiDS algorithm could be used for approximate dynamic programming by using inexact representations of the value function.

I-D Outline

We first describe in Sect. II the stochastic optimal control problem for managing local water storage. We also describe a stochastic model for daily water demand. Section III generalizes this into a broader class of reset control problems (i.e., stochastic optimal control with controlled resets to a single state and with a constraint on the maximum timespan in between system resets). Next we develop essential elements of the BiDS algorithm for solving reset control problems, and prove the correctness of the value function and control policy computed by BiDS. We conclude with Sect. IV, which uses BiDS to numerically solve the local water storage control problem for a specific set of numeric problem parameters.

II Description of Water Storage Model

In this section, we first describe the stochastic optimal control problem associated with managing a local water storage container. A key element of this model is the stochastic distribution used to model daily water demand. We discuss one common stochastic model for water demand in a household and provide an algorithm to generate random samples from this distribution, which is useful for numeric computations when solving the stochastic optimal control.

II-A Stochastic Optimal Control Formulation

The subscript $n$ denotes the corresponding variable for the $n$-th day. We define the following variables:

• •

$x_{n}$ - amount of water stored in the tank;

• •

$t_{n}$ - number of days since the tank was last emptied;

• •

$u_{n}$ - control of adding new water to the tank;

• •

$d_{n}$ - amount of water consumed;

• •

$f_{n}$ - control of emptying the entire tank;

We also have the following dynamics to describe the water storage tank:

[TABLE]

where $(x)^{+}=\max\{x,0\}$ and $(x)^{-}=\min\{x,0\}$. Let $q:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ be a strictly increasing nonnegative function, and suppose $c,p,c_{f}\in\mathbb{R}_{+}$ are fixed positive constants. The goal is to solve the following stochastic control problem

[TABLE]

The interpretation is that $p\cdot(\xi_{n}+u_{n}-d_{n})^{-}$ is a penalty for not having stored enough water to meet demand, $q(t_{n})\cdot(\xi_{n}+u_{n}-d_{n})^{+}$ measures the health cost due to storing excess water, $c$ is the cost of the water, $c_{f}$ is the waste penalty for flushing too much water, $\gamma\in[0,1)$ is a discount factor, and the constraint $\textstyle\sum_{j=1}^{k}f_{n+j}\geq 1$ for $n\geq 0$ ensures that the storage tank is fully emptied at least once every $k$ days to maintain safe water quality. Since the function $q$ is strictly increasing, the $q(t_{n})$ term indicates that the quality of water deteriorates as time passes between emptying the tank.

II-B Water Demand Model

A commonly used stochastic distribution for water demand $d_{n}$ is the Poisson Rectangular Pulse (PRP) model [30], and the simplest version of this model is that the number of usage events has a Poisson distribution with rate $\lambda$ while the duration of water usage is the sum of exponential distributions with rate $\mu$. The probability density function (PDF) of $d_{n}$ is given by

[TABLE]

where $\lambda>0$ is the usage event rate, $\mu>0$ is the usage duration rate, $\mathbf{1}[\cdot]$ is the indicator function, and $\delta(\cdot)$ is the Dirac delta function. An alternative expression [30] for the PDF of this distribution is

[TABLE]

where $I_{1}$ is the modified Bessel function of the first type. Though this is an unusual distribution, generating random numbers is relatively straightforward:

Proposition 1

Algorithm 1 generates random numbers whose distribution is that of the PRP model.

Proof:

This is true by construction of the algorithm and by the definition of the PRP model. ∎

III Reset Control Model and BiDS Algorithm

An important characteristic of the local water storage control problem is that choosing $f_{n}=1$ in effect causes the state of the system to be reset to $(x_{n},t_{n})=(0,0)$. In this section, we generalize the local water storage control problem into a broader class of stochastic optimal control problems with controlled resets and with constraints on the maximum time between rests. Then we develop a new algorithm that uses binary search to compute the value function and hence solve this particular stochastic optimal control problem.

III-A Stochastic Optimal Control with Controlled Resets

Let the subscript $n\in\mathbb{Z}$ denote time, and consider the discrete-time dynamical system

[TABLE]

where $x_{n}\times t_{n}\in\mathbb{R}^{p}\times\mathbb{Z}_{+}$ are states, $\xi_{n}\times\tau_{n}\in\mathbb{R}^{p}\times\mathbb{Z}_{+}$ are psuedo-states, $u_{n}\times r_{n}\in\mathbb{R}^{q}\times\{0,1\}$ are controls, $w_{n}\in\mathbb{R}^{m}$ are i.i.d. random variables, and $h:\mathbb{R}^{p}\times\mathbb{Z}_{+}\times\mathbb{R}^{q}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a deterministic function. The interpretation is that the control $r_{n}=1$ resets the system to a known initial state $\zeta\in\mathbb{R}^{p}$, the state $t_{n}$ keeps track of how many time steps have passed since the last system reset, and $h$ describes the dynamics when there is no reset.

Let $\gamma\in[0,1)$ be a discount factor, and suppose that our goal is to solve the stochastic control problem

[TABLE]

where $g:\mathbb{R}^{p}\times\mathbb{Z}_{+}\times\mathbb{R}^{q}\times\mathbb{R}^{m}\rightarrow\mathbb{R}_{+}$ is a nonnegative and continuous stage cost, $s:\mathbb{R}^{p}\times\mathbb{Z}_{+}\times\mathbb{R}^{m}\rightarrow\mathbb{R}_{+}$ is a nonnegative and continuous reset cost with $s(\zeta,t,w)\equiv 0$, the constraint $u_{n}\in\mathcal{U}$ restricts the possible control actions, and the constraint $t_{n}\leq k$ requires the system be reset at least once every $k$ time steps. This stochastic optimal control problem has special structure that can be used to develop a new approach for its solution.

We first characterize the optimal value function for (7). Let $J(x,t)$ be the optimal cost to go, and define $J_{0}=J(\zeta,0)$. Recall that $J(x,t)$ is defined as the minimum value of (7) for the initial conditions $x_{0}=x$ and $t_{0}=t$.

Proposition 2

The dynamic programming equations for (7) are given by

[TABLE]

where the middle $J(x,t)$ holds for all $x$ and $t=0,\ldots,k-1$.

Proof:

We prove this by induction on $t$. Recalling that $s(\zeta,0,w)=0$, we rewrite (7) for $(x_{0},t_{0})=(\zeta,0)$ as

[TABLE]

But since $(x_{0},t_{0})=(\zeta,0)$, we have $\xi_{0}=\zeta$ and $t_{1}=1$ for both $r_{0}=0$ and $r_{0}=1$. Thus the above simplifies to the expression for $J(\zeta,0)$ given in (8).

Next consider the case with $(x_{0},t_{0})=(x,k)$. We must have $r_{0}=1$ because otherwise $t_{1}=k+1$ which violates the constraint $t_{n}\leq k$. So $\xi_{0}=\zeta$, and we can rewrite (7) as

[TABLE]

Clearly we have $J(x,k)=J(\zeta,0)+\mathbb{E}[s(x,k,w)]$, since the reset cost $s(x,k,w)$ does not depend on the control $u_{0}$.

Now consider the case with $(x_{0},t_{0})=(x,k-1)$, and let $J_{0}(x,k-1)$ and $J_{1}(x,k-1)$ denote the optimal cost to go under $r=0$ and $r=1$, respectively. By Bellman’s principle of optimality we have $J(x,k-1)=\min\{J_{0}(x,k-1),J_{1}(x,k-1)\}$. For $J_{1}(x,k-1)$, we have $\xi_{0}=\zeta$ by definition, and so a similar argument to the one for $J(x,k)$ gives that $J_{1}(x,k-1)=J_{0}+\mathbb{E}[s(x,k-1,w)]$. Next, observe that we can express $J_{0}(x,k-1)$ as

[TABLE]

Simplifying gives $J_{0}(x,k-1)=\min\mathbb{E}[\textstyle g(x,0,u_{0},w_{0})+\gamma J(h(x,k-1,u_{0},w_{0}),k)]$. Hence $J(x,k-1)$ satisfies the middle equation for $J(x,t)$ in (8). Now suppose that the middle equation for $J(x,t)$ in (8) is satisfied up to $0<t+1\leq k-1$: We will show this implies that the middle equation for $J(x,t)$ in (8) is satisfied for $t$. Indeed, using the same argument as for the base case gives this desired step. The final result follows by applying induction. ∎

III-B Design of the BiDS Algorithm

To solve these dynamic programming equations and compute the optimal policy, we develop our Binary Dynamic Search (BiDS) algorithm that is presented in Algorithm 2. The key feature of BiDS is that unlike value iteration [7], which requires finding a fixed point over the space of functions by repeated application of a contraction mapping, BiDS performs a binary search in order to find a fixed (vector-valued) point of an appropriately constructed mapping.

Our main result about the BiDS algorithm concerns its finite termination and convergence properties. To show this, we first study the recursion calculations within BiDS.

Proposition 3

Consider the recurrence relation given by

[TABLE]

where the middle equation for $V(x,t,v)$ holds for all $x$ and $t=0,\ldots,k-1$. Then $\Upsilon(v)$ is monotone increasing in $v$, we have that $V(x,t,J_{0})=J(x,t)$, and $v=J_{0}$ is the unique (vector-valued) fixed point of $\Upsilon(v)$.

Proof:

To prove the first result, note that $V(x,k,v)$ is monotone increasing by construction. Next, suppose $V(x,t+1,v)$ is monotone increasing in $v$. The first term in the $\min$ defining $V(x,t,v)$ is monotone increasing by construction, and the second term in the $\min$ defining $V(x,t,v)$ is monotone increasing by assumption. Thus $V(x,t,v)$ is monotone increasing in $v$, and this result extends to all $t=0,\ldots,k-1$ by induction. Lastly, observe that $\Upsilon(v)$ is monotone increasing in $v$ since its first term is independent of $v$ and since the second term $V(h(\zeta,0,u,w),1,v)$ was shown to be monotone increasing in $v$.

To prove the second result, we assume $v=J_{0}$ and use induction on $t$. First note $V(x,k,J_{0})=J_{0}+\mathbb{E}(s(x,k,w))=J(x,k)$ by Proposition 2. Next suppose that $V(x,t+1,J_{0})=J(x,t+1)$, and note that by definition

[TABLE]

Comparing to Proposition 2 shows $V(x,t,J_{0})=J(x,t)$, and so the result for $V(x,t,J_{0})$ for all $x$ and $t=0,\ldots,k-1$ follows by induction. Lastly, inserting the relation $V(x,1,J_{0})=J(x,1)$ into the definition for $\Upsilon(J_{0})$ gives $\Upsilon(J_{0})=J_{0}$ by comparison with Proposition 2.

To show the third result, let $\nu$ be any value such that $\nu=\Upsilon(\nu)$. (Note that we have shown above that such a $\nu$ must exist.) Then we can rewrite (12) as

[TABLE]

Comparing these to Proposition 2 shows that $V(x,t,\nu)$ is a fixed point of the dynamic programming equations for (7). Thus $V(x,t,\nu)=J(x,t)$, and $J(x,t)$ is unique by Proposition 4.1.5 of [7]. As a result, we have $\nu=\Upsilon(\nu)=V(\zeta,0,\nu)=J_{0}$. This means that $\nu$ is uniquely defined. ∎

Before we prove our main result about the BiDS algorithm, we need an additional lemma.

Lemma 1

Suppose $|a^{\prime}-a|\leq\epsilon$ and $|b^{\prime}-b|\leq\epsilon$. Then we have that $|\min\{a^{\prime},b^{\prime}\}-\min\{a,b\}|\leq\epsilon$.

Proof:

Let $c^{\prime}=\min\{a^{\prime},b^{\prime}\}$ and $c=\min\{a,b\}$. We have four cases: (1) $c^{\prime}=a^{\prime}$ and $c=a$, then $|c^{\prime}-c|\leq\epsilon$; (2) $c^{\prime}=b^{\prime}$ and $c=b$, then $|c^{\prime}-c|\leq\epsilon$; (3) $c^{\prime}=a^{\prime}$ and $c=b$, then $c^{\prime}-c=a^{\prime}-b\geq a^{\prime}-a\geq-\epsilon$ and $c^{\prime}-c=a^{\prime}-b\leq b^{\prime}-b\leq\epsilon$; and (4) $c^{\prime}=b^{\prime}$ and $c=a$, then $c^{\prime}-c=b^{\prime}-a\geq b^{\prime}-b\geq-\epsilon$ and $c^{\prime}-c=b^{\prime}-a\leq a^{\prime}-a\leq\epsilon$. ∎

We can now show our main result about BiDS:

Theorem 1

The BiDS algorithm terminates in finite time when $\epsilon>0$, and the solution $v^{*}$ returned by BiDS is such that \max_{x}\big{|}V(x,t,v^{*})-J(x,t)\big{|}\leq\epsilon for all $t=0,\ldots,k$.

Proof:

Since the state costs and reset costs are nonnegative, Proposition 4.1.5 from [7] implies that the value function $J(x,t)$ is the unique fixed point of the dynamic programming equations (8) for (7). Our first intermediate claim is that $J_{0}$ is finite. In particular, define

[TABLE]

and observe that $\overline{J}$ is finite since the objective is nonnegative and the optimization problem is always feasible. Next note $J_{0}\geq 0$ since the stage costs and reset costs are nonnegative. Furthermore, using the dynamic programming equations (8) from Proposition 2 gives that $J(x,1)\leq J_{0}+\mathbb{E}(s(x,1,w))$. Again using the dynamic programming equations (8) gives

[TABLE]

Solving for $J_{0}$ gives $J_{0}\leq\overline{J}$. Since we showed $J_{0}\in[0,\overline{J}]$, this implies $J_{0}$ is finite. Next observe that the BiDS algorithm is a binary search over the range $[0,\overline{J}]$. Recalling that $V(\zeta,0,v)$ is monotone increasing in $v$ by Proposition 3, we must have that the binary search terminates in a finite number of steps when $\epsilon>0$. This proves the first part of the theorem.

We use induction to prove the second part of the theorem. First note that by Propositions 2 and 3 we have

[TABLE]

where the inequality follows by the termination condition of BiDS. Now suppose we have that $\max_{x}|V(x,t+1,v^{*})-J(x,t+1)|\leq\epsilon$. Then a similar calculation gives \max_{x}\big{|}v^{*}+\mathbb{E}(s(x,t,w))-J_{0}-\mathbb{E}\big{(}s(x,t,w)\big{)}\big{|}\leq\epsilon and that \max_{x,u}\big{|}g(x,t,u,w)+\gamma V(h(x,t,u,w),t+1,v^{*})-g(x,t,u,w)-\gamma J(h(x,t,u,w),t+1)\big{|}\leq\epsilon. And so Lemma 1 implies that $\max_{x}|V(x,t,v^{*})-J(x,t)|\leq\epsilon$ since expectation preserves distances. Thus $\max_{x}|V(x,t,v^{*})-J(x,t)|\leq\epsilon$ for all $t=0,\ldots,k$ by induction. ∎

We can use the BiDS algorithm to compute the optimal policy by taking the minimizing arguments of $V(x,t,v^{*})$.

Corollary 1

*Let $v^{*}$ be the value returned by BiDS for a given $\epsilon>0$. If $\textstyle(u,r)^{*}(x,t,\epsilon)\in\operatorname*{arg\,min}V(x,t,v^{*})$ for $(x,t)\neq(\zeta,0)$ and $\textstyle(u,r)^{*}(\zeta,0,\epsilon)\in\operatorname*{arg\,min}\Upsilon(v^{*})$, then $\textstyle\limsup_{\epsilon\downarrow 0}(u,r)^{*}(x,t,\epsilon)\subseteq\operatorname*{arg\,min}J(x,t)$ for $t=0,\ldots,k$. *

Proof:

Theorem 1 says $V(x,t,v^{*})$ converges uniformly to $J(x,t)$, and that $\Upsilon(v^{*})$ converges uniformly to $J(\zeta,0)$. Next recall Fatou’s lemma [31] implies lower semicontinuity (lsc) of nonnegative functions is preserved under expectation. This means $V(x,k,v^{*})$ is lsc in $x$. Now suppose $V(x,t+1,v^{*})$ is lsc in $x$. Then $v^{*}+\mathbb{E}(s(x,t,w))$ and $\mathbb{E}[g(x,t,u,w)+\gamma V(h(x,t,u,w),t+1,v)]$ are lsc in $(x,u)$ by Fatou’s lemma. Since lsc is preserved under minimization [32], this means $V(x,t,v^{*})$ is lsc in $x$. By induction, $V(x,t,v^{*})$ for $t=0,\ldots,k$ is lsc in $x$. Hence, Proposition 7.15 of [32] says the objectives of $V(x,t,v^{*})$ and $\Upsilon(v^{*})$ epi-converge. So the result follows by Theorem 7.33 of [32]. ∎

In words, the above result says that as $\epsilon\rightarrow 0$ then BiDS generates a sequence of $v^{*}$ such that the minimizers of $V(x,t,v^{*})$ and $\Upsilon(v^{*})$ converge to minimizers of $J(x,t)$.

IV Numerical Solution of Water Storage Model

In this section, we use BiDS to numerically solve the water storage model (3) from Sect. II. Computations were performed in MATLAB 2016b on a laptop computer with a 2.4GHz processor and 16GB of RAM. We generated random demand from the PRP model using Algorithm 1 with Poisson rate $\lambda=40$ (i.e., average of 40 water usage events for a family) and exponential rate $\mu=2/\text{L}$ (i.e., average water usage event consumes $500$mL of water). The costs in Indian Rupees were 0.25/L, 0.5/L, 0.5/L for the purchasing, shortage, and flushing costs of water, respectively. The holding cost function in Indian Rupees was $q(t_{n})=15t_{n}$/L, which satisfies the monotonicity assumption and is high compared to the other costs to represent a sharp deterioration in the quality of water.

The BiDS algorithm was used to solve the local water storage problem with a maximum of $k=6$ days between flushing and with a discount factor of $\gamma=0.8$, and the results are shown in Fig. 1. These results show the optimal policy has a clear structure, and this control policy is easily implementable via a lookup table that could be distributed through paper pamphlets or internet websites. At each $t$, there are at most three thresholds that separate the control policy into four regions. The first region is where the amount of water in the tank is so small that the flushing cost is negligible, and hence it is optimal to flush the tank and order a one step optimal quantity. Next is a region with enough water in the tank were it is optimal to order up to some quantity and still consume most of the water before it becomes unsafe to drink. Then is a do-nothing region where it is not optimal to flush the tank or order any more water since there is enough water in the tank to satisfy possible future consumption. Finally, there is a region where there is so much water in the tank that it is almost certain to never be consumed, and thus it is optimal to flush this water from the tank and reorder up to some optimal quantity.

V Conclusion

In this paper, we studied the problem of managing local water storage in the developing world. We first presented a stochastic optimal control formulation of this problem, and gave an algorithm for sampling from a PRP model that describes daily water demand. We showed this water storage problem can be generalized to stochastic optimal control problems with controlled resets to a single state and with constraints on the maximum time between resets. To solve this class of problems, we developed the BiDS algorithm that performs binary search to obtain the optimal value function and corresponding optimal policy. We showed this algorithm converges in finite time and is guaranteed to find the optimal value with $\epsilon$ precision. Finally, we concluded by applying the BiDS algorithm to numerical solve the problem of managing local water storage. Our results suggest that the optimal control policy for this problem has a well defined structure, which we leave to be studied in future work.