==== Planejamento Regional de Kibutzim - Exemplo de Programação Linear usando o Pyomo1 ====

Autor: Heitor Gabriel S. Monteiro

1) Apresentação da Questão e Formulação do Problema

1.1) A questão apresentada segue Hillier & Lieberman, Introdução à Pesquisa Operacional, 9a ed., capítulo 3, página 41.

A CONFEDERAÇÃO MERIDIONAL DE KIBUTZIM é um grupo de três kibutzim (comunidades agrícolas coletivas em Israel). O planejamento geral para esses grupos é feito em seu Centro Técnico de Coordenação. Esse escritório está planejando atualmente a produção agrícola para o próximo ano. A produção agrícola de cada kibutz é limitada tanto pela quantidade de área irrigável disponível como pela quantidade de água alocada para a irrigação pelo Comissariado de Recursos Hídricos (um órgão governamental). Esses dados são fornecidos na Tabela 3.8. Entre as plantações adequadas para essa região encontram-se beterraba, algodão e sorgo e são estas que estão sendo consideradas para o próximo período. Essas plantações diferem basicamente nos respectivos retornos líquidos esperados e consumo de água. Além disso, o Ministério de Agricultura tem uma cota máxima para a área total que pode ser dedicada a cada uma dessas plantações pela Confederação Meridional de Kibutzim, conforme ilustrado na Tabela 3.9. Em razão da limitada disponibilidade de água para irrigação, a Confederação Meridional de Kibutzim não será capaz de usar toda sua área irrigável para plantação de culturas na próxima temporada. Para garantir equilíbrio entre os três kibutzim, foi acordado que cada um deles vai plantar a mesma proporção de sua área irrigável. Por exemplo, se o kibutz 1 plantar 200 de seus 400 acres disponíveis, então o kibutz 2 terá de plantar 300 de seus 600 acres, ao passo que o kibutz 3 plantaria 150 de seus 300 acres. Entretanto, qualquer combinação das plantações pode ser cultivada no kibutzim. A tarefa que o Centro Técnico de Coordenação deve enfrentar é planejar quantos acres devem ser dedicados a cada plantação no respectivo kibutzim satisfazendo as dadas restrições. O objetivo é maximizar o retomo líquido total para a Confederação Meridional do Kibutzim como um todo.

Tabela 3.8: Dados de Recursos para a Confederação Meridional de Kibutzim

Kibutzim	Terra Utilizável (acres)	Locação de Água (acres pés)
1	400	600
2	600	800
3	300	375

Tabela 3.9: Dados de Plantações para a Confederação Meridional de Kibutzim

Plantação	Cota Máxima (acres)	Consumo de Água (acres pés/acres)	Retorno Líquido (US$/acre)
Beterraba	600	3	1000
Algodão	500	2	750
Sorgo	325	1	250

1.2) Definição Matemática do Problema

1.2.1) Variáveis de Decisão

• $b_{1}$: área plantada com beterraba no kibutz 1 (acres)
• $a_{1}$: área plantada com algodão no kibutz 1 (acres)
• $s_{1}$: área plantada com sorgo no kibutz 1 (acres)
• $b_{2}$: área plantada com beterraba no kibutz 2 (acres)
• $a_{2}$: área plantada com algodão no kibutz 2 (acres)
• $s_{2}$: área plantada com sorgo no kibutz 2 (acres)
• $b_{3}$: área plantada com beterraba no kibutz 3 (acres)
• $a_{3}$: área plantada com algodão no kibutz 3 (acres)
• $s_{3}$: área plantada com sorgo no kibutz 3 (acres)

1.2.2) Função Objetivo

• $\max Z = 1000 \sum^3_1 b_{i} + 750 \sum^3_1 a_{i} + 250 \sum^3_1 s_{i}$

1.2.3) Restrições

Terra utilizável:

• $b_{1} + a_{1} + s_{1} \leq 400$
• $b_{2} + a_{2} + s_{2} \leq 600$
• $b_{3} + a_{3} + s_{3} \leq 300$

Alocação de água:

• $3b_{1} + 2a_{1} + s_{1} \leq 600$
• $3b_{2} + 2a_{2} + s_{2} \leq 800$
• $3b_{3} + 2a_{3} + s_{3} \leq 375$

Área total para cada plantação:

• $b_{1} + b_{2} + b_{3} \leq 600$
• $a_{1} + a_{2} + a_{3} \leq 500$
• $s_{1} + s_{2} + s_{3} \leq 325$

Proporcionalidade entre os kibutzim:

$ \frac{b_{1} + a_{1} + s_{1}}{400} = \frac{b_{2} + a_{2} + s_{2}}{600} = \frac{a_{3} + s_{3} + b_{3}}{300} $

