2. Modelar componentes en Pyomo.

Como vimos en un apartado anterior, existen dos tipos básicos de modelos para definir en Pyomo, El modelo concreto y el abstracto. Con el primero se define el modelo que se quiere optimizar para unos valores de los parámetros concretos, mientras que con el segundo formato se declara de forma genérica el modelo y después se utiliza el mismo para los parámetros concretos sobre los cuales se quiere hacer la optimización. En este apartado se darán más detalles de cómo construir estos dos tipos de modelos.

Pyomo soporta un diseño de programación orientada a objetos para la construcción de los modelos a optimizar. En este sentido son cuatro los pasos básicos a seguir para construir este modelo:

1.- Crear el modelo y declarar las componentes.

2.- Instanciar al modelo.

3.- Aplicar el solver

4.- Estudiar los resultados obtenidos.

En la práctica estos son los cuatro pasos básicos que se deben ejecutar, con diferentes datos o diversas restricciones para entender mejor la naturaleza del problema con el que estamos trabajando.

Un modelo Pyomo consiste en una serie de componentes del modelo ( modeling components) que definen diferentes aspectos del modelo que se quiere construir. Con Pyomo se pueden tener las diferentes modeling components que generalmente son soportadas por los modelos de tipo AMLs. En concreto las componentes del modelo que se pueden construir en Pyomo, están soportadas por las siguientes clases de tipo Python:

  • Set. Conjunto de datos utilizados para definir una instancia del modelo.

  • Param. Conjunto de parámetros que son utilizados para definir la instancia del modelo.

  • Var. Variables de decsión, que conformarán los resultados del modelo.

  • Objetive. Es la expresión matemática que hay que maximizar o minimizar.

  • Constraint. Son el conjunto de restricciones a tener en cuenta en la solución el problema.

Pymo ofrece dos tipos de contextos para definir e inicializar las componentes del modelo: ConcreteModel y AbstractModel. Y así po ejemplo, para modelar componentes de un modelo abstracto, se utilizará la siguiente expresión:

model = AbstractModel()
model.I = Set()
model.p = Param(model.I)

En este código, se obtiene un objeto denominado model que es una instancia de la clase AbstractModel, y model.I será un objeto de la clase Set, al igual que model.p es una instancia de clase Param. Pyomo también soporta los componentes indexados ( indexed components ), que son una colección de componentes que tienen como referencia un determinado índice, y que se verán con más detalle en apartados posteriores. En el ejemplo anterior, el parámetro p es un componente indexado que tiene como indice el elemento I.

Para entender la construcción de los modelos ( tanto concreto como abstractos ) que se pueden crear con Pyom, nos basaremos en los modelos matemáticos que se definen a continuación.

Para el modelo abstracto, se utilizará la siguiente expresión:

\[\begin{split} \begin{array}{lll} \min & \sum_{j=1}^n c_j x_j & \\ s.t. & \sum_{j=1}^n a_ij x_j \geq b_i & \forall i = 1 \ldots m\\ & x_j \geq 0 & \forall j = 1 \ldots n \end{array} \end{split}\]

Mientras que para un modelo concreto, se usarán los datos concretos del modelo general anterior que se muestran a continuación:

\[\begin{split} \begin{array}{lll} \min & 4x_1 + 3x_2\\ s.t. & 3x_1 + 4x_2 \geq 1\\ & 2x_1 + 6x_2 \geq 2 \\ & x_1,x_2 \geq 0 \end{array} \end{split}\]

2.1. Tipos de parámetros o variables.

