python – linear optimization PV-battery system

I just learned PYOMO and trying to optimize an off grid PV-battery 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()

Leave a Comment