As igualdades abaixo tornam-se as três restrições a seguir:

• $3(b_{1} + a_{1} + s_{1}) - 2(b_{2} + a_{2} + s_{2}) = 0$
• $4(b_{1} + a_{1} + s_{1}) - 3(b_{3} + a_{3} + s_{3}) = 0$
• $2(b_{2} + a_{2} + s_{2}) - (b_{3} + a_{3} + s_{3}) = 0$

Não-negatividade:

• $b_{1}, a_{1}, s_{1}, b_{2}, a_{2}, s_{2}, b_{3}, a_{3}, s_{3} \geq 0$

2) Solução do Problema no Pyomo

import numpy as np
import pandas as pd
import pyomo.environ as pyo

2.1) Dados

# Tabela 3.8 em pandas dataframe
df1 = pd.DataFrame({
    'kibutz': ['k1', 'k2', 'k3'],
    'area': [400, 600, 300],
    'agua': [600, 800, 375]},
                   index=['k1', 'k2', 'k3'])

df2 = pd.DataFrame({
    'plantacao': ['beterraba', 'algodao', 'sorgo'],
    'cota':[600,500,325],
    'consumo':[3,2,1],
    'lucro':[1000,750,250]},
                   index=['beterraba', 'algodao', 'sorgo'])

df1

	kibutz	area	agua
k1	k1	400	600
k2	k2	600	800
k3	k3	300	375

df2

	plantacao	cota	consumo	lucro
beterraba	beterraba	600	3	1000
algodao	algodao	500	2	750
sorgo	sorgo	325	1	250

2.2) Declaração do modelo:

# Nome
kibutz = pyo.ConcreteModel()

# Indices
kibutz.I = pyo.Set(initialize = df1.kibutz.tolist())    # unidades produtivas
kibutz.J = pyo.Set(initialize = df2.plantacao.tolist()) # culturas

# Parametros
kibutz.area = pyo.Param(kibutz.I, initialize=df1.area.to_dict())
kibutz.agua = pyo.Param(kibutz.I, initialize=df1.agua.to_dict())
kibutz.cota = pyo.Param(kibutz.J, initialize=df2.cota.to_dict())
kibutz.consumo = pyo.Param(kibutz.J, initialize=df2.consumo.to_dict())
kibutz.lucro = pyo.Param(kibutz.J, initialize=df2.lucro.to_dict())

# Variaveis
kibutz.x = pyo.Var(kibutz.I, kibutz.J, domain=pyo.NonNegativeReals)

# Funcao objetivo
def obj_expression(kibutz):
    return sum(kibutz.lucro[j] * sum(kibutz.x[i,j] for i in kibutz.I) for j in kibutz.J)

kibutz.lucro_total = pyo.Objective(rule=obj_expression, sense=pyo.maximize)

# Restricoes
def area_rule(kibutz, i):
    return sum(kibutz.x[i,j] for j in kibutz.J) <= kibutz.area[i]

kibutz.area_rule = pyo.Constraint(kibutz.I, rule=area_rule)

def agua_rule(kibutz, i):
    return sum(kibutz.x[i,j]*kibutz.consumo[j] for j in kibutz.J) <= kibutz.agua[i]

kibutz.agua_rule = pyo.Constraint(kibutz.I, rule=agua_rule)

def cota_rule(kibutz, j):
    return sum(kibutz.x[i,j] for i in kibutz.I) <= kibutz.cota[j]

kibutz.cota_rule = pyo.Constraint(kibutz.J, rule=cota_rule)

kibutz.pprint()

3 Set Declarations
    I : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {'k1', 'k2', 'k3'}
    J : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {'beterraba', 'algodao', 'sorgo'}
    x_index : Size=1, Index=None, Ordered=True
        Key  : Dimen : Domain : Size : Members
        None :     2 :    I*J :    9 : {('k1', 'beterraba'), ('k1', 'algodao'), ('k1', 'sorgo'), ('k2', 'beterraba'), ('k2', 'algodao'), ('k2', 'sorgo'), ('k3', 'beterraba'), ('k3', 'algodao'), ('k3', 'sorgo')}

5 Param Declarations
    agua : Size=3, Index=I, Domain=Any, Default=None, Mutable=False
        Key : Value
         k1 :   600
         k2 :   800
         k3 :   375
    area : Size=3, Index=I, Domain=Any, Default=None, Mutable=False
        Key : Value
         k1 :   400
         k2 :   600
         k3 :   300
    consumo : Size=3, Index=J, Domain=Any, Default=None, Mutable=False
        Key       : Value
          algodao :     2
        beterraba :     3
            sorgo :     1
    cota : Size=3, Index=J, Domain=Any, Default=None, Mutable=False
        Key       : Value
          algodao :   500
        beterraba :   600
            sorgo :   325
    lucro : Size=3, Index=J, Domain=Any, Default=None, Mutable=False
        Key       : Value
          algodao :   750
        beterraba :  1000
            sorgo :   250