En Pyomo se pueden utilizar los siguientes tipos de parámetros o variables, y que en apartados posteriores se verán cómo se pueden utilizar ( todos estos parámetros se pueden ver en este enlace).

  • Any. Para cualquier valor.

  • Reals. Valores de punto flotante.

  • PositiveReals. Valores de punto flotante estrictamente positivos.

  • NonPositiveReals. Valores de punto flotante no-positivos.

  • NegativeReals. Valores de punto flotante estrictamente negativos.

  • NonNegativeReals. Valores de punto flotante no-negativos.

  • PercentFraction. Valores de punto flotante dentro del intervalo [0,1].

  • UnitInterval. Es un alias para PercentFraction.

  • Integers. Valores enteros

  • PositiveIntegers. Enteros positivos.

  • NonPositiveIntegers. Enteros no-positivos.

  • NegativeIntegers. Valores enteros negativos.

  • NonNegativeIntegers. Valores enteros no negativos.

  • Boolean. Valores booleanos, pueden estar representados por False/True, 0/1, ‘False’/’True’ y ‘F’/’T’.

  • Binary. Los números enteros \(\{0,1\}\)

2.2. Modelos concretos en Pyomo.

Esta es la opción más rápida par construir un modelo con unos parámetros concretos, pero que servirá para resolver ese modelo en concreto, si se quiere generalizar el modelo y que sirva para un conjunto de parámetros diferentes, se tendrá que utilizar un modelo abstracto que se desarrollará posteriormente.

Para definir un modelo concreto, se tendrá que utilizar la clase denominada ConcreteModel(), de la siguiente manera:

model = ConcreteModel()

A diferencia de lo que ocurre con un modelo abstracto, en estos modelos de tipo concreto, los modelos de componentes dentro del modelo son inmediatamente inicializadas y son añadidas a la isntancia del modelo.

Las variables de decisión para este ejemplo concreto se declarán de la siguiente manera:

model.x_1 = Var(within=NonNegativeReals)
model.x_2 = Var(within=NonNegativeReals)

El argumento within se utiliza para definir el rango de valores que soporta la variable. Los valores que se pueden usar para este argumento, se pueden ver en este apartado anterior.

Otro apartado imprescindible a la hora de construir los modelos en Pyomo es la función objetivo. Para hacer esto, se debe construir una expresión algebraica compuesta de constantes numéricas y variables que siguen a la palabra clave expr. Para el ejemplo anterior, se utilizará la siguiente expresión:

model.obj = Objective(expr=4*model.x_1 + 3*model.x_2)

Igualmente se necesitan definir las restricciones del modelo lineal, que para el ejemplo que se sigue, serían de la siguiente manera:

model.con1 = Constraint(expr=3*model.x_1 + 4*model.x_2 >= 1)
model.con2 = Constraint(expr=2*model.x_1 + 6*model.x_2 >= 2)

Uniendo todas estas expresiones, el modelo concreto que debería definirse para este caso, sería el siguiente:

from pyomo.environ import *

model = ConcreteModel()
model.x_1 = Var(within=NonNegativeReals)
model.x_2 = Var(within=NonNegativeReals)
model.obj = Objective(expr=4*model.x_1 + 3*model.x_2)
model.con1 = Constraint(expr=3*model.x_1 + 4*model.x_2 >= 1)
model.con2 = Constraint(expr=2*model.x_1 + 6*model.x_2 >= 2)

2.2.1. Variables y parámetros indexados.

El modelo anterior queda muy bien construido de la forma que lo hemos hecho en el anterior apartado, pero tiene el inconveniente de que es muy difícil de construir y mantener, en el caso de que se disponga de muchas variables. Para salvar esta dificultada, Pymo tiene previsto la utilización de variables y parámetros de tipo indexado, y su utilización y uso se pasa a explicar a continuación.

Por ejemplo para el caso anterior, si se utilizan variables indexadas, se utilizará una construcción similar a la siguiente:

model.x = Var([1,2], within=NonNegativeReals)

Los índices de las variables, se podrían utilizar para construir tanto para construir la función objetivo como las restricciones, de la siguiente manera.

