Update: potential answer added.
I am trying to write a multiobjective pyomo optimization script that optimizes the battery capacity and operational schedule based on minimizing the cost. It is supposed to answer what is the optimal size of a battery that reduces the electricity bill by charging when prices are low and discharging when prices are low, e.g. minimizing grid purchases when prices are high. I have been able to optimize the operational schedule when the capacity of the battery is given. However, I am having issues implementing the capacity optimization. The hourly load profile (electricity consumption) for a building and hourly electricity purchase prices are given (load[i] and price[i]). The cost is a constant with price per kwh for battery capacity.
I am using the ipopt solver.
I have tried different approaches. The first one does provide an optimal schedule and a capacity. However, I am unsure if this is the correct approach for optimal capacity since it is implemented in the same objective function that minimizes cost of battery + cost of purchased electricity - cost of energy discharged from battery.
m = AbstractModel()
m.cap = Var(domain=Reals, bounds=(10,1000)) #battery capacity
m.pe_c = Var(time,domain=NonNegativeReals,bounds=(0,1000)) #charge to battery
m.pe_d = Var(time,domain=NegativeReals,bounds=(-1000,0)) #discharge from battery
m.soc = Var(soc_time,bounds=(value(min_soc),value(max_soc))) #state of charge between 20% and 95%
def obj_expression(m):
return ((m.cap *cost)+sum([price[i]*(load[i]) for i in time])) + sum([price[i]*(m.pe_c[i]) for i in time]) + sum([price[i]*(m.pe_d[i]) for i in time])
m.OBJ = Objective(rule=obj_expression, sense=minimize)
def soc_constraint_rule(m,i):
return m.soc[i+1] - m.soc[i] <= float(100) * dt * (m.pe_d[i] + m.pe_c[i]*energy_efficiency + self_discharge_power) / m.cap
m.soc_constraints = Constraint(time, rule=soc_constraint_rule)
def soc_start_rule(m): #battery starts at initial value
return m.soc[0] == soc_init
m.soc_start = Constraint(rule=soc_start_rule)
def soc_end_rule(m): #battery returns to initial value at end of simulation
return m.soc[n] == m.soc[0]
m.soc_end = Constraint(rule=soc_end_rule)
def discharge_rule(m,i): # Cannot discharge more than what battery can hold
return m.pe_d[i] <= m.cap *C_rate #C_rate is an efficiency constant, in this case 0.9
m.discharge_rule = Constraint(time, rule=discharge_rule)
def charge_rule(m,i): #cannot charge more than battery can hold
return m.pe_c[i] <= m.cap *C_rate
m.charge_rule = Constraint(time, rule=charge_rule)
instance = m.create_instance()
solver = SolverFactory('ipopt')
solver.solve(instance,tee=True)
My second approach is not working properly to optimize the operational schedule to minimize cost of purchased energy from the grid, but it does return both an operational schedule and a capacity. In this approach I try to use state-of-energy (SOE) instead of SOC, making the capacity a limit of SOE which is supposed to avoid the non-convexity in the SOC constraint that may cause ipopt problems converging. This one does not return m.utility values and I am generally unsure about the objective function being correct and if it avoids charging and discharging occurring at the same time which is not possible in a real-life scenario.
m = AbstractModel()
m.cap = Var(domain=Reals, bounds=(10,1000)) #battery capacity variable
m.utility = Var(time,domain=NonNegativeReals) #purchase from grid
m.pe_c = Var(time,domain=NonNegativeReals,bounds=(0,1000)) #charge to battery
m.pe_d = Var(time,domain=NegativeReals,bounds=(-1000,0)) #discharge from battery
m.SOE =Var(soe_time, domain=NonNegativeReals,bounds=(0,1000)) #state of energy inside battery
def obj_expression(m):
return (m.cap*cost) + sum([price[i]*(m.utility[i]) for i in time])
m.OBJ = Objective(rule=obj_expression, sense=minimize)
def battery_rule(m,i): #to connect the states of the amount of energy stored at, energy discharged from, energy charged to, the battery system
return m.SOE[i+1] == m.SOE[i] + m.pe_c[i] + m.pe_d[i]
m.battery_rule = Constraint(time, rule=battery_rule)
def capacity_rule(m,i): # a physical constraint that limits the amount of energy stored in the battery system
return m.SOE[i+1] <= m.cap
m.capacity_rule = Constraint(time, rule=capacity_rule)
def charging_rule(m,i): # a physical constraint that limits the amount of energy charged to the battery system
return m.pe_c[i] <= m.cap * C_rate
m.charging_rule = Constraint(time, rule=charging_rule)
def discharging_rule(m,i): # a physical constraint that limits the amount of energy discharged from the battery system
return m.pe_d[i] <= m.cap * C_rate
m.discharging_rule = Constraint(time, rule=discharging_rule)
def soe_start_rule(m):
return m.SOE[0] == soe_init
m.soe_start_rule = Constraint(rule=soe_start_rule)
def soe_end_rule(m):
return m.SOE[n] == m.SOE[0]
m.soe_end = Constraint(rule=soe_end_rule)
def demand_rule(m,i): #a physical constraint, representing the law of conservation of energy
return m.pe_d[i] + m.pe_c[i] + m.utility[i] == load[i]
m.demand_rule = Constraint(time, rule=demand_rule)
instance = m.create_instance()
solver = SolverFactory('ipopt')
solver.solve(instance,tee=True)
My final approach, which is not working at all is using the PAO library to access bilevel optimization. This could be the best approach because it clearly defines an upper and lower level objective, but I am having issues with constructing iterating constraints with a submodel. The code below returns "TypeError: 'NoneType' object is not iterable". However, I am sure there are larger logical errors in this approach due to human error. I have also tried to have capacity * cost as the upper-level objective, by swapping the variables m. with m.sub; making m.sub.cap into m.cap and m.utility into m.sub.utility, etc. but that returns "AttributeError: 'NoneType' object has no attribute 'transpose'" when i run the solver.
m = AbstractModel()
m.sub = pao.SubModel()
m.sub.cap = Var(domain=Reals, bounds=(10,1000)) #battery capacity variable
m.utility = Var(time,domain=NonNegativeReals) #purchase from grid
m.pe_c = Var(time,domain=NonNegativeReals,bounds=(0,1000)) #charge to battery
m.pe_d = Var(time,domain=NegativeReals,bounds=(-1000,0)) #discharge from battery
m.SOE =Var(soe_time, domain=NonNegativeReals,bounds=(0,1000)) #state of energy inside battery
def obj_expression(m):
return sum([price[i]*(m.pe_c[i]) for i in time]) +sum([price[i]*(m.pe_d[i]) for i in time])
m.OBJ = Objective(rule=obj_expression, sense=minimize)
m.sub.obj = Objective(expr=m.sub.cap*cost, sense=minimize)
def battery_rule(m,i): #to connect the states of the amount of energy stored at, energy discharged from, energy charged to, the battery system
return m.SOE[i+1] == m.SOE[i] + m.pe_c[i] + m.pe_d[i]
m.battery_rule = Constraint(time, rule=battery_rule)
def capacity_rule(m,i): # a physical constraint that limits the amount of energy stored in the battery system
return m.SOE[i+1] <= m.sub.cap
m.capacity_rule = Constraint(time, rule=capacity_rule)
def charging_rule(m,i): # a physical constraint that limits the amount of energy charged to the battery system
return m.pe_c[i] <= m.sub.cap * C_rate
m.charging_rule = Constraint(time, rule=charging_rule)
def discharging_rule(m,i): # a physical constraint that limits the amount of energy discharged from the battery system
return m.pe_d[i] <= m.sub.cap * C_rate
m.discharging_rule = Constraint(time, rule=discharging_rule)
def soe_start_rule(m):
return m.SOE[0] == soe_init
m.soe_start_rule = Constraint(rule=soe_start_rule)
def soe_end_rule(m):
return m.SOE[n] == m.SOE[0]
m.soe_end = Constraint(rule=soe_end_rule)
def demand_rule(m,i): #a physical constraint, representing the law of conservation of energy
return m.pe_d[i] + m.pe_c[i] + m.utility[i] == load[i]
m.demand_rule = Constraint(time, rule=demand_rule)
nlp = pao.Solver('ipopt')
opt = pao.Solver('pao.pyomo.REG', nlp_solver=nlp)
results = opt.solve(m)
I am writing this post in hopes that someone has done a similar multi-objective two-Stage stochastic optimization and are able to see my human errors or give me feedback on what approach is better. Please ask me if something is unclear as I can provide more details.
Update: I am playing around with two new constraints that may or may not help the above models. They are below
def cost_reduction(m,i): #making sure the without-battery cost of electricity is higher than the with-battery cost of electricity.
return (m.pe_c[i]*price[i]) + (m.pe_d[i]*price[i]) + (load[i]*price[i]) <= (load[i]*price[i])
m.cost_reduction = Constraint(time, rule=cost_reduction)
def battery_cost(m,i): #make sure recover battery cost M.pe_d is a negative number and m.year is bounded between 0,20.
return (m.cap*cost) + (m.pe_d[i]*price[i]) *m.year ==0 #I am aware that this is not the correct way but it shows the idea of a constraint that attempts to ensure a payback
m.battery_cost = Constraint(time, rule=battery_cost)
I believe I have solved it using gurobi solver. Had to make changes to the demand_rule and add a binary_charge_constraint
m = ConcreteModel()
m.cap = Var(domain=Reals, bounds=(10,1000))#battery capacity variable
m.utility = Var(time,domain=NonNegativeReals) #purchase from grid
m.pe_c = Var(time,domain=NonNegativeReals,bounds=(0,1000)) #charge to battery
m.pe_d = Var(time,domain=NonNegativeReals,bounds=(0,1000)) #discharge from battery
m.SOE =Var(soe_time, domain=NonNegativeReals,bounds=(0,1000)) #state of energy inside battery
def obj_expression(m):
return (m.cap*cost) + sum([price[i]*(m.utility[i]) for i in time])
m.OBJ = Objective(rule=obj_expression, sense=minimize)
def battery_rule(m,i): #to connect the states of the amount of energy stored at, energy discharged from, energy charged to, the battery system
return m.SOE[i+1] == m.SOE[i] + m.pe_c[i] - m.pe_d[i]
m.battery_rule = Constraint(time, rule=battery_rule)
def capacity_rule(m,i): # a physical constraint that limits the amount of energy stored in the battery system
return m.SOE[i] <= m.cap
m.capacity_rule = Constraint(time, rule=capacity_rule)
def charging_rule(m,i): # a physical constraint that limits the amount of energy charged to the battery system
return m.pe_c[i] <= m.cap * C_rate
m.charging_rule = Constraint(time, rule=charging_rule)
def discharging_rule(m,i): # a physical constraint that limits the amount of energy discharged from the battery system
return m.pe_d[i] <= m.cap * C_rate
m.discharging_rule = Constraint(time, rule=discharging_rule)
def soe_start_rule(m):
return m.SOE[0] == soe_init
m.soe_start_rule = Constraint(rule=soe_start_rule)
def soe_end_rule(m):
return m.SOE[n] == m.SOE[0]
m.soe_end = Constraint(rule=soe_end_rule)
def demand_rule(m,i): #a physical constraint, representing the law of conservation of energy
return m.pe_d[i] - m.pe_c[i] + m.utility[i] == load[i]
m.demand_rule = Constraint(time, rule=demand_rule)
def binary_charge_constraint(m,i): #making sure it is not charging while discharging and viceversa
return m.pe_d[i]*m.pe_c[i]== 0
m.binary_charge_constraint = Constraint(time, rule=binary_charge_constraint)
instance = m.create_instance()
opt = SolverFactory("gurobi")
opt.options['NonConvex'] = 2
results = opt.solve(instance, tee=True)