Categoria: Pyomo

Pyomo Abstract exemplo 1

Exemplo de modelo abstrato no Pyomo.

Link para download.

A grande vantagem é que separa claramente os dados do modelo.

Uma desvantagem é que o formato de dados tem que carregar os índices das matrizes, de forma até redundante.

Portanto, creio que vale a pena para modelos grandes, onde a chance de infeasible é alta e uma organização melhor do problema se paga.

Dados do arquivo dat:

param : N : c :=
1 1
2 2;

param : M : b :=
1 1
2 2;

param a :=
1 1 3
1 2 4
2 1 2
2 2 5;

Modelo abstrato:

–– coding: utf-8 ––

“””
Created on Mon Jul 13 05:48:48 2020

@author: asgun
“””

from pyomo.environ import *

model = AbstractModel()

model.N = Set()
model.M = Set()
model.c = Param(model.N)
model.a = Param(model.M, model.N)
model.b = Param(model.M)

model.x = Var(model.N, within= NonNegativeReals)

def obj_rule(model):
return sum(model.c[i]*model.x[i] for i in model.N)
model.obj = Objective(rule = obj_rule)

def con_rule(model, m):
return sum(model.a[m,i]*model.x[i] for i in model.N) >= model.b[m]

model.con = Constraint (model.M, rule = con_rule)

opt = SolverFactory(‘glpk’)

instance = model.create_instance(“Abstract_data.dat”)
results = opt.solve(instance)

instance.display()

Pyomo concrete 02

Este é um exemplo do Pyomo, para um modelo concreto.

Porém, uma vantagem bastante grande é que os dados podem vir em listas e dicionários do Python, e as regras da função objetivo e restrições podem ser escritas como funções do Python.

Como os projetos normalmente serão de tamanho médio, eu acho este o método mais interessante, por manipular via Python o arquivo inteiro.

Projetos muito mais pesados teriam que ser resolvidos com solvers mais poderosos, como o CPLEX, e aí a coisa muda de nível – é melhor usar AIMMS ou algo mais robusto, até para fins de debug.

Link para download.

Código:

from pyomo.environ import *

#dados

N = [1,2]
M = [1,2]

c = {1:1, 2:2}
a = {(1,1):3, (1,2):4, (2,1):2, (2,2):5}
b = {1:1, 2:2}

#Model

model = ConcreteModel()

model.x = Var(N, within = NonNegativeReals)

def obj_rule(model):
return sum(c[i] + model.x[i] for i in N)

model.obj = Objective(rule = obj_rule)

def cons_1(model,m):
return sum(a[m,i]*model.x[i] for i in N) >= b[m]

model.con = Constraint(M, rule = cons_1)

opt = SolverFactory(‘cbc’) #Ou glpk
results = opt.solve(model)

print(results)
model.display()

Pyomo Exemplo 1

Link para download.

Exemplo mais simples possível do Pyomo.

É um modelo concreto, que resolve

Minimizar x1 + 2*x2

s.a.

2*x1+4*x2>=1

2*x1+5*x2>=2

Código:

from pyomo.environ import *

model = ConcreteModel()
model.x1 = Var(within=NonNegativeReals)
model.x2 = Var(within=NonNegativeReals)
model.obj= Objective(expr = model.x1 + 2*model.x2, sense=minimize)

model.con1= Constraint(expr = 2model.x1 +4model.x2 >=1)
model.con2 = Constraint(expr = 2model.x1 + 5model.x2>=2)

opt = SolverFactory(‘cbc’) #Ou glpk
results = opt.solve(model)

print(results)
model.display()

Resolvendo o Set Cover com OpenSolver e Pyomo

A ideia deste exercício é resolver o problema da cobertura de conjuntos, utilizando o Pyomo (http://www.pyomo.org/) e o OpenSolver (www.opensolver.org).

Link para download: https://github.com/asgunzi/SetCover.

O Set cover é o problema de cobertura de conjuntos.

Tenho uma quantidade de elementos, digamos 50.

E tenho uma quantidade de conjuntos, digamos 15.

Cada conjunto tem um preço (abaixo está como peso), e cada conjunto contém os elementos.

Algumas aplicações: imagine que há várias doenças, e uma série de vacinas. Cada vacina cobre um conjunto diferente de doenças, e cada vacina tem um custo. Minimizar o custo total, cobrindo todas as doenças.

Open Solver:

Variável de decisão binária: utilizo ou não o conjunto.

Função objetivo: se utilizo o conjunto, é o preço deste.

Restrição: devo cobrir todos os elementos.

Resultado: FO de 223 e as escolhas pintadas em rosa.

Formulação Pyomo:

Leitura de dados utilizando o módulo pandas.

    df1 = pd.read_excel(io=’SetCover_Instancia01.xlsm’, sheet_name=’Base’, usecols = “B:P”, skiprows = 1, nrows = 1,header = None)

    df2 = pd.read_excel(io=’SetCover_Instancia01.xlsm’, sheet_name=’Base’, usecols = “B:P”, skiprows = 7, header = None)

    Nsets =  len(df1.columns)

    Nelementos =  len(df2)

Crio o model com o primeiro comando abaixo. Defino os índices (conjuntos e elementos).

   #—————- Modelo Pyomo ————#     

    model = ConcreteModel()

    model.S = list(range(0,Nsets)) #Num sets

    model.E = list(range(0,Nelementos)) #Num elementos

Defino a variável de decisão com o comando Var. O nome da variável é AtivaSet, ela é binária e inicializa com o valor zero.

    model.AtivaSet = Var(model.S, within=Binary, initialize=0)

Função objetivo.

Soma do produto entre a variável de decisão e o peso de cada conjunto.

def obj_rule(model):

        return sum(model.AtivaSet[s]*pesos[s] for s in model.S)

model.obj = Objective(rule = obj_rule, sense=minimize)

Restrição.

Para cada elemento, ele deve pertencer a um conjunto, ou seja, a soma do produto da matriz de alocação e a variável de decisão devem ser maior do que 1.

    def constrAlocacao(model, e):

        return sum(model.AtivaSet[s]*Malocacao[e][s] for s in model.S) >= 1

    model.restr1 = Constraint(model.E, rule = constrAlocacao)      

Algumas configurações do solver, e depois mando resolver

    SOLVER_NAME = ‘cbc’

    TIME_LIMIT = 60*1   #Em segundos

    solver = SolverFactory(SOLVER_NAME)

    solver.options[‘seconds’] = TIME_LIMIT

    results  = solver.solve(model)

    results.write() #para mostrar o status da solucao

O resultado obtido foi o mesmo do OpenSolver.

Além do CBC, outro solver free é o GLPK.

Todos os outros são pagos. O Pyomo pode chamar CPLEX, Gurobi, se tiver a licença destes.

Toda a manipulação de entrada de dados, matrizes, e saída de dados, é Python puro.

Vantagens: é free, e muito melhor que o Open Solver para manipular grande quantidade de dados.

Desvantagens: requer um grau maior de conhecimento que o OpenSolver. Não é tão fácil de manipular ou tão poderoso quanto um AIMMS + CPLEX.