Search code examples
pythonpandasoptimizationpyomo

Pyomo battery optimisation. If condition

I am currently building a simple optimisation tool for battery storage. At time t, the battery is loaded by Q_charge[t], which is the electricity drawn from the electric grid Q_el[t] multiplied by an efficiency factor eta_charge.

  1. Q_el[t]*eta_charge=Q_charge[t]

The discharge power Q_discharge[t] is multiplied by an efficiency factor eta_discharge and gives the demand Q_demand[t]. Demand profile is known for every hour of the year.

  1. Q_discharge[t]*eta_discharge=Q_demand[t]

The battery SOC(t) is defined as state of charge at previous time step + power input - power output. No losses considered.

  1. SOC(t)=SOC(t-1)+Q_charge[t]-Q_demand[t]

The state of charge at t=0 should be equal to the last hour of the period considered.

  1. SOC(t=0)=SOC(t=end)

Charge and discharge cannot occurr simultaneously.

  1. Charge[t]+Discharge[t] <= 1

Charge and discharge have min and max limits.

Knowing the price of electricity Price[t] for every hour of the year, the objective function will decide when to charge, minimizing the cost of purchase. The demand is used internally, no electricity sold.

Demand[t] is known (imported from Excel). Price[t] is known (imported from Excel).

Here the code.

model = pyo.ConcreteModel()

##Parameters

tot_time=8760 #Hours in a year
model.T=pyo.Set(initialize=[t for t in range(0,tot_time)])
model.EtaCharge=pyo.Param(initialize=1.0, mutable=True) #-    
model.PMinCharge=0 #MW
model.PMaxCharge=2 #MW
model.EtaDischarge=pyo.Param(initialize=1.0, mutable=True) #-
model.PMinDischarge=2 #MW
model.PMaxDischarge=4 #MW
model.S0=pyo.Param(initialize=0.0, mutable=True) #MWh SOC at t=0

##Variables

model.QEl = pyo.Var(model.T, within=NonNegativeReals, initialize=0)
model.QCharge= pyo.Var(model.T, within=NonNegativeReals)
model.QDischarge= pyo.Var(model.T, within=NonNegativeReals)
model.Charge=pyo.Var(model.T, within = pyo.Binary)
model.Discharge=pyo.Var(model.T, within = pyo.Binary)
model.SOC=pyo.Var(model.T, bounds=(0,100)) #Size of the battery limited to 100 MWh

##Constraints

model.C1 = pyo.ConstraintList()
model.C2 = pyo.ConstraintList()
model.C3 = pyo.ConstraintList()
model.C4 = pyo.ConstraintList()
model.C5 = pyo.ConstraintList()
model.C6 = pyo.ConstraintList()
model.C7 = pyo.ConstraintList()

for t in model.T:
    
    model.C1.add(model.QEl[t] * model.EtaCharge ==  model.QCharge[t])   

    model.C2.add(model.QDischarge[t] * model.EtaDischarge ==  Demand[t])   
          
    model.C3.add(model.Charge[t] + model.Discharge[t] <= 1)
    
    model.C4.add(model.QCharge[t] >= model.PMinCharge*model.Charge[t])
                 
    model.C5.add(model.QCharge[t] <= model.PMaxCharge*model.Charge[t])
    
    model.C6.add(model.QDischarge[t] >= model.PMinDischarge*model.Discharge[t])
    
    model.C7.add(model.QDischarge[t] <= model.PMaxDischarge*model.Discharge[t])  


def SOC_storage(model,t):

    if t == 0:
        return (model.SOC[t] == model.S0+ model.QCharge[t]-model.QDischarge[t]) 
    else:
        return (model.SOC[t] == model.SOC[t-1] + model.QCharge[t] - model.QDischarge[t])

model.C8 = pyo.Constraint(model.T,rule = SOC_storage)

model.C9 = pyo.Constraint(expr = model.S0 == model.SOC[len(tot_time)-1])

##Objective

model.objective = pyo.Objective(expr = sum(model.QEl[t]*Price[t] for t in model.T), sense = pyo.minimize)

