Basic setup

Optimize the size of a reservoir over time under uncertainty.

Suppose we have the following inputs:

$q_t$: precipitation or inflow (volume / year)
$d_t$: water demand (volume / year)

Static parameters:

$e$: Evaporation or other losses (portion / year)
$T$: timestep (years)

Choice variables:

$w_t$: withdrawals (volume / year)
$x_t$: supersource water (volume / year)
$k_t$: investment in increasing capacity (volume / year)

State space variables:

$v_t$: reservoir level (volume)
$c_t$: reservoir capacity (volume)

Following reservoir dynamics: $$v_{t+1} = (1 - e)^T v_t + q_t - w_t$$ $$c_{t+1} = c_t + k_t$$ $v_0 = c_0 = 0$

Constrained such that,

$w_t + x_t \ge d_t$: Must satisfy demand
$v_t \le c_t$: Cannot over-fill reservoir
$w_t \le v_{t-1} + q_t$: Cannot over-withdrawal reservoir

And the objective is to minimize costs: $$\sum_t (A x_t + k_t) e^{-\delta t}$$

Handling of uncertainty

We have N Monte Carlo-like scenarios of the inputs, $q_{it}$ and $d_{it}$. In the example below, only precipitation will vary by scenario. Each scenario has a probability, $p_i$, s.t. $\sum p_i = 1$.

We estimate separate timeseries of the immediate choice variables, $w_{it}$ and $x_{it}$. These do not require long-term planning, and so do not necessitate a larger state-space. However, there is just one timeseries for investment choices, $k_t$.

The objective function now takes the expected value of these costs: $$\sum_t \left(A \left(\sum_i p_i x_{it}\right) + k_t\right) e^{-\delta t}$$

Handling of optimization under uncertainty

I propose that AWASH not do proper optimization under uncertainty. Under proper optimization, the state-space transition in each period is uncertain, which both means that optimization is taken under the knowledge that more information is provided in each period, and that the state-space needs to be enumerated.

Instead, I propose that an entire timeseries of long-term investments will be made in period $t_0$, not knowing which scenario will be entered. This investment will then be taken, the uncertainty resolved, and a new timeseries of investments will be optimized for in period $t_1$. This continues iteratively. Furthermore, it should be iteratively estimated for each of the N scenarios.

Under this plan, the number of variables only needs to multiply by the number of scenarios, and only for those variables that do not reflect long-term investments. The main loss is that the rational optimizer will design strategies that wait for uncertainty to be resolved, since it will never be resolved (in the future periods considered by the optimization, the optimizer still will not "know" which scenario it is in, as far as investment decisions are concerned). However, this is preferable to having scenario-specific investment decision variables, where the optimizer will always wait one more period to know which scenario it is in.

Single Optimization Example

There is uncertainty between two scenarios. Under scenario 1, there is always enough water (a demand of 0.5 under a supply of 1.0), while the other scenario goes into periods of drought. We optimize a single timeseries of investment, in time 0.



In [1]:

    
using JuMP
using Clp
using Gadfly
using DataFrames



In [2]:

    
time = 0:49
precip1 = repeat([1], inner=length(time))
precip2 = (cos.(2*pi*time / 25) + 1) / 3 + 1/3
demand = repeat([.5], inner=length(time));

Here is what the demand and supply look like.



In [3]:

    
plot(x=repmat(time, 3), y=[precip1; precip2; demand], color=repeat(["Precip 1", "Precip 2", "Demand"], inner=50), Geom.line)









    Out[3]:

We set up the model...



In [4]:

    
m = Model(solver=ClpSolver())

@variable(m, kk[time] >= 0) # Investment
@variable(m, xx1[time] >= 0) # Supersource under scenario 1
@variable(m, ww1[time] >= 0) # Withdrawals under scenario 1
@variable(m, xx2[time] >= 0) # Supersource under scenario 2
@variable(m, ww2[time] >= 0) # Withdrawals under scenario 2









    Out[4]:





$$ ww2_{i} \geq 0 \quad\forall i \in \{0,1,\dots,48,49\} $$



In [5]:

    
# Always satisfy demand
@constraint(m, ww1[time] + xx1[time] .>= demand[time + 1])
@constraint(m, ww2[time] + xx2[time] .>= demand[time + 1])









    Out[5]:





