Distributionally robust SDDP

Abstract

We study a version of stochastic dual dynamic programming (SDDP) with a distributionally robust objective. The classical SDDP algorithm uses a finite (nominal) probability distribution for the random outcomes at each stage. We modify this by defining a distributional uncertainty set in each stage to be a Euclidean neighbourhood of the nominal probability distribution. We derive a formula for the worst-case expectation of future costs over this set that can be applied in the backward pass of SDDP. We verify the correctness of this algorithm, show its almost sure convergence under standard assumptions, and illustrate it by applying it to a model of the New Zealand hydrothermal electricity system.

Appendices

Appendix A

This appendix contains proofs of all the results in the paper. We begin by proving some technical lemmas. Recall the problem

$$\begin{aligned} \begin{array}{rl} \rho _{r}(z)=\max &{}\displaystyle \sum _{i=1}^{m}z_{i}p_{i} \\ \text {s.t.} &{} \displaystyle \sum _{i=1}^{m}p_{i}=1, \\ &{} \left\| p-q\right\| _{2}\le r, \\ &{} p\ge 0. \end{array} \end{aligned}$$

Lemma 1

Given $z_{i}$, $i=1,2,\ldots ,m$, $\rho _{r}(z)$ is a continuous and nondecreasing function of $r\ge 0$.

Proof

Observe that given $r\ge 0$ and $z_{i}$, $i=1,2,\ldots ,m$, P is a convex optimization problem with a compact feasible region, so it has an optimal solution, with optimal value denoted $\rho _{r}(z)$. Moreover for fixed z, the optimal value function is a concave increasing function of r. By Rockafellar (1972, Theorem 10.1) $\rho _{r}(z)$ is therefore continuous on $\{r:r>0\}$. It is also easy to show that $\lim _{_{r\rightarrow 0}}\rho _{r}(z)=\sum _{i=1}^{m}q_{i}z_{i}$, so $\rho _{r}(z)$ is also continuous at $r=0$. $\square$

Consider the problem

$$\begin{aligned} \begin{array}{rl} \text {P(}n\text {):}\min &{}\displaystyle -\sum _{i=1}^{n}z_{i}y_{i} \\ \text {s.t.} &{} \displaystyle \sum _{i=1}^{n}y_{i}=a, \\ &{} \displaystyle \sum _{i=1}^{n}y_{i}^{2}\le b^{2}. \end{array} \end{aligned}$$

Lemma 2

P(n) has an optimal solution if and only if $a^{2}\le nb^{2}$.

Proof

The objective function is continuous and the feasible region is compact, and so P(n) has an optimal solution if and only if the feasible region is nonempty. If $a^{2}\le nb^{2}$ then $y_{i}=\frac{a}{n}$ is feasible as it satisfies

$$\begin{aligned} \sum _{i=1}^{n}y_{i}^{2}=\sum _{i=1}^{n}\left( \frac{a}{n}\right) ^{2}=\frac{ a^{2}}{n}\le b^{2}. \end{aligned}$$

If $a^{2}>nb^{2}$ then the feasible region is empty, for any feasible y satisfies $\sum _{i=1}^{n}y_{i}=a$ so

$$\begin{aligned} \left( \sum _{i=1}^{n}y_{i}\right) ^{2}>nb^{2}\ge n\sum _{i=1}^{n}y_{i}^{2} \end{aligned}$$

But this contradicts

$$\begin{aligned} n\sum _{i=1}^{n}y_{i}^{2}\ge \left( \sum _{i=1}^{n}y_{i}\right) ^{2} \end{aligned}$$

which is true because

$$\begin{aligned} n\sum _{i=1}^{n}y_{i}^{2}-\left( \sum _{i=1}^{n}y_{i}\right) ^{2}= & {} n\sum _{i=1}^{n}\left( y_{i}-\frac{\sum _{i=1}^{n}y_{i}}{n}\right) ^{2} \\\ge & {} 0. \square \end{aligned}$$

Lemma 3

Suppose $a^{2}\le nb^{2}$. The optimal solution to P(n) is

$$\begin{aligned} y_{i}=\frac{a}{n}+\sqrt{nb^{2}-a^{2}}\frac{z_{i}-\bar{z}}{ns}. \end{aligned}$$

Proof

Since P(n) is a convex program we can solve it by minimizing the Lagrangian, giving

