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?
Great.
So, two questions:
- 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. - 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?
-
@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 thedual.( -
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.
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.
- 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.
- 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?
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.
Maybe better to ask the questions step by step.
- in the link you shared, under the SDDP.evaluate section, is the outgoing reservoir volume set to 10 or increased by 10 units?
- 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?
- In the literature, such as one of your working papers, 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.
- 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)?
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.
1 Like
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, 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.