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()