50-element Array{JuMP.ConstraintRef{JuMP.Model,JuMP.GenericRangeConstraint{JuMP.GenericAffExpr{Float64,JuMP.Variable}}},1}:
 ww2[0] + xx2[0] ≥ 0.5  
 ww2[1] + xx2[1] ≥ 0.5  
 ww2[2] + xx2[2] ≥ 0.5  
 ww2[3] + xx2[3] ≥ 0.5  
 ww2[4] + xx2[4] ≥ 0.5  
 ww2[5] + xx2[5] ≥ 0.5  
 ww2[6] + xx2[6] ≥ 0.5  
 ww2[7] + xx2[7] ≥ 0.5  
 ww2[8] + xx2[8] ≥ 0.5  
 ww2[9] + xx2[9] ≥ 0.5  
 ww2[10] + xx2[10] ≥ 0.5
 ww2[11] + xx2[11] ≥ 0.5
 ww2[12] + xx2[12] ≥ 0.5
 ⋮                      
 ww2[38] + xx2[38] ≥ 0.5
 ww2[39] + xx2[39] ≥ 0.5
 ww2[40] + xx2[40] ≥ 0.5
 ww2[41] + xx2[41] ≥ 0.5
 ww2[42] + xx2[42] ≥ 0.5
 ww2[43] + xx2[43] ≥ 0.5
 ww2[44] + xx2[44] ≥ 0.5
 ww2[45] + xx2[45] ≥ 0.5
 ww2[46] + xx2[46] ≥ 0.5
 ww2[47] + xx2[47] ≥ 0.5
 ww2[48] + xx2[48] ≥ 0.5
 ww2[49] + xx2[49] ≥ 0.5



In [6]:

    
for tt1 in time
    # Do not overfill reservoirs (v_t < V_t)
    @constraint(m, sum(precip1[tt+1] - ww1[tt] for tt in 0:tt1) <= sum(kk[tt] for tt in 0:tt1))
    @constraint(m, sum(precip2[tt+1] - ww2[tt] for tt in 0:tt1) <= sum(kk[tt] for tt in 0:tt1))

    # Do not withdraw more than available (w_t < v_t-1 + p_t)
    if (tt1 == 0)
        @constraint(m, ww1[tt1] <= precip1[tt1+1])
        @constraint(m, ww2[tt1] <= precip2[tt1+1])
    else
        @constraint(m, ww1[tt1] <= sum(precip1[tt+1] - ww1[tt] for tt in 0:(tt1-1)) + precip1[tt1+1])
        @constraint(m, ww2[tt1] <= sum(precip2[tt+1] - ww2[tt] for tt in 0:(tt1-1)) + precip2[tt1+1])
    end
end



In [7]:

    
@objective(m, Min, sum(((100xx1[tt] + 100xx2[tt])/2 + kk[tt]) * exp(-.05tt) for tt in time))









    Out[7]:





