How to run this simple optimization program with Pyomo?

Question

I am trying to run a 'simple' optimation programm without any constraints to get the programming running. After getting the program running I am planning on adding constraints and other input variables. However, I cannot seem to get the model to work. I have filterd the used dateframe down to just 2 days for computational feasability. The dataframe named df contains dates on a 15 minute interval (datetime) and two column with building load and PV generation, both floats(See attached figure).

I used the following code just to see if the optimzation problem works. Note that no constraints are added (Pgrid_max is commented out). Eventually I will add a data set with data on EV charging to optimize bi-directional EV charging.

start_date = '2023-06-01 00:00:00'
end_date = '2023-06-02 23:59:59'

# Pgrid_max = (df['Building load [kWh]'].max() + 24*18*0.25) * 1.5

model = ConcreteModel()

# Time intervals
T = df[start_date:end_date].index.tolist() # Create a list of the timesteps

# Sets
model.T = Set(initialize=T) # Time interval

# Parameters
model.Pbuilding = Param(model.T, initialize=df[start_date:end_date]['Building load [kWh]'].to_dict())
model.Ppv = Param(model.T, initialize=df[start_date:end_date]['PV load [kWh]'].to_dict())


# Variables
model.Pgrid = Var(model.T, domain=Reals)
# model.Pev = Var(model.T, domain=PositiveReals, bounds=(0, 24*0.25)) 


# Energy balance equation
def energy_balance(model, t):
    return model.Pgrid[t] == model.Pbuilding[t] - model.Ppv[t]

# Objective function
def objective_rule(model):
    return sum(model.Pgrid[t] for t in model.T)

model.obj = Objective(rule=objective_rule, sense=minimize)

# Solve the model (assuming you have a solver like glpk installed)
solver = SolverFactory('glpk')
solver.solve(model, tee=True)

# Extract results directly into a pandas DataFrame
results = pd.DataFrame(index=T)
results['Pgrid'] = [model.Pgrid[t].value for t in T]

for t in model.T:
    results['Pgrid'][t] = model.Pgrid[t].value

print(results)

When running this code I keep getting the message that the problem has no solution. How can this be the case sinde Pgrid can take any real value. Printing the results gives only "None" values. Anyone that knows what I am doing wrong?

Answer 1

You have a couple issues in your code.

The most significant (and why you are probably getting zeros) is that you didn't construct the constraint. You need to call the function as a rule as shown in the code below.

A couple other tips/ideas:

• ALWAYS pprint the model to QA it. If you did, you would have noticed the missing constraint
• ALWAYS check to make sure the result is optimal or the results are junk. You can either inspect the solver's verbiage or put in an assertion as I show below
• When things get confusing, put away the pandas unless you are SUPER comfortable with it. Dictionaries as shown are easy to troubleshoot with small data
• If you are going to put things back into a DF when done, I'd suggest doing it like I show where you are using a dictionary to ensure the indexing is correct.

On your model:

• You probably want to make Pgrid a non-negative real and make the balance constraint greater than or equal (GTE) as shown to cover the cases where the pv > demand, unless it is the case that you can push power back to the grid.
• I see from your start/stop that you are treating the datetime things in your df as strings.... that is dangerous in many ways. If you have those type of things, ensure that the data type for that series (your index) is datetime and construct datetime objects for your cutoffs....go ALL IN on datetime or not at all.

Code:

"""
Simple Power Grid Model
"""
from pyomo.environ import *
import pandas as pd

### DATA

# take pandas out of the equation for simple model...
load = { 0:  8.2,
         15: 9.5,
         30: 7.7,
         45: 4.3}

pv = {   0:  4.2,
         15: 5.1,
         30: 6.0,
         45: 8.8}

# start_date = '2023-06-01 00:00:00'
# end_date = '2023-06-02 23:59:59'

model = ConcreteModel()

# Time intervals
# T = df[start_date:end_date].index.tolist() # Create a list of the timesteps

# Sets
model.T = Set(initialize=load.keys()) # Time interval

# Parameters
model.Pbuilding = Param(model.T, initialize=load) # df[start_date:end_date]['Building load [kWh]'].to_dict())
model.Ppv = Param(model.T, initialize=pv) # df[start_date:end_date]['PV load [kWh]'].to_dict())