$$\begin{aligned} \min \nolimits _{y} -\sum _{i=1}^{n}z_{i}y_{i}+\mu \left( \sum _{i=1}^{n}y_{i}-a\right) +\lambda \left( \sum _{i=1}^{n}y_{i}^{2}-b^{2}\right) . \end{aligned}$$

Differentiating gives

$$\begin{aligned} -z_{i}+\mu +2\lambda y_{i}=0 \end{aligned}$$

so

$$\begin{aligned} y_{i}=\frac{z_{i}-\mu }{2\lambda } \end{aligned}$$

Since

$$\begin{aligned} \sum _{i=1}^{n}y_{i}=a \end{aligned}$$

we have

$$\begin{aligned} \sum _{i=1}^{n}z_{i}-n\mu =2\lambda a. \end{aligned}$$

If $a=0$ then

$$\begin{aligned} \mu =\bar{z} \end{aligned}$$

so

$$\begin{aligned} y_{i}=\frac{z_{i}-\bar{z}}{2\lambda } \end{aligned}$$

This gives

$$\begin{aligned} \sum _{i=1}^{n}y_{i}^{2}=\frac{ns^{2}}{4\lambda ^{2}}=b^{2} \end{aligned}$$

so

$$\begin{aligned} 2\lambda =\sqrt{n}\frac{s}{b} \end{aligned}$$

and

$$\begin{aligned} y_{i}=\frac{b\left( z_{i}-\bar{z}\right) }{s\sqrt{n}} \end{aligned}$$

as required.

If $a\ne 0$ then substituting for $\lambda$ gives

$$\begin{aligned} y_{i}=\frac{-z_{i}+\mu }{-\sum _{i=1}^{n}z_{i}+n\mu }a \end{aligned}$$

Now

$$\begin{aligned} \sum _{i=1}^{n}y_{i}^{2}=b^{2} \end{aligned}$$

so

$$\begin{aligned}&\sum _{i=1}^{n}\left( \frac{-z_{i}+\mu }{-\sum _{i=1}^{n}z_{i}+n\mu }\right) ^{2}a^{2}=b^{2}\nonumber \\&a^{2}\sum _{i=1}^{n}\left( -z_{i}+\mu \right) ^{2}=b^{2}\left( -\sum _{i=1}^{n}z_{i}+n\mu \right) ^{2}\nonumber \\&a^{2}\left( \sum _{i=1}^{n}z_{i}^{2}-2\mu \sum _{i=1}^{n}z_{i}+n\mu ^{2}\right) =b^{2}\left( \left( \sum _{i=1}^{n}z_{i}\right) ^{2}-2n\mu \sum _{i=1}^{n}z_{i}+n^{2}\mu ^{2}\right) \nonumber \\ \end{aligned}$$

(10)

Let $z=\sum _{i=1}^{n}z_{i}$, and $d=\sum _{i=1}^{n}z_{i}^{2}$, so

$$\begin{aligned} dn-z^{2}= & {} n\sum _{i=1}^{n}z_{i}^{2}-\left( \sum _{i=1}^{n}z_{i}\right) ^{2} \\= & {} n^{2}s^{2}. \end{aligned}$$

Equation (10) is a quadratic in $\mu$,

$$\begin{aligned} \left( a^{2}n-n^{2}b^{2}\right) \mu ^{2}-\left( 2a^{2}z-2znb^{2}\right) \mu +\left( a^{2}d-z^{2}b^{2}\right) =0, \end{aligned}$$

that has roots

$$\begin{aligned} \mu \in \left\{ \bar{z}-\frac{a}{\sqrt{nb^{2}-a^{2}}}s,\bar{z}+\frac{a}{\sqrt{ nb^{2}-a^{2}}}s\right\} . \end{aligned}$$

The first root is

$$\begin{aligned} \mu =\bar{z}-\frac{a}{\sqrt{nb^{2}-a^{2}}}s. \end{aligned}$$

This gives

$$\begin{aligned} y_{i}= & {} \frac{-z_{i}+\mu }{-\sum _{i=1}^{n}z_{i}+n\mu }a \\= & {} \frac{-z_{i}+\bar{z}-\frac{a}{\sqrt{nb^{2}-a^{2}}}s}{ -\sum _{i=1}^{n}z_{i}+n\bar{z}-\frac{na}{\sqrt{nb^{2}-a^{2}}}s}a \\= & {} \frac{a}{n}+\sqrt{nb^{2}-a^{2}}\frac{z_{i}-\bar{z}}{ns} \end{aligned}$$

