Search code examples
pythonloopsoptimizationpyomo

Pyomo:For loop in Pyomo constraints


I am struggling on making constraints using for loop in Python Pyomo. My code concept is the following:

model.Q = Set()
model.sbase=Param(model.Q, within = PositiveReals)

My code can run properly with one Sbase. But if I get 3 different numbers of Sbase. And what I want to do is to output 3 different image by using these 3 numbers.

set Q := 1 2 3;
param sbase:=
1 800
2 1000 
3 1200;

What i have done so far is that:(really don't know how to do this..)

from pyomo.environ import *
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import random

# create a model
model = AbstractModel()

# declare decision variables
model.Q = Set()
model.J = Set()
model.A = Param(model.J, within = PositiveReals)
model.B = Param(model.J, within = PositiveReals)
model.C = Param(model.J, within = PositiveReals)
model.P_min = Param(model.J, within = PositiveReals)
model.P_max= Param(model.J, within = PositiveReals)
model.P = Var(model.J, within = NonNegativeReals)
model.sbase=Param(model.Q, within = PositiveReals)

# declare objective
def obj_expression(model):
    return sum(model.A[j] * model.P[j]**2 + model.B[j] * model.P[j] + model.C[j] for j in model.J)
model.profit = Objective(rule = obj_expression, sense=minimize)

# declare constraints
def lower_bound(model,j):
    return model.P_min[j] <= model.P[j]
model.laborA = Constraint(model.J, rule = lower_bound)

def upper_bound(model, j):
    return model.P[j] <= model.P_max[j]
model.laborB = Constraint(model.J, rule = upper_bound)

for q in range(2):    #from this line ,really confused about this
    def sum_labor(model,i,q):???
        if q==1:??
            return sum(model.P[j] for j in model.J) >= (model.sbase[Q] for q in model.Q)??
        elif q==2:??
            reture sum(model.P[j] for j in model.J) > = (model.sbase[Q] for q in model.Q)
        
model.laborC = Constraint(model.sbase,rule = sum_labor)

instance = model.create_instance("E:\pycharm_project\PD\pd.dat")
opt = SolverFactory('Ipopt')
results = opt.solve(instance)



plt.figure(figsize=(8,6))
plt.subplot(121)
for i in instance.J:
    plt.bar(i,value(instance.P[i]))
plt.xlabel('Generators')
plt.ylabel('Power')
plt.title('Power distribution') 

plt.subplot(122)
x=[0.1]
plt.bar(x,value(instance.profit),alpha=0.7,width=0.015)
plt.xlim(0,0.2)
plt.xlabel('')
plt.ylabel('Cost')
plt.title('minimum cost') 
plt.show()
print('\n\n---------------------------')
print('Cost: ',value(instance.profit))
#DATA 
set Q := 1 2 3;
set J := g1 g2 g3 g4 g5 g6 g7 g8 g9 g10  ;

param : A B C P_min P_max:=
g1   0.0148 12.1  82  80  200
g2   0.0289 12.6  49  120 320
g3   0.0135 13.2  100 50  150
g4   0.0127 13.9  105 250 520
g5   0.0261 13.5  72  80   280
g6   0.0212 15.4  29  50   150
g7   0.0382 14.0  32  30   120
g8   0.0393 13.5  40  30   110
g9   0.0396 15.0  25  20   80
g10  0.0510 14.3  15  20   60
;
param sbase:=
1 800
2 1000 
3 1200;

Do you have any suggestion on how to do it?

Thanks!


Solution

  • Here is a toy example that I think answers your question. You can use something like this or still read your sbase variable in from data and then just index it to set the value in a similar way within the loop.

    from pyomo.environ import *
    
    m = ConcreteModel()
    
    m.I = Set(initialize=[1,])      # mixed type set for demo only.
    
    m.sbase = Param(mutable=True)   # for the purposes of each instantiation, sbase is FIXED, so no index reqd.
    
    m.x = Var(m.I, domain=NonNegativeReals )
    
    ## constraint
    m.C1 = Constraint(expr=sum(m.sbase*m.x[i] for i in m.I) <= 10 )
    
    ## objective
    m.OBJ = Objective(expr=sum(m.sbase*m.x[i] for i in m.I), sense=maximize)
    
    solver = SolverFactory('cbc')
    
    sbase_values = {1, 2, 5}
    for sbase in sbase_values:
        m.sbase = sbase
        results = solver.solve(m)
        print(f'\n ***** solved for sbase = {sbase} *****')
        m.display()