$$ 50 xx1_{0} + 50 xx2_{0} + kk_{0} + 47.5614712250357 xx1_{1} + 47.5614712250357 xx2_{1} + 0.951229424500714 kk_{1} + 45.241870901797974 xx1_{2} + 45.241870901797974 xx2_{2} + 0.9048374180359595 kk_{2} + 43.03539882125289 xx1_{3} + 43.03539882125289 xx2_{3} + 0.8607079764250578 kk_{3} + 40.936537653899094 xx1_{4} + 40.936537653899094 xx2_{4} + 0.8187307530779818 kk_{4} + 38.94003915357025 xx1_{5} + 38.94003915357025 xx2_{5} + 0.7788007830714049 kk_{5} + 37.04091103408589 xx1_{6} + 37.04091103408589 xx2_{6} + 0.7408182206817178 kk_{6} + 35.23440448593567 xx1_{7} + 35.23440448593567 xx2_{7} + 0.7046880897187134 kk_{7} + 33.51600230178197 xx1_{8} + 33.51600230178197 xx2_{8} + 0.6703200460356393 kk_{8} + 31.88140758108866 xx1_{9} + 31.88140758108866 xx2_{9} + 0.6376281516217732 kk_{9} + 30.326532985631673 xx1_{10} + 30.326532985631673 xx2_{10} + 0.6065306597126334 kk_{10} + 28.847490519024333 xx1_{11} + 28.847490519024333 xx2_{11} + 0.5769498103804866 kk_{11} + 27.44058180470132 xx1_{12} + 27.44058180470132 xx2_{12} + 0.5488116360940264 kk_{12} + 26.102288838050804 xx1_{13} + 26.102288838050804 xx2_{13} + 0.522045776761016 kk_{13} + 24.829265189570474 xx1_{14} + 24.829265189570474 xx2_{14} + 0.49658530379140947 kk_{14} + 23.618327637050733 xx1_{15} + 23.618327637050733 xx2_{15} + 0.4723665527410147 kk_{15} + 22.466448205861077 xx1_{16} + 22.466448205861077 xx2_{16} + 0.44932896411722156 kk_{16} + 21.370746597436334 xx1_{17} + 21.370746597436334 xx2_{17} + 0.42741493194872665 kk_{17} + 20.328482987029954 xx1_{18} + 20.328482987029954 xx2_{18} + 0.4065696597405991 kk_{18} + 19.337051172725058 xx1_{19} + 19.337051172725058 xx2_{19} + 0.3867410234545012 kk_{19} + 18.393972058572118 xx1_{20} + 18.393972058572118 xx2_{20} + 0.36787944117144233 kk_{20} + 17.496887455557765 xx1_{21} + 17.496887455557765 xx2_{21} + 0.3499377491111553 kk_{21} + 16.643554184903977 xx1_{22} + 16.643554184903977 xx2_{22} + 0.3328710836980795 kk_{22} + 15.831838468952657 xx1_{23} + 15.831838468952657 xx2_{23} + 0.31663676937905316 kk_{23} + 15.0597105956101 xx1_{24} + 15.0597105956101 xx2_{24} + 0.301194211912202 kk_{24} + 14.325239843009504 xx1_{25} + 14.325239843009504 xx2_{25} + 0.2865047968601901 kk_{25} + 13.62658965170063 xx1_{26} + 13.62658965170063 xx2_{26} + 0.2725317930340126 kk_{26} + 12.962013032294575 xx1_{27} + 12.962013032294575 xx2_{27} + 0.2592402606458915 kk_{27} + 12.329848197080322 xx1_{28} + 12.329848197080322 xx2_{28} + 0.24659696394160643 kk_{28} + 11.72851440468988 xx1_{29} + 11.72851440468988 xx2_{29} + 0.23457028809379762 kk_{29} + 11.156508007421492 xx1_{30} + 11.156508007421492 xx2_{30} + 0.22313016014842982 kk_{30} + 10.612398691337152 xx1_{31} + 10.612398691337152 xx2_{31} + 0.21224797382674304 kk_{31} + 10.094825899732768 xx1_{32} + 10.094825899732768 xx2_{32} + 0.20189651799465538 kk_{32} + 9.602495431037706 xx1_{33} + 9.602495431037706 xx2_{33} + 0.1920499086207541 kk_{33} + 9.134176202636732 xx1_{34} + 9.134176202636732 xx2_{34} + 0.18268352405273464 kk_{34} + 8.688697172522255 xx1_{35} + 8.688697172522255 xx2_{35} + 0.1737739434504451 kk_{35} + 8.264944411079327 xx1_{36} + 8.264944411079327 xx2_{36} + 0.16529888822158653 kk_{36} + 7.86185831568138 xx1_{37} + 7.86185831568138 xx2_{37} + 0.1572371663136276 kk_{37} + 7.478430961131751 xx1_{38} + 7.478430961131751 xx2_{38} + 0.14956861922263504 kk_{38} + 7.113703579325678 xx1_{39} + 7.113703579325678 xx2_{39} + 0.14227407158651356 kk_{39} + 6.766764161830635 xx1_{40} + 6.766764161830635 xx2_{40} + 0.1353352832366127 kk_{40} + 6.436745179390209 xx1_{41} + 6.436745179390209 xx2_{41} + 0.12873490358780418 kk_{41} + 6.122821412649095 xx1_{42} + 6.122821412649095 xx2_{42} + 0.1224564282529819 kk_{42} + 5.824207888674849 xx1_{43} + 5.824207888674849 xx2_{43} + 0.11648415777349697 kk_{43} + 5.540157918116694 xx1_{44} + 5.540157918116694 xx2_{44} + 0.11080315836233387 kk_{44} + 5.269961228093217 xx1_{45} + 5.269961228093217 xx2_{45} + 0.10539922456186433 kk_{45} + 5.012942186140186 xx1_{46} + 5.012942186140186 xx2_{46} + 0.1002588437228037 kk_{46} + 4.76845811077748 xx1_{47} + 4.76845811077748 xx2_{47} + 0.09536916221554961 kk_{47} + 4.535897664470624 xx1_{48} + 4.535897664470624 xx2_{48} + 0.09071795328941247 kk_{48} + 4.314679324968525 xx1_{49} + 4.314679324968525 xx2_{49} + 0.0862935864993705 kk_{49} $$



In [8]:

    
status = solve(m)
println(status)
println("Objective value: ", getobjectivevalue(m))









    



Optimal
Objective value: 0.7059239279866208

The investment quickly builds a reservoir to handle the risk of drought. Even under the high-precip scenario, the reservoir is filled up, but that excess is never used.



