# Marginal water value while using SDDP.jl

**URL:** <https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169>\
**Category:** Optimization (Mathematical)\
**Tags:** question, sddp\
**Created:** [March 27, 2024, 6:03am UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169 "2024-03-27T06:03:13Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [March 27, 2024, 6:03am UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/1 "2024-03-27T06:03:13Z")

</div>

I have a multi-stage stochastic hydro-electric scheduling problem that SDDP.jl models and solves.  
The model includes a set of constraints related to water level control.  
So now I’d like to determine the reservoirs’ marginal water value for each simulation.  
How can I do that? How can I calculate the shadow price of that set of constraints in SDDP.jl?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 27, 2024, 6:18am UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/2 "2024-03-27T06:18:02Z")

</div>

See

- [An introduction to SDDP.jl · SDDP.jl](https://sddp.dev/stable/tutorial/first_steps/#Custom-recorders)
- [An introduction to SDDP.jl · SDDP.jl](https://sddp.dev/stable/tutorial/first_steps/#Extracting-the-marginal-water-values)

---

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [March 27, 2024, 5:04pm UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/3 "2024-03-27T17:04:41Z")

</div>

Great.  
So, two questions:

1. How can we name a set of constraints like this:  
@constraint(model, i in 1:10, x[i] \<= y[i])  
I tried creating a name with index i, but it didn’t work.
2. Assume I define daily blocks for each weekly “node” in my model. Then, I wrote the water level constraints within each node on a daily basis. Can I still use the custom recorders for those constrains?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 27, 2024, 7:10pm UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/4 "2024-03-27T19:10:46Z")

</div>

1. `@constraint(model, my_con[i in 1:10], x[i] <= y[i])`, with the recorder `:price => (sp::JuMP.Model) -> JuMP.dual.(sp[:my_con]),` Note the `dual.(`

2. You can have a custom recorder record anything that is accessible in the model. You just need to use an appropriate data structure. If you get stuck, please post a reproducible example of the code you currently have.

---

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [April 1, 2024, 5:06pm UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/5 "2024-04-01T17:06:27Z")

</div>

Great!  
I started with a small problem. Everything made sense.  
Then I applied custom recorders to my case.  
There are two behaviors that do not appear rational.

1. Because a few blocks are included within nodes, I separated the water control constraints in each node into three sets.
  - Water balance between node input and next block
  - Water balance between blocks
  - Water balance between last block and node output.

The marginal values for the water balance across blocks within each are same for the majority of the them, which makes no sense to me.

1. The marginal values for one reservoir are in the 10,000 range, but at the node before the final node, they jump to 100,000.

I know it’s impossible to guess the reason without a reproducible code.  
Perhaps I misunderstand the concept of marginal value. Does this calculate how much the expected cost-to-go will decrease if we add one unit of water to the reservoir at specific node or blocks?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 1, 2024, 9:22pm UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/6 "2024-04-01T21:22:13Z")

</div>

> I know it’s impossible to guess the reason without a reproducible code.

Correct.

> Does this calculate how much the expected cost-to-go will decrease if we add one unit of water to the reservoir at specific node or blocks?

I don’t know what you calculated 😄

But yes, the outputs from `ValueFunction` are the change in expected cost to go with respect to the change in the outgoing level of the reservoir.

---

<div class="post-metadata">