opt=SolverFactory('glpk')
results = opt.solve(model, tee = True)
model.write("myfile_lp.lp", io_options={'symbolic_solver_labels':True})
results.solver.status
results.solver.termination_condition
model.pprint()
value(model.objective)

The code, as it is, runs smoothly. However, how can I insert the condition to NOT buy electricity (Charge[t]=0) only during specific hours? Let's say the battery cannot be charged between 08.00 and 15.00 every day? I tried with:

day_stop=np.arange(8,15,1) # Index of hours When we don't charge
week_stop= np.concatenate([day_stop+24*j for j in range(0,7)]) # Index of hours When we don't charge for the first week
year_stop = np.concatenate([week_stop+168*j for j in range(0,52)]) # Index of hours When we don't charge for the all year

def NoCharge(model,t):
    if t in year_stop :
        return model.Charge[t]==0
    else:
        return Constraint.Skip

model.C10 = pyo.Constraint(model.T,rule = NoCharge)

All the variables after running the code are equal to zero or None. Am I formulating the constraint in a wrong way?

Any idea how to formulate the latest constraint in a better way? Thanks so much for your help!

Solution

  • Your approach with the NoCharge function appears to be correct. It is a bit difficult to troubleshoot without the data. So... I made some up.

    Some recommendations:

    • It is wayyyy difficult to troubleshoot a model with 8000+ time periods, so I strongly recommend you chop it down (like I did below) to something very small so you can pprint and display the model and see what is going on and x-check it
    • Your model is clearly written, which is nice!
    • Your model is very easy to make infeasible because of the charge/discharge limits and the limit that the end state==start state, so be careful and you must check the solver status before assessing your model. I'd bet $5 that when you added in the constraint that the model was INFEASIBLE. So, you may want to relax the charge/discharge limits while working that out

    Edit... A couple other things occurred to me... Get it working first, then consider a couple of these...

    • You can probably eliminate one of your charge/discharge binary variables because they are mutually exclusive (enforced by your constraint). So charge = 1 - discharge is true and you allow that charge could equal zero, so "idle" is included. You could just include 1 of the 2 binaries and use the "1 minus" construct where needed.
    • You may want to include another power source and just cost it for situations where the demand > what can be produced and your battery can't keep up. It would be easy to add an "overhead" power variable and add it to your OBJ with some small penalty. This would really help keep your model feasible by making it more elastic.
    • If your model becomes slow to execute, you should probably increase your time steps to 2, 4, ... hours or just bin the "night hours" into a large chunk to reduce size

    Anyhow, the below works and produces a result. I used your C10 constraint and you can toggle it and C11 on/off to see that it is working in the output. C11 is redundant with C10 and just shows an alternate construct... You can use either.

    CODE:

    import pyomo.environ as pyo

model = pyo.ConcreteModel()

##Parameters

tot_time=10 #Hours in a year
Demand = [0, 2, 3, 0, 0, 3, 3, 0, 0, 3]
Price = [0.5 for t in range(tot_time)]

model.T=pyo.Set(initialize=[t for t in range(0,tot_time)])
model.EtaCharge=pyo.Param(initialize=1.0, mutable=True) #-    
model.PMinCharge=0 #MW
model.PMaxCharge=5 #MW
model.EtaDischarge=pyo.Param(initialize=1.0, mutable=True) #-
model.PMinDischarge=2 #MW
model.PMaxDischarge=4 #MW

model.S0=pyo.Param(initialize=10, mutable=True) #MWh SOC at t=0

##Variables

model.QEl = pyo.Var(model.T, within=pyo.NonNegativeReals, initialize=0)
model.QCharge= pyo.Var(model.T, within=pyo.NonNegativeReals)
model.QDischarge= pyo.Var(model.T, within=pyo.NonNegativeReals)
model.Charge=pyo.Var(model.T, within = pyo.Binary)
model.Discharge=pyo.Var(model.T, within = pyo.Binary)
model.SOC=pyo.Var(model.T, bounds=(0,100)) #Size of the battery limited to 100 MWh

##Constraints