In [9]:

    
plot(x=repmat(time, 3), y=[getvalue(kk)[time]; getvalue(ww1)[time]; getvalue(ww2)[time]], color=repeat(["Investment", "Withdrawals 1", "Withdrawals 2"], inner=50), Geom.line)









    Out[9]:

Progressive Optimization

We do the entire optimization during each period, taking the size of the reservoir as given. In each period, the level of the reservoir and the actual precipitation come from scenario 2 (the droughty scenario).



In [10]:

    
kksaved = kk









    Out[10]:





$$ kk_{i} \geq 0 \quad\forall i \in \{0,1,\dots,48,49\} $$



In [11]:

    
volume = Float64[0]
invested = Float64[0]
for tt0 in time
    println(tt0)
    subtime = time[time .>= tt0]
    subm = Model(solver=ClpSolver())

    @variable(subm, kk[subtime] >= 0) # Investment
    @variable(subm, xx1[subtime] >= 0) # Supersource under scenario 1
    @variable(subm, ww1[subtime] >= 0) # Withdrawals under scenario 1
    @variable(subm, xx2[subtime] >= 0) # Supersource under scenario 2
    @variable(subm, ww2[subtime] >= 0) # Withdrawals under scenario 2
        
    for tt1 in subtime
        # Always satisfy demand
        @constraint(subm, ww1[tt1] + xx1[tt1] >= demand[tt1 + 1])
        @constraint(subm, ww2[tt1] + xx2[tt1] >= demand[tt1 + 1])

        # Do not overfill reservoirs (v_t < V_t)
        @constraint(subm, volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))
        @constraint(subm, volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))

        # Do not withdraw more than available (w_t < v_t-1 + p_t)
        if (tt1 == tt0)
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + precip2[tt1+1])
        else
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:(tt1-1)) + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:(tt1-1)) + precip2[tt1+1])
        end
    end
    
    @objective(subm, Min, sum(((100xx1[tt] + 100xx2[tt])/2 + kk[tt]) * exp(-.05tt) for tt in subtime))
    
    status = solve(subm)
    push!(invested, getvalue(kk)[tt0])
    push!(volume, volume[end] + precip2[tt0+1] - getvalue(ww2)[tt0])
end

The investment does not change (the single optimization and the iterative optimization lie exactly overtop each other). The is a relief: I had worried that without the committment of a single optimization, the optimizer might keep putting off building the reservoir.



In [12]:

    
plot(x=repmat(time, 2), y=[getvalue(kksaved)[time]; invested[2:51]], color=repeat(["Single", "Iterative"], inner=50), Geom.line)









    Out[12]:

Progressive Optimization with Scenario Updating

Finally, we allow Bayesian updating of the scenario we are believed to be in, as new information comes in. Below, we do this both for the case where the true scenario is scenario 1 and for scenario 2.

Note that the probabilities are strongly weighted toward scenario 1 (99% to 1%), while they weren't in previous cases. This is because drawing from the supersource is so expensive, so even a small weight on scenario 2 results in the same results as before. I changed the probabilities to keep things interesting.



In [13]:

    
volume = Float64[0]
invested = Float64[0]
prob1 = .99
prob2 = .01
for tt0 in time
    println(tt0)
    subtime = time[time .>= tt0]
    subm = Model(solver=ClpSolver())

    @variable(subm, kk[subtime] >= 0) # Investment
    @variable(subm, xx1[subtime] >= 0) # Supersource under scenario 1
    @variable(subm, ww1[subtime] >= 0) # Withdrawals under scenario 1
    @variable(subm, xx2[subtime] >= 0) # Supersource under scenario 2
    @variable(subm, ww2[subtime] >= 0) # Withdrawals under scenario 2
        
    for tt1 in subtime
        # Always satisfy demand
        @constraint(subm, ww1[tt1] + xx1[tt1] >= demand[tt1 + 1])
        @constraint(subm, ww2[tt1] + xx2[tt1] >= demand[tt1 + 1])

        # Do not overfill reservoirs (v_t < V_t)
        @constraint(subm, volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))
        @constraint(subm, volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))

        # Do not withdraw more than available (w_t < v_t-1 + p_t)
        if (tt1 == tt0)
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + precip2[tt1+1])
        else
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:(tt1-1)) + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:(tt1-1)) + precip2[tt1+1])
        end
    end
    
    @objective(subm, Min, sum(((prob1 * 100xx1[tt] + prob2 * 100xx2[tt]) + kk[tt]) * exp(-.05tt) for tt in subtime))
    
    status = solve(subm)
    push!(invested, getvalue(kk)[tt0])
    push!(volume, volume[end] + precip2[tt0+1] - getvalue(ww2)[tt0])
