This is just a basic Linear Programming problem I'm trying to solve to understand how pyomo works. I will add more constraints and turn it into a MILP later on. My index t represents one day in 0.25 timesteps. Using ConcreteModel(), and CPLEX as a solver.
The actual function for Dth(t) is not interesting for my problem it's only there for the sake of completeness.
model.t = RangeSet(0, 24, 0.25) #timesteps
model.Pth_gh = Var(model.t) #variable
Dth = {}
for i in model.i:
Dth[model.t[i]] = 2 + 4 * exp(-(model.t[i] - 9) ** 2 / 3 ** 2) + 2 * exp(-(model.t[i] - 17) ** 2 / 3 ** 2)
model.Dth = Param(model.t, initialize=Dth #actual parameter
Printing out these parameter/variable yields:
Pth_gh : Size=97, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
0.25 : None : None : None : False : True : Reals
0.5 : None : None : None : False : True : Reals
0.75 : None : None : None : False : True : Reals
1.0 : None : None : None : False : True : Reals
Dth : Size=97, Index=t, Domain=Any, Default=None, Mutable=False
Key : Value
0 : 2.0004936392163692
0.25 : 2.000808241156262
0.5 : 2.0013050898153586
0.75 : 2.002078298728981
1.0 : 2.0032639513411756
I don't know if you can do it easier and I feel like the timesteps could be pretty useful if I want to export my solutions later on. So I'm pretty happy with this part.
A generic constraint in words:
For each index t the said variable/combination of variables should have the value of the parameter on the right side.
This specific constraint in words:
Pth_gh should have the same value as Dth in each timestep.
Constraint in code:
def ThPowerBalance(model):
for t in model.t:
return model.Pth_gh[t] == model.Dth[t]
model.ThPowerBalanceEq = Constraint(model.t, rule=ThPowerBalance)
Printing the constraint yields the following:
Comparing this to my variable Dth, this pretty much is what I want to see, isn't it?
ThPowerBalanceEq : Size=97, Index=t, Active=True
Key : Lower : Body : Upper : Active
0 : 2.0004936392163692 : Pth_gh[0] : 2.0004936392163692 : True
0.25 : 2.000808241156262 : Pth_gh[0.25] : 2.000808241156262 : True
0.5 : 2.0013050898153586 : Pth_gh[0.5] : 2.0013050898153586 : True
0.75 : 2.002078298728981 : Pth_gh[0.75] : 2.002078298728981 : True
Now here's the catch. After trying to actually solve the problem the variables don't get the value assigned that was defined by my constraint. They stay empty. And that doesn't make sense to me.
Objective function
model.OBJ = Objective(expr=p_gas * dt * 1 / gh_eff * summation(model.Pth_gh)\
+ summation(model.pel, model.Pel_buy) * dt)
Note: Pel_buy is a variable. Analog definition to Pth_gh. pel is a parameter. Analog definition to Dth.
Solver settings
opt = SolverFactory('cplex')
opt.solve(model)
Printing Pth_gh after solving the model
Pth_gh : Size=97, Index=t
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : None : None : None : False : True : Reals
0.25 : None : None : None : False : True : Reals
0.5 : None : None : None : False : True : Reals
0.75 : None : None : None : False : True : Reals
1.0 : None : None : None : False : True : Reals
Note: Shouldn't the objective function just return a value? Instead it returns this:
4.11764705882353*(Pth_gh[0] + Pth_gh[0.25] + Pth_gh[0.5] + Pth_gh[0.75] + Pth_gh[1.0]...
All of that intro and you left out the code snippet that shows how you tried to print? :)
It takes a little getting used to how to extract values in pyomo. Here are a couple things to try that will get you the objective function value after solving.
# whatever solver you use above here ...
result = solver.solve(model)
print(result) # will show solver info, including the OBJ value
model.display() # will show you the values of ALL of the model variables and expressions, including the OBJ value
model.OBJ.display() # will show you the evaluation of the OBJ function
print(model.OBJ.expr()) # will print the evaluation of the OBJ function, and give direct access to the value
Try those. Comment back if you are still stuck!