with objective

$$\begin{aligned} -\sum _{i=1}^{n}z_{i}y_{i}= & {} -\sum _{i=1}^{n}z_{i}\frac{a}{n}-\sum _{i=1}^{n} \sqrt{nb^{2}-a^{2}}\frac{z_{i}^{2}-\bar{z}z_{i}}{ns} \\= & {} -a\bar{z}-s\sqrt{nb^{2}-a^{2}} \end{aligned}$$

The other root of (10) gives

$$\begin{aligned} y_{i}=\frac{a}{n}-\sqrt{nb^{2}-a^{2}}\frac{c_{i}-\bar{c}}{ns} \end{aligned}$$

which has value $-a\bar{z}+s\sqrt{nb^{2}-a^{2}}$, a local maximum of P(n).) $\square$

Lemma 6

Let $\bar{z}=\frac{\sum _{i=1}^{n}z_{i}}{n}$. For all i,

$$\begin{aligned} \left| z_{i}-\bar{z}\right| \le \sqrt{\frac{(n-1)}{n} \sum _{j=1}^{n}\left( z_{j}-\bar{z}\right) ^{2}}. \end{aligned}$$

Proof

Without loss of generality choose $i=1$. Then

$$\begin{aligned}&(n-1)\sum _{j=2}^{n}\left( z_{j}-\frac{1}{n}\sum _{i=1}^{n}z_{i}\right) ^{2}-\left( z_{1}-\frac{1}{n}\sum _{i=1}^{n}z_{i}\right) ^{2} \\&\quad =(n-1)\sum _{j=2}^{n}\left( z_{j}-\bar{z}\right) ^{2}-\left( z_{1}-\bar{z} \right) ^{2} \\&\quad =(n-1)\sum _{j=2}^{n}\left( z_{j}^{2}-2\bar{z}z_{j}+\bar{z}^{2}\right) -\left( z_{1}-\bar{z}\right) ^{2} \\&\quad =(n-1)\sum _{j=2}^{n}z_{j}^{2}-2(n-1)\bar{z}\sum _{j=2}^{n}z_{j}+(n-1)^{2} \bar{z}^{2}-z_{1}^{2}+2\bar{z}z_{1}-\bar{z}^{2} \\&\quad =(n-1)\sum _{j=2}^{n}z_{j}^{2}-2(n-1)\bar{z}(n\bar{z}-z_{1})+(n-1)^{2}\bar{z }^{2}-z_{1}^{2}+2\bar{z}z_{1}-\bar{z}^{2} \\&\quad =(n-1)\sum _{j=2}^{n}z_{j}^{2}-\left( z_{1}-\bar{z}n\right) ^{2} \\&\quad =(n-1)\sum _{j=2}^{n}z_{j}^{2}-\left( \sum _{j=2}^{n}z_{j}\right) ^{2}. \end{aligned}$$

This expression is $(n-1)^{2}$ times the variance of the quantities $z_{2},z_{3},\ldots ,z_{n}$ which is nonnegative, so it follows that

$$\begin{aligned} n\left( z_{1}-\frac{1}{n}\sum _{i=1}^{n}z_{i}\right) ^{2}\le (n-1)\left( \sum _{j=2}^{n}\left( z_{j}-\frac{1}{n}\sum _{i=1}^{n}z_{i}\right) ^{2}+\left( z_{1}-\frac{1}{n}\sum _{i=1}^{n}z_{i}\right) ^{2}\right) \end{aligned}$$

from which the result follows. $\square$

Lemma 4

If $r\le \sqrt{\frac{m}{m-1}}\min _{i}\{q_{i}\}$ then P has optimal solution

$$\begin{aligned} p_{i}=q_{i}+\frac{z_{i}-\bar{z}}{\sqrt{m}}\frac{r}{s}. \end{aligned}$$

with optimal value $\sum _{i=1}^{m}q_{i}z_{i}+\left( \sqrt{m}s\right) r$.

Proof

We consider a solution in which we drop the constraint $p\ge 0$. This gives $p_{i}=q_{i}+y_{i}$ where y solves