1 Var Declarations
    x : Size=9, Index=x_index
        Key                 : Lower : Value : Upper : Fixed : Stale : Domain
          ('k1', 'algodao') :     0 :  None :  None : False :  True : NonNegativeReals
        ('k1', 'beterraba') :     0 :  None :  None : False :  True : NonNegativeReals
            ('k1', 'sorgo') :     0 :  None :  None : False :  True : NonNegativeReals
          ('k2', 'algodao') :     0 :  None :  None : False :  True : NonNegativeReals
        ('k2', 'beterraba') :     0 :  None :  None : False :  True : NonNegativeReals
            ('k2', 'sorgo') :     0 :  None :  None : False :  True : NonNegativeReals
          ('k3', 'algodao') :     0 :  None :  None : False :  True : NonNegativeReals
        ('k3', 'beterraba') :     0 :  None :  None : False :  True : NonNegativeReals
            ('k3', 'sorgo') :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    lucro_total : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 1000*(x[k1,beterraba] + x[k2,beterraba] + x[k3,beterraba]) + 750*(x[k1,algodao] + x[k2,algodao] + x[k3,algodao]) + 250*(x[k1,sorgo] + x[k2,sorgo] + x[k3,sorgo])

3 Constraint Declarations
    agua_rule : Size=3, Index=I, Active=True
        Key : Lower : Body                                              : Upper : Active
         k1 :  -Inf : 3*x[k1,beterraba] + 2*x[k1,algodao] + x[k1,sorgo] : 600.0 :   True
         k2 :  -Inf : 3*x[k2,beterraba] + 2*x[k2,algodao] + x[k2,sorgo] : 800.0 :   True
         k3 :  -Inf : 3*x[k3,beterraba] + 2*x[k3,algodao] + x[k3,sorgo] : 375.0 :   True
    area_rule : Size=3, Index=I, Active=True
        Key : Lower : Body                                          : Upper : Active
         k1 :  -Inf : x[k1,beterraba] + x[k1,algodao] + x[k1,sorgo] : 400.0 :   True
         k2 :  -Inf : x[k2,beterraba] + x[k2,algodao] + x[k2,sorgo] : 600.0 :   True
         k3 :  -Inf : x[k3,beterraba] + x[k3,algodao] + x[k3,sorgo] : 300.0 :   True
    cota_rule : Size=3, Index=J, Active=True
        Key       : Lower : Body                                                : Upper : Active
          algodao :  -Inf :       x[k1,algodao] + x[k2,algodao] + x[k3,algodao] : 500.0 :   True
        beterraba :  -Inf : x[k1,beterraba] + x[k2,beterraba] + x[k3,beterraba] : 600.0 :   True
            sorgo :  -Inf :             x[k1,sorgo] + x[k2,sorgo] + x[k3,sorgo] : 325.0 :   True

13 Declarations: I J area agua cota consumo lucro x_index x lucro_total area_rule agua_rule cota_rule

resultado = pyo.SolverFactory('glpk').solve(kibutz)
resultado.write()

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 633333.333333333
  Upper bound: 633333.333333333
  Number of objectives: 1
  Number of constraints: 10
  Number of variables: 10
  Number of nonzeros: 28
  Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.0028123855590820312
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

kibutz.x.pprint()

x : Size=9, Index=x_index
    Key                 : Lower : Value            : Upper : Fixed : Stale : Domain
      ('k1', 'algodao') :     0 :            300.0 :  None : False : False : NonNegativeReals
    ('k1', 'beterraba') :     0 :              0.0 :  None : False : False : NonNegativeReals
        ('k1', 'sorgo') :     0 :              0.0 :  None : False : False : NonNegativeReals
      ('k2', 'algodao') :     0 :            200.0 :  None : False : False : NonNegativeReals
    ('k2', 'beterraba') :     0 : 133.333333333333 :  None : False : False : NonNegativeReals
        ('k2', 'sorgo') :     0 :              0.0 :  None : False : False : NonNegativeReals
      ('k3', 'algodao') :     0 :              0.0 :  None : False : False : NonNegativeReals
    ('k3', 'beterraba') :     0 :            125.0 :  None : False : False : NonNegativeReals
        ('k3', 'sorgo') :     0 :              0.0 :  None : False : False : NonNegativeReals

pd.DataFrame({i: [kibutz.x[i,j].value for j in kibutz.J] for i in kibutz.I}, index=kibutz.J)

	k1	k2	k3
beterraba	0.0	133.333333	125.0
algodao	300.0	200.000000	0.0
sorgo	0.0	0.000000	0.0