# Variables
model.Pgrid = Var(model.T, domain=NonNegativeReals) # probably can't "push" back to grid...
# model.Pev = Var(model.T, domain=PositiveReals, bounds=(0, 24*0.25))


# Energy balance equation
def energy_balance(model, t):
    # make this a GTE constraint:  must pull enough from grid to make up deficit, if there is one
    return model.Pgrid[t] >= model.Pbuilding[t] - model.Ppv[t]
# you were MISSING the contraint construction call w/ a rule:
model.balance = Constraint(model.T, rule=energy_balance)

# Objective function
def objective_rule(model):
    return sum(model.Pgrid[t] for t in model.T)
model.obj = Objective(rule=objective_rule, sense=minimize)

#  INSPECT THE MODEL!!!
model.pprint()

# Solve the model (assuming you have a solver like glpk installed)
solver = SolverFactory('cbc')  # or glpk or ...
results = solver.solve(model, tee=True)  # catch the results object

#  Ensure that it was OPTIMAL or results are junk
status = results.Solver()['Termination condition'].value
assert status == 'optimal', f'error occurred, status: {status}.  Check model!'

# Extract results directly into a pandas DataFrame
results = pd.DataFrame(index=model.T)

# providing the item pairs is "safer" as there is no risk of issues
# with the order of the values pulled from the set model.T
results['Pgrid'] = pd.Series({k:v.value for k, v in model.Pgrid.items()}) # [model.Pgrid[t].value for t in model.T]

print(results)

Output:

1 Set Declarations
    T : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    4 : {0, 15, 30, 45}

2 Param Declarations
    Pbuilding : Size=4, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   8.2
         15 :   9.5
         30 :   7.7
         45 :   4.3
    Ppv : Size=4, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   4.2
         15 :   5.1
         30 :   6.0
         45 :   8.8

1 Var Declarations
    Pgrid : Size=4, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :  None : False :  True : NonNegativeReals
         15 :     0 :  None :  None : False :  True : NonNegativeReals
         30 :     0 :  None :  None : False :  True : NonNegativeReals
         45 :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : Pgrid[0] + Pgrid[15] + Pgrid[30] + Pgrid[45]

1 Constraint Declarations
    balance : Size=4, Index=T, Active=True
        Key : Lower              : Body      : Upper : Active
          0 :  3.999999999999999 :  Pgrid[0] :  +Inf :   True
         15 :                4.4 : Pgrid[15] :  +Inf :   True
         30 : 1.7000000000000002 : Pgrid[30] :  +Inf :   True
         45 : -4.500000000000001 : Pgrid[45] :  +Inf :   True

6 Declarations: T Pbuilding Ppv Pgrid balance obj
Welcome to the CBC MILP Solver 
Version: 2.10.10 
Build Date: Jul 24 2023 

command line - /opt/homebrew/opt/cbc/bin/cbc -printingOptions all -import /var/folders/7l/f196n6c974x3yjx5s37t69dc0000gn/T/tmp28wkyeg0.pyomo.lp -stat=1 -solve -solu /var/folders/7l/f196n6c974x3yjx5s37t69dc0000gn/T/tmp28wkyeg0.pyomo.soln (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 0 (-4) rows, 0 (-4) columns and 0 (-4) elements
Statistics for presolved model


Problem has 0 rows, 0 columns (0 with objective) and 0 elements
Column breakdown:
0 of type 0.0->inf, 0 of type 0.0->up, 0 of type lo->inf, 
0 of type lo->up, 0 of type free, 0 of type fixed, 
0 of type -inf->0.0, 0 of type -inf->up, 0 of type 0.0->1.0 
Row breakdown:
0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, 
0 of type E other, 0 of type G 0.0, 0 of type G 1.0, 
0 of type G other, 0 of type L 0.0, 0 of type L 1.0, 
0 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, 
0 of type Free 
Presolve 0 (-4) rows, 0 (-4) columns and 0 (-4) elements
Empty problem - 0 rows, 0 columns and 0 elements
Optimal - objective value 10.1
After Postsolve, objective 10.1, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 10.1 - 0 iterations time 0.002, Presolve 0.00
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

    Pgrid
0     4.0
15    4.4
30    1.7
45    0.0

Process finished with exit code 0