model.obj = Objective(expr=4*model.x[1] + 3*model.x[2])
model.con1 = Constraint(expr=3*model.x[1] + 4*model.x[2]>=1)
model.con2 = Constraint(expr=2*model.x[1] + 6*model.x[2]>=2)

El modelo conjunto con este tipo de planteamiento sería el siguiente:

from pyomo.environ import *

model = ConcreteModel()
model.x = Var([1,2], within=NonNegativeReals)
model.obj = Objective(expr=4*model.x[1] + 3*model.x[2])
model.con1 = Constraint(expr=3*model.x[1] + 4*model.x[2]>=1)
model.con2 = Constraint(expr=2*model.x[1] + 6*model.x[2]>=2)

2.2.2. Aportación de datos de forma externa.

Como puede verse con el anterior modelo, se mejora el indicado en primer lugar, sin embargo adolece de un impportante problema: que los datos están inscrustados en la definición del modelo. Para solucionar esto, existe la posibilidad de definir por un lado los datos de forma independiente, después usar estos para definir el modelo. Con esto se gana en claridad en el código y permite un mejor mantenimiento del mismo.

Para conseguir esto, se pueden definir los datos a utilizar de forma independiente de la la siguiente manera (observar el uso muy extendido en este tipo de problemas de los diccionarios):

N = [1,2]
c = {1:4, 2:3}
a = {(1,1):3, (2,1):4, (1,2):2, (2,2):6}
b = {1:1, 2:2}

Entonces utilizando estos datos, la función objetivo se puede definir de la siguiente manera:

model.obj = Objective(expr=c[1]*model.x[1] + c[2]*model.x[2])

expresión anterior que es perfectamente válida pero que adole de un serio inconveniente: cuando ay muchos variables y parámetros, se hace inmanejable a medida que dichos parámetros y variables se incrementan. Para resolver esto, se puede utilizar la función sum de la siguiente manera:

model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))

De forma similar a la anterior, se plantearían las restricciones, por lo que el modelo que nos quedará en esta caso sería el siguiente:

from pyomo.environ import *
# defino los valores
N = [1,2]
c = {1:4, 2:3}
a = {(1,1):3, (2,1):4, (1,2):2, (2,2):6}
b = {1:1, 2:2}
#defino el modelo
model = ConcreteModel()
model.x = Var(N, within=NonNegativeReals)
model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))
model.con1 = Constraint(expr=sum(a[i,1]*model.x[i] for i in N) >= b[1])
model.con2 = Constraint(expr=sum(a[i,2]*model.x[i] for i in N) >= b[2])

2.2.3. Modelo con definición restricciones de tipo Rules

Aún se puede mejorar la construcción de los modelos en Pymo, utilizando lo que en términos anglosajones de denominan rules. En concreto, utilizando estas herramientas, el modelo anterior se puede construir de la siguiente manera:

from pyomo.environ import *

N = [1,2]
M = [1,2]
c = {1:4, 2:3}
a = {(1,1):3, (2,1):4, (1,2):2, (2,2):6}
b = {1:1, 2:2}

model = ConcreteModel()
model.x = Var(N, within=NonNegativeReals)
model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))
def con_rule(model, m): #definimos la rule
    return sum(a[i,m]*model.x[i] for i in N) >= b[m]
#Ahora utilizamos esa rule
model.con = Constraint(M, rule=con_rule)

Se puede generalizar el uso de rules, para construir los modelos den Pyomo, como se puede ver en el siguiente ejemplo:

from pyomo.environ import *

model = ConcreteModel()

N = [1,2]

model.N = Set(initialize= N)
model.M = Set(initialize=[1,2])

model.c = Param(model.N, initialize={1:4, 2:3})
model.a = Param(model.N, model.M,initialize={(1,1):3, (2,1):4, (1,2):2, (2,2):6})
model.b = Param(model.M, initialize={1:1, 2:2})
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[i,m]*model.x[i] for i in model.N) \
    >= model.b[m]