$$\begin{aligned} \begin{array}{cc} \min &{}\displaystyle -\sum _{i=1}^{m}z_{i}y_{i} \\ \text {s.t.} &{} \displaystyle \sum \nolimits _{i=1}^{m}y_{i}=0, \\ &{} \displaystyle \sum \nolimits _{i=1}^{m}y_{i}^{2}\le r^{2}. \end{array} \end{aligned}$$

Thus Lemma 3 with $a=0$ gives

$$\begin{aligned} y_{i}=\sqrt{mr^{2}}\frac{z_{i}-\bar{z}}{ms} \end{aligned}$$

whence

$$\begin{aligned} p_{i}=q_{i}+\frac{z_{i}-\bar{z}}{\sqrt{m}}\frac{r}{s}\text {.} \end{aligned}$$

Now by Lemma 6

$$\begin{aligned} \bar{z}-z_{i}\le & {} \sqrt{\frac{(m-1)}{m}\sum _{i}\left( z_{i}-\bar{z} \right) ^{2}} \\= & {} \sqrt{m-1}s \end{aligned}$$

so

$$\begin{aligned} p_{i}= & {} q_{i}+\frac{\bar{z}-z_{i}}{\sqrt{m}}\frac{r}{s} \\\ge & {} q_{i}+\frac{\sqrt{m-1}r}{\sqrt{m}} \\\ge & {} 0 \end{aligned}$$

as $r\le \sqrt{\frac{m}{m-1}}\min _{i}\{q_{i}\}$ by assumption.

We also have

$$\begin{aligned} \rho (Z(x))= & {} \sum _{i=1}^{m}p_{i}z_{i} \\= & {} \sum _{i=1}^{m}q_{i}z_{i}+\sum _{i=1}^{m}\frac{z_{i}-\bar{z}}{\sqrt{m}} \frac{r}{s}z_{i} \\= & {} \sum _{i=1}^{m}q_{i}z_{i}+\frac{r}{\sqrt{m}s}\sum _{i=1}^{m}\left( z_{i}^{2}-\bar{z}z_{i}\right) \\= & {} \sum _{i=1}^{m}q_{i}z_{i}+\frac{r}{\sqrt{m}s}\left( \sum _{i=1}^{m}z_{i}^{2}-m\bar{z}^{2}\right) \\= & {} \sum _{i=1}^{m}q_{i}z_{i}+\frac{r}{\sqrt{m}s}\left( ms^{2}\right) \\= & {} \sum _{i=1}^{m}q_{i}z_{i}+\left( \sqrt{m}s\right) r.\square \end{aligned}$$

Lemma 5

If $r_{1}<r_{2}$ then $k(r_{1})\le k(r_{2}).$

Proof

Recall for any fixed k and r that if $r\le \sqrt{\frac{1}{m(m-1)}}$,

$$\begin{aligned} p_{i}(r)=\frac{1}{m}+\frac{r}{s\sqrt{m}}(z_{i}-\bar{z}),\quad i=1,\ldots ,m, \end{aligned}$$

and otherwise