model.C1 = pyo.ConstraintList()
model.C2 = pyo.ConstraintList()
model.C3 = pyo.ConstraintList()
model.C4 = pyo.ConstraintList()
model.C5 = pyo.ConstraintList()
model.C6 = pyo.ConstraintList()
model.C7 = pyo.ConstraintList()

for t in model.T:
    
    model.C1.add(model.QEl[t] * model.EtaCharge == model.QCharge[t])   

    model.C2.add(model.QDischarge[t] * model.EtaDischarge == Demand[t])   
          
    model.C3.add(model.Charge[t] + model.Discharge[t] <= 1)
    
    model.C4.add(model.QCharge[t] >= model.PMinCharge*model.Charge[t])
                 
    model.C5.add(model.QCharge[t] <= model.PMaxCharge*model.Charge[t])
    
    model.C6.add(model.QDischarge[t] >= model.PMinDischarge*model.Discharge[t])
    
    model.C7.add(model.QDischarge[t] <= model.PMaxDischarge*model.Discharge[t])  


def SOC_storage(model,t):

    if t == 0:
        return (model.SOC[t] == model.S0 + model.QCharge[t] - model.QDischarge[t]) 
    else:
        return (model.SOC[t] == model.SOC[t-1] + model.QCharge[t] - model.QDischarge[t])

model.C8 = pyo.Constraint(model.T,rule = SOC_storage)

model.C9 = pyo.Constraint(expr = model.S0 == model.SOC[tot_time-1])

def NoCharge(model,t):
    if t in no_charge_zone :
        return model.Charge[t]==0
    else:
        return pyo.Constraint.Skip

no_charge_zone = [7, 8]
model.C10 = pyo.Constraint(model.T,rule = NoCharge)

# The below C11 is redundant with C10...  It just shows a different approach.
# I think this is cleaner....  Just make a subset from your "prohibited" time periods
# and pass that in...
def NoCharge2(model, t):
    return model.Charge[t] == 0

model.C11 = pyo.Constraint(no_charge_zone, rule=NoCharge2)

##Objective

model.objective = pyo.Objective(expr = sum(model.QEl[t]*Price[t] for t in model.T), sense = pyo.minimize)

opt=pyo.SolverFactory('glpk')
results = opt.solve(model, tee = True)
#model.write("myfile_lp.lp", io_options={'symbolic_solver_labels':True})
results.solver.status
results.solver.termination_condition
#model.pprint()
model.Charge.display()
model.SOC.display()
#print(pyo.value(model.objective))

    Output:

    INTEGER OPTIMAL SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.1 Mb (87515 bytes)
Writing MIP solution to '/var/folders/rt/5gc2md7s7y5bw8ddt74tz5vw0000gp/T/tmpi82nq7rg.glpk.raw'...
154 lines were written
Charge : Size=10, Index=T
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      0 :     0 :   1.0 :     1 : False : False : Binary
      1 :     0 :   0.0 :     1 : False : False : Binary
      2 :     0 :   0.0 :     1 : False : False : Binary
      3 :     0 :   1.0 :     1 : False : False : Binary
      4 :     0 :   1.0 :     1 : False : False : Binary
      5 :     0 :   0.0 :     1 : False : False : Binary
      6 :     0 :   0.0 :     1 : False : False : Binary
      7 :     0 :   0.0 :     1 : False : False : Binary
      8 :     0 :   0.0 :     1 : False : False : Binary
      9 :     0 :   0.0 :     1 : False : False : Binary
SOC : Size=10, Index=T
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      0 :     0 :  14.0 :   100 : False : False :  Reals
      1 :     0 :  12.0 :   100 : False : False :  Reals
      2 :     0 :   9.0 :   100 : False : False :  Reals
      3 :     0 :  14.0 :   100 : False : False :  Reals
      4 :     0 :  19.0 :   100 : False : False :  Reals
      5 :     0 :  16.0 :   100 : False : False :  Reals
      6 :     0 :  13.0 :   100 : False : False :  Reals
      7 :     0 :  13.0 :   100 : False : False :  Reals
      8 :     0 :  13.0 :   100 : False : False :  Reals
      9 :     0 :  10.0 :   100 : False : False :  Reals