model.con = Constraint(model.M, rule=con_rule)

2.2.4. Usando una función para construir el modelo.

Una baza importante a jugar con Python y poder usar con mayor agilidad las ventajas que ofrece Pyomo es poder utilizar las funciones de Python, para poder ser utilizadas de forma flexible en el código. A continuación se muestra un ejemplo de un posible uso, para el ejemplo sobre el que estamos trabajando.

from pyomo.environ import *

def create_model(N=[], M=[], c={}, a={}, b={}):
    model = ConcreteModel()
    model.x = Var(N, within=NonNegativeReals)
    model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))
    def con_rule(model, m):
        return sum(a[i,m]*model.x[i] for i in N) >= b[m]
    model.con = Constraint(M, rule=con_rule)
    return model

model = create_model(N = [1,2], M = [1,2], c = {1:4, 2:3},
    a = {(1,1):3, (2,1):4, (1,2):2, (2,2):6},
    b = {1:1, 2:2})

2.2.5. Resolución de modelos concretos.

Para resolver este tipo de modelos, existen dos posibilidades. Una es mediante la inclusión de código dentro del mismo entorno que donde está definido el modelo, como se muestra en el ejemplo que sigue.

from pyomo.environ import *

model = ConcreteModel()
model.x = Var([1,2], within=NonNegativeReals)
model.obj = Objective(expr=4*model.x[1] + 3*model.x[2])
model.con1 = Constraint(expr=3*model.x[1] + 4*model.x[2]>=1)
model.con2 = Constraint(expr=2*model.x[1] + 6*model.x[2]>=2)

model.pprint()

opt = SolverFactory('glpk')
opt.solve(model) 
1 Set Declarations
    x_index : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

1 Var Declarations
    x : Size=2, Index=x_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  None :  None : False :  True : NonNegativeReals
          2 :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 4*x[1] + 3*x[2]

2 Constraint Declarations
    con1 : Size=1, Index=None, Active=True
        Key  : Lower : Body            : Upper : Active
        None :   1.0 : 3*x[1] + 4*x[2] :  +Inf :   True
    con2 : Size=1, Index=None, Active=True
        Key  : Lower : Body            : Upper : Active
        None :   2.0 : 2*x[1] + 6*x[2] :  +Inf :   True

5 Declarations: x_index x obj con1 con2
{'Problem': [{'Name': 'unknown', 'Lower bound': 1.0, 'Upper bound': 1.0, 'Number of objectives': 1, 'Number of constraints': 3, 'Number of variables': 3, 'Number of nonzeros': 5, 'Sense': 'minimize'}], '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.14016246795654297}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

Otra forma e mediante la línea de comandos, utilizando el comando pyomo solve my_model.py –solver=’glpk’. A continuación, se muestra cómo hacer eso con nuestro ejemplo. Tener en cuenta que al estar trabajando sobre jupyter notebook para escribir todo esto, se utiliza el comando %%writefile codigos/model2.py para guardar el código, como si fuera un fichero python, y posteriormente se ejecuta su código con la instrucción !pyomo solve codigos/model2.py –solver=glpk.

%%writefile codigos/model2.py
from pyomo.environ import *

model = ConcreteModel()
model.x = Var([1,2], within=NonNegativeReals)
model.obj = Objective(expr=4*model.x[1] + 3*model.x[2])
model.con1 = Constraint(expr=3*model.x[1] + 4*model.x[2]>=1)
model.con2 = Constraint(expr=2*model.x[1] + 6*model.x[2]>=2)
Overwriting codigos/model2.py
!pyomo solve codigos/model2.py --solver=glpk
[    0.00] Setting up Pyomo environment
[    0.00] Applying Pyomo preprocessing actions
[    0.02] Creating model
[    0.02] Applying solver
[    0.06] Processing results
    Number of solutions: 1
    Solution Information
      Gap: 0.0
      Status: feasible
      Function Value: 0.9999999999999989
    Solver results file: results.json
[    0.06] Applying Pyomo postprocessing actions
[    0.06] Pyomo Finished

2.3. Modelos abstractos en Pyomo.

Una vez vistas las diferentes modalidades de modelos concretos ya presentadas en los apartado anteriores, el modelo de tipo abstracto para trabajar con Pyomo es fácil de definir, en concreto el ejemplo ya mostrado anteriormente se definiría de la siguiente manera:

from pyomo.environ import *

model = AbstractModel() 

model.N = Set()
model.M = Set() 

model.c = Param(model.N)
model.a = Param(model.N, model.M)
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[i,m]*model.x[i] for i in model.N) \
        >= model.b[m]