$$\begin{aligned} p_{i}(r)=\left\{ \begin{array}{ll} 0 &{}\quad i=1,\ldots ,k(r), \\ \frac{1}{(m-k(r))}+\frac{\sqrt{(m-k(r))r^{2}-\frac{k(r)}{m}}\frac{z_{i}-\bar{ z}}{s}}{(m-k(r))} &{}\quad i=k(r)+1,\ldots ,m. \end{array} \right. \end{aligned}$$

If we assume that k(r) is constant at value k as r varies then

$$\begin{aligned} \frac{dp_{i}(r)}{dr}=\left\{ \begin{array}{ll} 0 &{}\quad i=1,\ldots ,k, \\ \frac{z_{i}-\bar{z}}{s}\frac{r}{\sqrt{(m-k)r^{2}-\frac{k}{m}}} &{}\quad i=k+1,\ldots ,m. \end{array} \right. \end{aligned}$$

(11)

Since $z_{i}\le z_{i+1}$, for a fixed k and r we have

$$\begin{aligned} p_{i}(r)\le p_{i+1}(r),\quad i=k+1,\ldots ,m-1. \end{aligned}$$

Also $\frac{dp_{k+1}(r)}{dr}\le \frac{dp_{i}(r)}{dr}$ for $i=k+2,\ldots ,m$ and $\frac{dp_{k+1}(r)}{dr}\le 0$. This means that as r increases $p_{k+1}(r)$ becomes negative before (or at that same r as) any $p_{i}(r)$, $i=k+2,\ldots ,m$. In other words as r increases k(r) remains constant and then increases at some value of r. If the $z_{i}$ are distinct then k(r) increases by 1 at each step. We therefore have $k(0)=0$, and k(r) is piecewise constant and nondecreasing with r. $\square$

Proposition 1

Suppose that $Z(x,\omega )$ is a convex function of x for each $\omega \in \varOmega$, and that $g(\tilde{x},\omega )$ is a subgradient of $Z(x,\omega )$ at $\tilde{x}$. Then $\mathbb {E}_{\mathbb {P}^{*}}[g(\tilde{x},\omega )]$ is a subgradient of $\max _{\mathbb {P}\in \mathcal {P}}\mathbb {E}_{\mathbb {P} }[Z(x,\omega )]$ at $\tilde{x}$, where $\mathbb {P}^{*}\in \arg \max _{ \mathbb {P}\in \mathcal {P}}\mathbb {E}_{\mathbb {P}}[Z(\tilde{x},\omega _{m})]$.

Proof

For any x,

$$\begin{aligned} \max _{\mathbb {P}\in \mathcal {P}}\mathbb {E}_{\mathbb {P}}[Z(x,\omega )]= & {} \mathbb {E}_{\mathbb {P}^{*}}[Z(x,\omega )] \\\ge & {} \mathbb {E}_{\mathbb {P}^{*}}[Z(\tilde{x},\omega )+g(\tilde{x} ,\omega )^{\top }(x-\tilde{x})] \\= & {} \mathbb {E}_{\mathbb {P}^{*}}[Z(\tilde{x},\omega )]+\left( \mathbb {E}_{ \mathbb {P}^{*}}[g(\tilde{x},\omega )]\right) ^{\top }(x-\tilde{x}) \end{aligned}$$

which demonstrates that $\mathbb {E}_{\mathbb {P}^{*}}[g(\tilde{x},\omega )]$ is a subgradient at $\tilde{x}$. $\square$

Proposition 2

If for any $x_{t}\in \mathcal {X}_{t}(\omega _{t})$, $h_{t+1,k}-\bar{\pi } _{t+1,k}^{\top }H_{t+1}x_{t}\le \mathbb {E}_{\mathbb {P}_{t}^{*}}[Q_{t+1}(x_{t},\omega _{t+1})]$ for every $k=1,2,\ldots ,\nu$, then

$$\begin{aligned} \tilde{Q}_{t}(x_{t-1},\omega _{t})\le Q_{t}(x_{t-1},\omega _{t}). \end{aligned}$$

Proof

For any $x_{t}\in \mathcal {X}_{t}(\omega _{t})$ the optimal choice of $\theta _{t+1}$ satisfies

$$\begin{aligned} c_{t}^{\top }x_{t}+\theta _{t+1}= & {} c_{t}^{\top }x_{t}+\max _{k}\left\{ h_{t+1,k}- \bar{\pi }_{t+1,k}^{\top }H_{t+1}x_{t}\right\} \\\le & {} c_{t}^{\top }x_{t}+\mathbb {E}_{\mathbb {P}_{t}^{*}}[Q_{t+1}(x_{t},\omega _{t+1})] \end{aligned}$$

by hypothesis. It follows that

$$\begin{aligned} \tilde{Q}_{t}(x_{t-1},\omega _{t})= & {} \min _{x_{t}\in \mathcal {X}_{t}(\omega _{t})}\text { }\left\{ c_{t}^{\top }x_{t}+\max _{k}\{h_{t+1,k}-\bar{\pi } _{t+1,k}^{\top }H_{t+1}x_{t}\}\right\} \\\le & {} \min _{x_{t}\in \mathcal {X}_{t}(\omega _{t})}\text { }\left\{ c_{t}^{\top }x_{t}+\mathbb {E}_{\mathbb {P}_{t}^{*}}\left[ Q_{t+1}(x_{t},\omega _{t+1})\right] \right\} \\= & {} Q_{t}(x_{t-1},\omega _{t}) \end{aligned}$$

giving the desired result. $\square$

Appendix B

This appendix gives a formulation of a risk-neutral version of the hydrothermal scheduling model that we solve in Sect. 5. More details of this model can be found in Philpott and Pritchard (2013). The risk-neutral hydrothermal scheduling model we solve seeks a policy of electricity generation over T weeks that meets demand and minimizes the expected fuel cost of thermal generation. All data are deterministic except for weekly inflows that are assumed to be stagewise independent. This results in a large-scale stochastic programming model which is defined as follows. Let $\bar{x}_{j}$ denote the known storage in reservoir j at the start of the first week, and $x_{j}\left( t\right)$ denote the storage in reservoir j at the end of week t, and let $C_{T}(x)$ be the known expected future cost of the system at the end of week T when reservoir levels are x. The risk-neutral hydrothermal scheduling model can then be formulated as follows

$$\begin{aligned} \begin{array}{lll} \text {P(}\bar{x}\text {):} &{} \text {min} &{}\displaystyle {\mathbb {E}}\left[ \sum _{t=1}^{T}\sum _{i \in \mathcal {N}}\sum _{m\in \mathcal {F}(i)}\phi _{m}\sum _{b}T(b,t)f_{m}(b,t)+C_{T}(x(T))\right] \\ &{} \text {s.t.} &{} \displaystyle g_{i}(y(b,t))+\sum _{m\in \mathcal {F}(i)}f_{m}(b,t)+\sum _{m \in \mathcal {H}(i)}\gamma _{m}h_{m}(b,t)\\ &{}&{}\quad =D_{i}(b,t),\quad i\in \mathcal {N}, \\ &{} &{}\displaystyle x(t)=x(t-1)-S\sum _{b}T(b,t)(Ah(b,t)+As(b,t)-\omega (t)), \\ &{} &{} 0\le f_{m}(t)\le a_{m},\quad m\in \mathcal {F}(i), \quad i\in \mathcal {N}, \\ &{} &{} 0\le h_{m}(t)\le b_{m},\quad 0\le s_{m}(t)\le c_{m}, \quad m\in \mathcal {H}(i), \\ &{} &{} x(0)=\bar{x},\quad 0\le x_{j}(t)\le r_{j},\quad j\in \mathcal {J}, \quad i\in \mathcal {N},\quad y\in Y. \end{array} \end{aligned}$$

This description uses the following indices:

Index Refers to

t :

Index of week

i :

Node in transmission network

b :

Index of load block

m :

Index of plant

j :

Index of reservoir

$\mathcal {N}$ :

Set of nodes in transmission network

$\mathcal {F}(i)$ :

Set of thermal plants at node i

$\mathcal {H}(i)$ :

Set of hydro plants at node i

$\mathcal {J}$ :

Set of reservoirs

and the following notation for parameters and decision variables:

Symbol	Meaning	Units
$\phi _{m}$	Short-run marginal cost of thermal plant m	$/MWh
$f_{m}(b,t)$	Generation of thermal plant m in load block b in week t	MW
$x_{j}(t)$	Storage in reservoir j at end of week t	$\text {m}^{3}$
$\bar{x}_{j}$	Known storage in reservoir j at start of week t	$\text {m}^{3}$
h(b, t)	Vector of hydro releases in block b, week t	$\text {m}^{3}/\text {s}$
s(b, t)	Vector of hydro spills in block b, week t	$\text {m}^{3}/\text {s}$
$\omega (t)$	Inflow (assumed constant over the week)	$\text {m}^{3}/\text {s}$
$\gamma _{m}$	Conversion factor for water flow into energy	$\text {MWs}/\text {m}^{3}$
$D_{i}(b,t)$	Electricity demand in node i in block b, week t	MW
T(b, t)	Number of hours in load block b in week t	h
S	Number of seconds per hour (3600)
y(b, t)	Flow in transmission lines in load block b in week t	MW

$g_{i}(y)$	Sum of flow into node i when transmission flows are y	MW
$a_{m}$	Thermal plant capacity	MW
$b_{m}$	Hydro plant capacity	$\text {m}^{3}/\text {s}$
$c_{m}$	Spillway capacity	$\text {m}^{3}/\text {s}$
$r_{j}$	Reservoir capacity	$\text {m}^{3}$
Y	Feasible set of transmission flows
A	Node-arc incidence matrix of water network