end

volume = Float64[0]
invested = Float64[0]
prob1 = .99
prob2 = .01
for tt0 in time
    println(tt0)
    subtime = time[time .>= tt0]
    subm = Model(solver=ClpSolver())

    @variable(subm, kk[subtime] >= 0) # Investment
    @variable(subm, xx1[subtime] >= 0) # Supersource under scenario 1
    @variable(subm, ww1[subtime] >= 0) # Withdrawals under scenario 1
    @variable(subm, xx2[subtime] >= 0) # Supersource under scenario 2
    @variable(subm, ww2[subtime] >= 0) # Withdrawals under scenario 2
        
    for tt1 in subtime
        # Always satisfy demand
        @constraint(subm, ww1[tt1] + xx1[tt1] >= demand[tt1 + 1])
        @constraint(subm, ww2[tt1] + xx2[tt1] >= demand[tt1 + 1])

        # Do not overfill reservoirs (v_t < V_t)
        @constraint(subm, volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))
        @constraint(subm, volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))

        # Do not withdraw more than available (w_t < v_t-1 + p_t)
        if (tt1 == tt0)
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + precip2[tt1+1])
        else
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:(tt1-1)) + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:(tt1-1)) + precip2[tt1+1])
        end
    end
    
    @objective(subm, Min, sum(((prob1 * 100xx1[tt] + prob2 * 100xx2[tt]) + kk[tt]) * exp(-.05tt) for tt in subtime))
    
    status = solve(subm)
    push!(invested, getvalue(kk)[tt0])
    push!(volume, volume[end] + precip2[tt0+1] - getvalue(ww2)[tt0])
    prob1 = .999 * prob1
    prob2 = 1 - prob1
end

kksaved = invested

invested2 = invested

volume = Float64[0]
invested = Float64[0]
prob1 = .99
prob2 = .01
for tt0 in time
    println(tt0)
    subtime = time[time .>= tt0]
    subm = Model(solver=ClpSolver())

    @variable(subm, kk[subtime] >= 0) # Investment
    @variable(subm, xx1[subtime] >= 0) # Supersource under scenario 1
    @variable(subm, ww1[subtime] >= 0) # Withdrawals under scenario 1
    @variable(subm, xx2[subtime] >= 0) # Supersource under scenario 2
    @variable(subm, ww2[subtime] >= 0) # Withdrawals under scenario 2
        
    for tt1 in subtime
        # Always satisfy demand
        @constraint(subm, ww1[tt1] + xx1[tt1] >= demand[tt1 + 1])
        @constraint(subm, ww2[tt1] + xx2[tt1] >= demand[tt1 + 1])

        # Do not overfill reservoirs (v_t < V_t)
        @constraint(subm, volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))
        @constraint(subm, volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:tt1) <= sum(invested) + sum(kk[tt] for tt in tt0:tt1))

        # Do not withdraw more than available (w_t < v_t-1 + p_t)
        if (tt1 == tt0)
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + precip2[tt1+1])
        else
            @constraint(subm, ww1[tt1] <= volume[tt0 + 1] + sum(precip1[tt+1] - ww1[tt] for tt in tt0:(tt1-1)) + precip1[tt1+1])
            @constraint(subm, ww2[tt1] <= volume[tt0 + 1] + sum(precip2[tt+1] - ww2[tt] for tt in tt0:(tt1-1)) + precip2[tt1+1])
        end
    end
    
    @objective(subm, Min, sum(((prob1 * 100xx1[tt] + prob2 * 100xx2[tt]) + kk[tt]) * exp(-.05tt) for tt in subtime))
    
    status = solve(subm)
    push!(invested, getvalue(kk)[tt0])
    push!(volume, volume[end] + precip1[tt0+1] - getvalue(ww1)[tt0])
    prob2 = .999 * prob2
    prob1 = 1 - prob1
end

invested1 = invested









    



0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49






    Out[13]:





51-element Array{Float64,1}:
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0198672
 0.084339 
 0.0      
 0.0      
 0.0      
 0.0      
 ⋮        
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0      
 0.0

As evidence mounts for scenario 1 (over-supply), the reservoir is first delayed, and then only partially built. When the true scenario is 2 (droughty), the reservoir is almost completely built initially (dispite the high probability against it), and then finished during the second drought.



In [16]:

    
plot(x=repmat(time, 3), y=[kksaved[2:51]; invested1[2:51]; invested2[2:51]], color=repeat(["Single", "Iterative 1", "Iterative 2"], inner=50), Geom.line)









    Out[16]: