Lately, I’ve been working on a hydropower expansion problem using the SDDP method and the SDDP.jl tools.
Here are the most important things to understand about this problem:
The horizon plan is for 5 years.
At the beginning of each year, decisions about expansion are made (at the start of years 1, 2, 3, and 4). The expansion decisions are yearly decisions.
There’s a one-year delay, meaning if we increase power capacity at the beginning of year 1, it becomes usable from the second year on. That’s why we do not have expansion decision at the beginning of year 5.
The total available power capacity in, let’s say, year 4 is the sum of expansion in year 1, expansion in year 2, and expansion in year 3.
Alongside expansion choices, I also have many operational decisions to make, like how much power to generate and whether to satisfy demand, all of which are determined “every hour” within a year (in contrast to expansion decisions which are taken “yearly”).
The exogenous variable is demand. At the start of each year, we only know the demand for that year hour by hour. But we don’t know the demand for the years after, so we make a few scenarios for those.
I’ve figured out that my setup fits into a linear policy graph using SDDP.jl.
Now, I have a question: Should I create nodes just for the expansion years (so I’ll have only 5 or 4 nodes), or do I need to define nodes for every single hourly decision (5 * 365 nodes)?
Thanks.
Now, I have another question.
The water level in reservoirs are also state variables (SDDP.State) which I called water_level. In the first node (year) I initialized this state variable.
The point is that water_level.out is not calculated explicitly by water_level.in, but it should be calculated by the the operational decisions we get hourly within the year.
Then, whenever I try to calculate this water_level.out by these operational decisions I got the error:
Unable to retrieve solution from node 2.
Termination status : INFEASIBLE_OR_UNBOUNDED
Primal status : NO_SOLUTION
Dual status : NO_SOLUTION.
And once I relax this constraint (calculating water_level.out) the error is solved.
Is it possible that I get the error because I am not calculating water_level.out by water_level.in? Or should I search for the error in another place?
It’s hard to offer specific advice without a reproducible example of the code, but my hunch would be that your constraints are wrong. Most probably, you have bounds on water level, and the update constraints are trying to either force negative water level, or you haven’t accounted for spill, and the water level exceeds the maximum volume.
Thank you.
I checked that link.
So are you saying that logically I should be able to use SDDP.jl package and update water_level.out by operational decisions and not water_level.in state variable?
Thanks. The reason was exactly what you said. The relatively complete recourse assumption had been violated.
By fixing this issue, I solved the problem. The numerical stability report is as follows:
matrix range [4e-03, 2e+00]
objective range [1e+00, 1e+06]
bounds range [1e+04, 4e+04]
rhs range [2e-01, 3e+06]
Then I tried to play with the parameters.
Therefore, by just changing the coefficient of objective function (and not any constraints), the solver became stalled and failed to solve the problem, without showing any errors.
I should also mention that the numerical stability report did not change after changing the objective function parameters.
I know it is hard to diagnose the issue without having the detail of the code, but do you have any guess or suggestion?
Yes I am using Gurobi.
I have read the document. Also, I read the Gurobi’s page which is referred by the document.
Could you please tell me the suitable ranges?
In Gurobi’s page it is written that:
rhs: [1 1e+4]
Matrix coefficient: [1e-3 1e+6]
Bound range: [1e 1e+4]
Do you agree with those range? how about objective function coefficient range?
Yes. But you should see these as a general guide, not a fixed rule.
Perhaps as insight on why big ranges are a problem: SDDP uses a a cutting plane approximation to build policies. If you have an objective with a coefficient of 1e+00 on one variable, and a coefficient of 1e+06 on another, then the cutting planes will be 1,000,000 times steeper in one direction than another. This means that the “flatter” direction essentially doesn’t matter. Ideally, all your coefficients should be of a similar magnitude so that they contribute equally to the objective.
If you have numerical issues with Gurobi, it normally means that there’s something else going on that needs fixing, but it’s hard to guess what or offer suggestions without a reproducible example.
Okay. Here’s the problem. It is taking a long time to solve, not because of numerical problems, but because you have 3 node problems, each with 632576 variables. And you’re then going to need to repeatedly re-solve these problems again and again. What machine is this running on? It’ll need a lot of memory.
I guess from the 8760 that you’re trying to solve an hourly dispatch problem over a year.
To step back a bit:
From a modeling perspective, does it make sense to try and solve a multi-year problem with uncertainty at an hourly resolution? Can you claim with any certainty that what happens on 30 September at 02:00 in the morning in two years time is important?
In general, if you want to solve multi-year problems with uncertainty, you need to give up some fidelity at the short time-scales. That might mean daily/weekly dispatch problems, or it might mean using a subset of hours, appropriately weighted.
(I realize at the start of this thread that I suggested solving a year of operations within each node, but I didn’t really think through details.)
Thanks.
I see. I am running into this problem on my personal laptop, which has 16 GB of RAM. So, do you mean that by changing the coefficient of objective function, the problem got more complicated; therefore the code used almost all the memory, and it gets stalled?
Regarding the numerous operational nodes, you are correct. I am planning to implement sampling to reduce the size of operational computations.
Also could you tell me what “VariableRef in MOI.~” rows mean?
Thank you.
Actually I have three minor questions about SDDP package. I don’t know if I have to add new topic or I can ask here.
I will ask the questions, but please let me know if I have to ask somewhere else.
We have SDDP.train() command in the package. Does “iteration limit” argument specify the number of backward recursions we do? e.g if iteration limit is 3, we move backward on the entire of policy graph 3 times and add Kelly’s cuts?
After solving the problem, can I have access to the cutting planes the SDDP package generated? I ask because I see model.cuts.jason is created, however I cannot see anything related to cuts in that.
There is no way to get upper bound (in the minimization problem) per iteration using the SDDP package?
Hi dear Odow,
As you suggested, I changed the operational nodes’ time scales from “hour” to “day”.
Now, my model is running much faster.
However, when I use simulation (both in-sample and out-of sample simulation), the result does not seem optimal. I mean, e.g. for a 5-year problem, the model decides to use most of the initial water in the reservoirs in year 1 and 2, and generate more than demand to sell electricity with 10$/MWh fee. As a consequence, in year 3, 4, and 5 we don’t have enough water to generate electricity and satisfy the demand, then we have to pay much higher cost (100$/MWh) to purchase electricity.
I should mention that I defined the level of water in the reservoirs as state variables, so the model should be able to manage the level of the waters within the years.
I am completely confused. I even let the model run for an hour (300 iterations), so the lower bound (my model is to minimized) changes very very little in the last minutes of an hour (so we can say it converged). However, the optimal policy SDDP.jl gives me is very poor.
Can It be because of many operational nodes (365) that I have in a year?
Or my model has not been converged yet and I should let it run for more than an hour?
P.S: I have 6 hydro-reservoir plants + 11 run of river hydro plants. I use Gurobi to solve the problem, and set SDDP.jl to single cut.
