python – linear optimization – Stack Overflow

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


Leave a Comment