I just learned PYOMO and trying to optimize a system. However when l run the model l seem not to get a feasible solution. The results of my model is zero. I seem to have written the objective function and constrains the right way but l don’t know why the solution is showing zero values when run in pyomo.
import numpy as np
import pyomo.environ as pyo
def pyomoresults(mod): # define a function to print optimization results
""""Print results of the optimization of model 'mod' including
objective value at optimum, optimal decision variable values and
dual/slack information on constraints"""
print('')
print('=========================')
print('Solution')
print('=========================')
print('')
print('Objective Function Value at Optimum')
# loop over objective objects in the model
for f in mod.component_objects(pyo.Objective, active=True):
print ('Objective', f) # print the name of the objective
# loop over subobjectives (in this course we will only have one)
for index in f:
print (' ', index, pyo.value(f[index]))
print('')
print('Decision Variable Values at Optimum')
# loop over decision variable objects in model
for v in mod.component_objects(pyo.Var, active=True):
print ('Variable', v) # print variable name
for index in v: # loop over index of variable object
# print index/optimal values
print (' ', index, pyo.value(v[index]))
print('')
print ('Duals/Shadow Prices at Optimum')
# loop over constraint objects in model
for d in mod.component_objects(pyo.Constraint, active=True):
print ('Constraint',d) # print constraint name
for index in d: # loop over index of constraint object
# print dual value of constraint
print (' ', index, mod.dual[d[index]])
print('')
print ("Slack at Optimum")
# loop over constraint objects in model
for c in mod.component_objects(pyo.Constraint, active=True):
print (" Constraint",c) # print constraint name
for index in c: # loop over index of constraint object
# print constraint slack information
print (" ", index, c[index].lslack(), ' - ', c[index].uslack())
Load = load = np.array([0.580416667,0.539066667,0.390116667,0.232033333,
0.204533333,0.194716667,0.194633333,0.209233333,
0.247266668,0.407916668,0.537349998,0.576983332,
0.580216667,0.520566667,0.485200003,0.4197,
0.424300002,0.448333332,0.546983333,0.840733333,
1.320233332,0.856422014,0.921716667,0.720283335])
solar_rad = np.array([0,0,0,0,
0.846573268,6.670823882,22.34096457,48.40323145,
95.10129002,161.7686087,236.9894473,293.9150696,
305.3854497,294.6843366,251.7269744,182.2991627,
123.210826,73.11869927,33.55642336,9.910144956,
1.621109317,0.008980831,0,0]*5)
T= len(load)
# create a concrete model
model = pyo.ConcreteModel()
model.Batsize =pyo.Var(within=pyo.NonNegativeReals)
model.PV= pyo.Var(within=pyo.NonNegativeReals)
model.Bstate= pyo.Var(within=pyo.NonNegativeReals)
model.Pdischarge= pyo.Var(within=pyo.NonNegativeReals)
model.Pcharge_a=pyo.Var(within=pyo.NonNegativeReals)
model.Pcharge=pyo.Var(within=pyo.NonNegativeReals)
model.PV_gen =pyo.Var(within=pyo.NonNegativeReals)
model.solar_rad=pyo.Var(within=pyo.NonNegativeReals)
model.Pflow=pyo.Var(within=pyo.NonNegativeReals)
#Objective function
model.cost = pyo.Objective(expr= model.Batsize *1000 + model.PV* 1000)
#Constrains
def Charging_rule(model,A):
return model.Bstate[A] == model.Bstate[A-1] + model.Pdischarge[A] + model.Pcharge_a[A]
model.Charging_rule = pyo.Constraint(T, rule=Charging_rule)
# PV generation
def PV_generation(model,B):
return model.PV_gen[B] == model.PV*0.2*solar_rad[B]/1000
model.PV_generation = pyo.Constraint(T, rule=PV_generation)
# Pflow is the energy flow reuired *from the battery* to meet the load
# Negative if load greater than PV, positive if PV greater than load
def Pflow(model,C):
return model.Pflow[C] == model.PV_gen[C]-load[C]
model.Pflow = pyo.Constraint(T, rule=Pflow)
# Discharge cannot exceed available charge in battery
# Discharge is negative
def Discharge(model,D):
return model.Pdischarge[D] >= (-1)*model.Bstate[D-1]
model.Discharge = pyo.Constraint(T, rule=Discharge)
# Ensures that energy flow rePflowuired is satisifed by charge and discharge flows
def Energylaw(model,E):
return model.Pflow[E] == model.Pcharge[E] + model.Pdischarge[E]
model.Energylaw = pyo.Constraint(T, rule=Energylaw)
# If Pflow is negative (discharge), then it will at least ePflowual discharge rePflowuired
# If Pflow is positive (charge), then Pdischarge (discharge rePflowuired will ePflowual 0)
def Dischargelimit(model,F):
return model.Pdischarge[F] <= 0
return model.Pdischarge[F] <= Pflow[F]
model.Dischargelimit = pyo.Constraint(T, rule=Dischargelimit)
# Given the below, it will push Pflow available for charge to zero or to to or greater than excess PV
def chargelimit(model,G):
return model.Pcharge[G] >= 0
return model.Pdischarge[G] >= model.Pflow[G]
model.chargelimit = pyo.Constraint(T, rule=chargelimit)
def SOC(model,H):
return 0<=model.Bstate[H]<=model.Batsize
model.SOC=pyo.Constraint(T, rule=SOC)
# Limit amount charge delivered by the available space in the battery
def Chargespace(model,I):
return model.Pcharge_a[I] >= 0
return model.Pcharge_a[I] <= model.Pcharge[I]
return model.Pcharge_a[I] <= model.Batsize- model.Bstate[I-1]
model.SOC=pyo.Constraint(T, rule=Chargespace)
model.pprint()
model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
# set solver
opt = pyo.SolverFactory('glpk')
# solve optimization problem
results = opt.solve(model)
# display model results
model.display()