I am looking to simulate the hydropower generation in Norway, taking into account all the cascade dependencies (for an hourly resolution). To beginn with I have set up an optimisation problem in Pyomo for a set of three plants and three reservoirs for one day. For the moment I want to maximise profit over all hours and all plant/reservoirs calculated as: electricity price[$/MWh]*turbine output[MW]*3600[s=1h] + volume[m3]*energy equivalent[MWh/m3]*electricity price[$/MWh]
My constraints are:
Hourly load = generated power must equal load at each hour
Water balance = storage levels updated properly at each timestep
Start and end storage levels = at 60% of max storage a beginning and end of simulation period (this is in the water balance constraint)
QUESTIONS:
1) why does the print statement print out 0.0 for cascade_inflow and discharged_flow: because this is printed while the model is created and not while solving it?
2A) Termination condition is optimal and I have a value for the objective function, but number of solutions is 0: I think the problem lies in the constraint water balance, I will post the results further down for 6 hours. Do I need to somehow set upper and lower bound contraints for the water balance? edit: if I run the script from the command line with pyomo solve --solver=glpk script.py input.dat I get a solution displayed..?!
2B) The water constraint function does not behave properly. If I look at the results then the jump in volume from time-step 1 to 2 is not feasible. What is wrong there? Is the way I am adding the cascade flows incorrect, or does the variable m.volume just do what it wants?
3) Might it be a better idea to create a network flow problem? There are example codes for these kinds of problems in the Pyomo Gallery. But I am not sure yet how to model the nodes as reservoirs. (I will most probably make a new post for this as soon as I have tried to implement a script).
4) (This is my first post: anything I did wrong or should do better next time?)
Code (reading parameters left out)
from pyomo.environ import *
from pyomo.opt import SolverFactory
opt = SolverFactory("glpk")
# Initiate model instance
m = AbstractModel()
# Define the variable power for each time step
m.turbine = Var(m.stage, m.res, initialize=0, bounds=turbine_bounds, within=NonNegativeReals)
# Define the variable volume for each time step
m.volume = Var(m.stage, m.res, initialize=volume_Init, bounds=volume_bounds, within=NonNegativeReals)
# Define the variable spill flow for each time step
m.spilledFlow = Var(m.stage, m.res, initialize=0, bounds=spill_bounds, within=NonNegativeReals)
# Constrain total power generated over day
def hourly_load_rule(model, t):
return model.P_load*model.hourly_demand[t] <= sum(model.turbine[t, res] for res in model.res) <= model.P_load*model.hourly_demand[t]
m.hourly_load = Constraint(m.stage, rule=hourly_load_rule)
# Water balance equation
def water_balance_rule(model, t, r):
if t == model.T: # final volume is same as initial volume at 60%
return model.volume[t, r] == model.max_Vol[r]*0.6
elif t > 1:
cascade_inflow = 0
for stor in model.res:
# connectMat is a matrix that has a 1 where there is a connection and 0 where there is not
cascade_inflow = cascade_inflow + model.connectMat[stor, r]*(model.turbine[t, stor]/model.slope[stor]+model.spilledFlow[t, stor])
if model.connectMat[stor, r] == 1:
print(stor, r, t, value(cascade_inflow)) # this always prints out 0.0 for cascade_inflow
discharged_flow = model.turbine[t, r]/model.slope[r] # model.turbine is in MW: divide by slope to get discharge flow [m3/s]
print(value(discharged_flow)) # this always prints out 0.0
return model.volume[t, r] == model.volume[t-1, r]+(cascade_inflow+model.inflow[t, r]-model.spilledFlow[t, r]-discharged_flow)*model.secPerTimeStep
else:
# start volume constrained to 60% of max volume
return model.volume[t, r] == model.max_Vol[r]*0.6
m.water_balance = Constraint(m.stage, m.res, rule=water_balance_rule)
# Revenue = Objective function (sum over all hours and all plants/reservoirs)
def revenue_rule(model):
return sum(sum(model.el_price[i]*model.turbine[i, j]*model.secPerTimeStep+model.volume[i, j]*model.energy_stored[j]*model.el_price[i] for i in model.stage) for j in model.res)
m.obj = Objective(rule=revenue_rule, sense=maximize)
instance = m.create("three_input_stack.dat")
results = opt.solve(instance)
instance.display()
results.write()
Results
0.0
(1, 2, 2, 0.0)
0.0
(2, 3, 2, 0.0)
0.0
0.0
(1, 2, 3, 0.0)
0.0
(2, 3, 3, 0.0)
0.0
0.0
(1, 2, 4, 0.0)
0.0
(2, 3, 4, 0.0)
0.0
0.0
(1, 2, 5, 0.0)
0.0
(2, 3, 5, 0.0)
0.0
Model unknown
Variables:
turbine : Size=18, Index=turbine_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(1, 1) : 0 : 3.31 : 3.31 : False : False : NonNegativeReals
(1, 2) : 0 : 3.71 : 5.9 : False : False : NonNegativeReals
(1, 3) : 0 : 0.0 : 9.0 : False : False : NonNegativeReals
(2, 1) : 0 : 0.8 : 3.31 : False : False : NonNegativeReals
(2, 2) : 0 : 5.9 : 5.9 : False : False : NonNegativeReals
(2, 3) : 0 : 0.0 : 9.0 : False : False : NonNegativeReals
(3, 1) : 0 : 0.0 : 3.31 : False : False : NonNegativeReals
(3, 2) : 0 : 0.242202133966 : 5.9 : False : False : NonNegativeReals
(3, 3) : 0 : 6.31779786603 : 9.0 : False : False : NonNegativeReals
(4, 1) : 0 : 0.0 : 3.31 : False : False : NonNegativeReals
(4, 2) : 0 : 0.0 : 5.9 : False : False : NonNegativeReals
(4, 3) : 0 : 6.5 : 9.0 : False : False : NonNegativeReals
(5, 1) : 0 : 0.0 : 3.31 : False : False : NonNegativeReals
(5, 2) : 0 : 1.9665 : 5.9 : False : False : NonNegativeReals
(5, 3) : 0 : 4.6535 : 9.0 : False : False : NonNegativeReals
(6, 1) : 0 : 3.31 : 3.31 : False : False : NonNegativeReals
(6, 2) : 0 : 3.83 : 5.9 : False : False : NonNegativeReals
(6, 3) : 0 : 0.0 : 9.0 : False : False : NonNegativeReals
volume : Size=18, Index=volume_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(1, 1) : 0.0 : 39600000.0 : 66000000.0 : False : False : NonNegativeReals
(1, 2) : 0.0 : 10020000.0 : 16700000.0 : False : False : NonNegativeReals
(1, 3) : 0.0 : 1260000.0 : 2100000.0 : False : False : NonNegativeReals
(2, 1) : 0.0 : 32149783.0468 : 66000000.0 : False : False : NonNegativeReals
(2, 2) : 0.0 : 16684216.9532 : 16700000.0 : False : False : NonNegativeReals
(2, 3) : 0.0 : 2100000.0 : 2100000.0 : False : False : NonNegativeReals
(3, 1) : 0.0 : 32167783.0468 : 66000000.0 : False : False : NonNegativeReals
(3, 2) : 0.0 : 16700000.0 : 16700000.0 : False : False : NonNegativeReals
(3, 3) : 0.0 : 2100000.0 : 2100000.0 : False : False : NonNegativeReals
(4, 1) : 0.0 : 32185783.0468 : 66000000.0 : False : False : NonNegativeReals
(4, 2) : 0.0 : 16700000.0 : 16700000.0 : False : False : NonNegativeReals
(4, 3) : 0.0 : 2100000.0 : 2100000.0 : False : False : NonNegativeReals
(5, 1) : 0.0 : 32203783.0468 : 66000000.0 : False : False : NonNegativeReals
(5, 2) : 0.0 : 16700000.0 : 16700000.0 : False : False : NonNegativeReals
(5, 3) : 0.0 : 2100000.0 : 2100000.0 : False : False : NonNegativeReals
(6, 1) : 0.0 : 39600000.0 : 66000000.0 : False : False : NonNegativeReals
(6, 2) : 0.0 : 10020000.0 : 16700000.0 : False : False : NonNegativeReals
(6, 3) : 0.0 : 1260000.0 : 2100000.0 : False : False : NonNegativeReals
spilledFlow : Size=18, Index=spilledFlow_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(1, 1) : 0.0 : 0 : 10000.0 : False : True : NonNegativeReals
(1, 2) : 0.0 : 0 : 10000.0 : False : True : NonNegativeReals
(1, 3) : 0.0 : 0 : 10000.0 : False : True : NonNegativeReals
(2, 1) : 0.0 : 2069.67087236 : 10000.0 : False : False : NonNegativeReals
(2, 2) : 0.0 : 213.332062039 : 10000.0 : False : False : NonNegativeReals
(2, 3) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(3, 1) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(3, 2) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(3, 3) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(4, 1) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(4, 2) : 0.0 : 5.0 : 10000.0 : False : False : NonNegativeReals
(4, 3) : 0.0 : 4.22222222222 : 10000.0 : False : False : NonNegativeReals
(5, 1) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(5, 2) : 0.0 : 0.0 : 10000.0 : False : False : NonNegativeReals
(5, 3) : 0.0 : 5.86355555556 : 10000.0 : False : False : NonNegativeReals
(6, 1) : 0.0 : 0 : 10000.0 : False : True : NonNegativeReals
(6, 2) : 0.0 : 0 : 10000.0 : False : True : NonNegativeReals
(6, 3) : 0.0 : 0 : 10000.0 : False : True : NonNegativeReals
Objectives:
obj : Size=1, Index=None, Active=True
Key : Active : Value
None : True : 39250323.6272
Constraints:
hourly_load : Size=6
Key : Lower : Body : Upper
1 : 7.02 : 7.02 : 7.02
2 : 6.7 : 6.7 : 6.7
3 : 6.56 : 6.56 : 6.56
4 : 6.5 : 6.5 : 6.5
5 : 6.62 : 6.62 : 6.62
6 : 7.14 : 7.14 : 7.14
water_balance : Size=18
Key : Lower : Body : Upper
(1, 1) : 39600000.0 : 39600000.0 : 39600000.0
(1, 2) : 10020000.0 : 10020000.0 : 10020000.0
(1, 3) : 1260000.0 : 1260000.0 : 1260000.0
(2, 1) : 0.0 : 2.14204192162e-08 : 0.0
(2, 2) : 0.0 : -2.23517417908e-08 : 0.0
(2, 3) : 0.0 : -5.82076609135e-10 : 0.0
(3, 1) : 0.0 : 0.0 : 0.0
(3, 2) : 0.0 : 1.55250745593e-08 : 0.0
(3, 3) : 0.0 : -7.13669123797e-09 : 0.0
(4, 1) : 0.0 : 0.0 : 0.0
(4, 2) : 0.0 : -7.35411731512e-11 : 0.0
(4, 3) : 0.0 : -6.39488462184e-12 : 0.0
(5, 1) : 0.0 : 0.0 : 0.0
(5, 2) : 0.0 : -5.49960077478e-10 : 0.0
(5, 3) : 0.0 : -1.79056769412e-10 : 0.0
(6, 1) : 39600000.0 : 39600000.0 : 39600000.0
(6, 2) : 10020000.0 : 10020000.0 : 10020000.0
(6, 3) : 1260000.0 : 1260000.0 : 1260000.0
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 39250323.6272
Upper bound: 39250323.6272
Number of objectives: 1
Number of constraints: 25
Number of variables: 49
Number of nonzeros: 89
Sense: maximize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Error rc: 0
Time: 0.0750000476837
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
Input File
param secPerTimeStep:=3600;
param T:=6;
param numReservoirs:=3;
param connectMat:=
1 1 0
1 2 1
1 3 0
2 1 0
2 2 0
2 3 1
3 1 0
3 2 0
3 3 0;
param el_price:=
1 242.16
2 242.09
3 239.3
4 231.52
5 224.25
6 219.77;
param inflow:=
1 1 5
2 1 5
3 1 5
4 1 5
5 1 5
6 1 5
1 2 5
2 2 5
3 2 5
4 2 5
5 2 5
6 2 5
1 3 5
2 3 5
3 3 5
4 3 5
5 3 5
6 3 5;
param min_Vol:=
1 0.0
2 0.0
3 0.0;
param max_Vol:=
1 66000000.0
2 16700000.0
3 2100000.0;
param min_Turb_gen:=
1 0
2 0
3 0;
param max_Turb_gen:=
1 3.31
2 5.9
3 9.0;
param min_spill:=
1 0.0
2 0.0
3 0.0;
param max_spill:=
1 10000.0
2 10000.0
3 10000.0;
param min_discharge:=
1 0.0
2 0.0
3 0.0;
param max_discharge:=
1 20.0
2 15.0
3 8.0;
param slope:=
1 0.1655
2 0.3933
3 1.125;
param energy_stored:=
1 0.000046
2 0.000109
3 0.00031;
param hourly_demand:=
1 0.0351
2 0.0335
3 0.0328
4 0.0325
5 0.0331
6 0.0357;
param P_load:=200;
1) Zeros are always printed because m.tubine is initialized to 0. During model construction (which is when the print statements are executed) the expressions evaluate to 0.
2A) The result is automatically loaded back into the model after calling the solver, so the command instance.display() is printing the results. The results object is only used if you want to delay loading the solutions (which is an option you can set when calling the solver).
2B) Not sure. You can always run the command instance.water_balance.pprint() to confirm that the constraint expressions are what you expect them to be.