model.con = Constraint(model.M, rule=con_rule)

Ahora el modelo abstracto necesita alimentarse de una serie de datos, los cuales deben ser suministrados mediante un fichero que tiene el siguiente formato (llamaremos al fichero datos.dat):

set N := 1 2 ;
set M := 1 2 ;
param c :=
1 4
2 3 ;
param a :=
1 1 3
2 1 4
1 2 2
2 2 6 ;
param b :=
1 1
2 2 ;
```

Ahora para resolverlo desde las instrucciones del propio programa habría que crear una instancia, a la que se le indica además el fichero de datos con el que debe trabajar, y después la ejecución sería similar al la del modelo concreto ya expuesto en apartados posteriores. En concreto, las instrucciones serían las siguientes:

model.pprint()

instance = model.create_instance('dat/datos.dat')


print("______________________________________________________________")
opt = SolverFactory('glpk')
opt.solve(instance) 
3 Set Declarations
    M : Size=0, Index=None, Ordered=Insertion
        Not constructed
    N : Size=0, Index=None, Ordered=Insertion
        Not constructed
    a_index : Size=0, Index=None, Ordered=True
        Not constructed

3 Param Declarations
    a : Size=0, Index=a_index, Domain=Any, Default=None, Mutable=False
        Not constructed
    b : Size=0, Index=M, Domain=Any, Default=None, Mutable=False
        Not constructed
    c : Size=0, Index=N, Domain=Any, Default=None, Mutable=False
        Not constructed

1 Var Declarations
    x : Size=0, Index=N
        Not constructed

1 Objective Declarations
    obj : Size=0, Index=None, Active=True
        Not constructed

1 Constraint Declarations
    con : Size=0, Index=M, Active=True
        Not constructed

9 Declarations: N M c a_index a b x obj con
______________________________________________________________
{'Problem': [{'Name': 'unknown', 'Lower bound': 1.0, 'Upper bound': 1.0, 'Number of objectives': 1, 'Number of constraints': 3, 'Number of variables': 3, 'Number of nonzeros': 5, 'Sense': 'minimize'}], '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.03125}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

Si la ejecución es desde la línea de comandos, primero creariamos un fichero con extensión .py

%%writefile codigos/model3.py
from pyomo.environ import *

model = AbstractModel() 

model.N = Set()
model.M = Set() 

model.c = Param(model.N)
model.a = Param(model.N, model.M)
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[i,m]*model.x[i] for i in model.N) \
        >= model.b[m]
model.con = Constraint(model.M, rule=con_rule)
Overwriting codigos/model3.py

Y después lo ejecutamos desde la linea de comandos con la siguiente instrucción:

!pyomo solve codigos/model3.py dat/datos.dat --solver=glpk
[    0.00] Setting up Pyomo environment
[    0.00] Applying Pyomo preprocessing actions
[    0.11] Creating model
[    0.13] Applying solver
[    0.19] Processing results
    Number of solutions: 1
    Solution Information
      Gap: 0.0
      Status: feasible
      Function Value: 0.9999999999999989
    Solver results file: results.json
[    0.19] Applying Pyomo postprocessing actions
[    0.19] Pyomo Finished