**Author:** ![Sorooshsa](https://avatars.discourse-cdn.com/v4/letter/s/977dab/32.png) [@Sorooshsa](https://discourse.julialang.org/u/Sorooshsa)\
**Post date:** [April 2, 2024, 12:10am UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/7 "2024-04-02T00:10:32Z")

</div>

Maybe better to ask the questions step by step.

1. in the link you shared, under the SDDP.evaluate section, is the outgoing reservoir volume set to 10 or increased by 10 units?
2. Costum recorder shows what happens when we increase or decrease one unit on any constriant’s RHS. However, the SDDP.evaluate can only be used for testing changes in state variables. Is this true? So, I believe both custom recorders and SDDP.evaluate can determine the gradient of the cost-to-go function with respect to state variables. Am I wrong?
3. In the literature, [such as one of your working papers,](https://arxiv.org/pdf/2312.08556.pdf) the unit of water value is money per megawatt hour. However, if we use SDDP.evaluate for water level variables or a custom recorder for water level constraints, the unit of the output is currency/volume. So, if feasible, could you let me know how you calculated it.
4. In this reproducible code, I created a small network with two hydropowers, A and B. A has a reservoir. There are four nodes in the problem. Each node contains four blocks. After training, I played with different points in SDDP.evaluate, both positive and negative. In every situation, the cost of changing is negative. How can a shadow price be negative if we calculate the price at the negative point (which means decreasing the water level)?

```julia
Blocks = [1, 2, 3, 4]
Powerplant = ["A", "B"]
Inflow = [30 10 2 1; 15 15 15 15; 10 15 3 25; 5 18 12 3]
Demand = [63 52 20 31; 26 29 50 32; 42 53 29 22; 10 25 30 10]

OMEGA = Dict{Any, Any}()
for i in 1:4
    OMEGA[i] = [(v) for v in [0.7*Inflow[i,:], 1.0*Inflow[i,:], 1.3*Inflow[i,:]]]
end
P = fill(1/3, 3)

graph = SDDP.LinearGraph(4)
function builder(model, node)

    @variables(model, begin
        0 <= gen_hyd[i in Powerplant, j in Blocks]  
        0 <= ls[i in Blocks] <= 100
        0 <= wd[i in Powerplant, j in Blocks] # water discharge with production
        0 <= ww[i in Powerplant, j in Blocks] # water spillage - goes to the next power plant 
        0 <= pe[i in Blocks] <= 10
        water_inf[i in Powerplant, j in Blocks]
        1 <= ws[i in ["A"], j in Blocks] <= 15
    end)

    @variable(model, 1 <= ws_block[i in ["A"]] <= 15, SDDP.State, initial_value = 5)

    SDDP.parameterize(model, OMEGA[node], P) do ω
        for i in Powerplant
            for j in Blocks
                JuMP.fix(water_inf[i,j], ω[j])
            end
        end
    end

    @constraints(model, begin
        [i in Powerplant, j in Blocks], gen_hyd[i, j] <= 20
        [i in Powerplant, j in Blocks], gen_hyd[i, j] <= wd[i, j]
        [j in Blocks], sum(gen_hyd[i, j] for i in Powerplant) + ls[j] + pe[j] == Demand[node, j]
    end) # Generation -> operational node

    @constraints(model, begin
        cons_lvl1[i in ["A"]], ws[i, 2] == ws_block[i].in + (water_inf[i, 1] - ww[i, 1] - wd[i, 1])
        cons_lvl2[i in ["A"], j in Blocks[2:end - 1]], ws[i, j + 1] == ws[i, j] + (water_inf[i, j] - ww[i, j] - wd[i, j])
        cons_lvl3[i in ["A"]], ws_block[i].out == ws[i, end] + (water_inf[i, end] - ww[i, end] - wd[i, end])

        [i in ["B"], j in Blocks], 0 == water_inf[i, j] - ww[i, j] -wd[i, j] + (ww["A", j] + wd["A", j])

        [i in ["A"]], 1 <= ws_block[i].out <= 15
    end) # Hydropower balances
   
   
    @stageobjective(model, sum(sum(2*gen_hyd[i, j] for i in Powerplant)+ 20*pe[j] + 100*ls[j] for j in Blocks))
    return model
end
   
main_model = SDDP.PolicyGraph(
builder,
graph;
sense = :Min,
lower_bound = 0,
optimizer = Gurobi.Optimizer,
)
SDDP.numerical_stability_report(main_model)
SDDP.train(main_model; stopping_rules = [SDDP.TimeLimit(5)], cut_type = SDDP.MULTI_CUT)
V = SDDP.ValueFunction(main_model; node = 1)
SDDP.evaluate(V, Dict("ws_block[A]" => +10))
SDDP.evaluate(V, Dict("ws_block[A]" => -10))

```

Thanks for your time, and I’m sorry it’s too long.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 2, 2024, 12:24am UTC](https://discourse.julialang.org/t/marginal-water-value-while-using-sddp-jl/112169/8 "2024-04-02T00:24:54Z")

</div>

> in the link you shared, under the SDDP.evaluate section, is the outgoing reservoir volume set to 10 or increased by 10 units?

The outgoing reservoir volume is 10 units.

> Custum recorder shows what happens when we increase or decrease one unit on any constraint’s RHS. However, the SDDP.evaluate can only be used for testing changes in state variables. Is this true? So, I believe both custom recorders and SDDP.evaluate can determine the gradient of the cost-to-go function with respect to state variables. Am I wrong?

There’s a very subtle difference. The `SDDP.evaluate` with the `ValueFunction` records the gradient of the expected cost function. The custom recorders in `simulate` will record the gradient of the expected cost function _plus_ the current stage objective.

> In the literature, [such as one of your working papers,](https://arxiv.org/pdf/2312.08556.pdf) the unit of water value is money per megawatt hour. However, if we use SDDP.evaluate for water level variables or a custom recorder for water level constraints, the unit of the output is currency/volume. So, if feasible, could you let me know how you calculated it.

I just assumed that 1 unit of water was equal to 1 MWh of output. In reality, there’d be a production curve.

> In this reproducible code, I created a small network with two hydropowers, A and B. A has a reservoir. There are four nodes in the problem. Each node contains four blocks. After training, I played with different points in SDDP.evaluate, both positive and negative. In every situation, the cost of changing is negative. How can a shadow price be negative if we calculate the price at the negative point (which means decreasing the water level)?

> `SDDP.evaluate(V, Dict("ws_block[A]" => -10))`

This evaluates the value function at `-10`, not if the level was decreased by `-10`. In reality, the point is infeasible, but the value function doesn’t know that. It just looks at the future cost function.