Extracción minera - Solución#

Tags: Parámetro binario, Mezcla, Recurso

Hide code cell source 
import os
# Por precaución, cambiamos el directorio activo de Python a aquel que contenga este notebook
if "optimizacion" in os.listdir():
    os.chdir(r"optimizacion/Formulaciones/7. Extraccion minera/")

Enunciado - Punto 1#

La minería de cielo abierto es una actividad industrial que consiste en remover grandes cantidades de suelo para extraer el mineral deseado. Este tipo de minas son comunes en la extracción de materiales como: arena, arcilla, cobre y carbón. Las directivas de una empresa de extracción minera desean explotar un conjunto de zonas (\(M\)), de las cuales se puede extraer un conjunto de diferentes tipos de carbón (\(K:\{1\): Antracita, \(2\): Hulla, \(3\): Turba, \(4\): Lignito\(\})\). Se sabe que de la zona \(i\in M\) sólo puede extraerse el carbón tipo \(j\in K\). Para ello, suponga que \(a_{ij}\) es un parámetro binario que toma el valor de 1 si en la zona \(i\in M\) se puede extraer carbón del tipo \(j\in K\) y 0 en el caso contrario. La Figura 1 presenta un ejemplo en el cual hay 36 zonas y 4 tipos de carbón. Por ejemplo, para la zona uno se tiene que \(a_{11}=0\), \(a_{12}=1\), \(a_{13}=0\) y \(a_{14}=0\) ya que de la zona 1 sólo puede extraerse el tipo de carbón 2 (hulla).

\[\begin{split} \bbox[5px, border:2px solid black]{ \begin{array} {cccccc} \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{\ 1\ } & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{\ 2\ } & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{\ 3\ } & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{\ 4\ } & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{\ 5\ } & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{\ 6\ } \\ \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{\ 7\ } & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{\ 8\ } & \bbox[Plum, 5px, border:2px solid Plum]{\ 9\ } & \bbox[Plum, 5px, border:2px solid Plum]{10} & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{11} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{12} \\ \bbox[Plum, 5px, border:2px solid Plum]{13} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{14} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{15} & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{16} & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{17} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{18} \\ \bbox[Plum, 5px, border:2px solid Plum]{19} & \bbox[LightGoldenRodYellow, 5px, border:2px solid LightGoldenRodYellow]{20} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{21} & \bbox[Plum, 5px, border:2px solid Plum]{22} & \bbox[Plum, 5px, border:2px solid Plum]{23} & \bbox[Plum, 5px, border:2px solid Plum]{24} \\ \bbox[Plum, 5px, border:2px solid Plum]{25} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{26} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{27} & \bbox[LightGoldenRodYellow, 5px, border:2px solid LightGoldenRodYellow]{28} & \bbox[LightGoldenRodYellow, 5px, border:2px solid LightGoldenRodYellow]{29} & \bbox[LightGoldenRodYellow, 5px, border:2px solid LightGoldenRodYellow]{30} \\ \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{31} & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{32} & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{33} & \bbox[Plum, 5px, border:2px solid Plum]{34} & \bbox[LightGoldenRodYellow, 5px, border:2px solid LightGoldenRodYellow]{35} & \bbox[Plum, 5px, border:2px solid Plum]{36} \\ \end{array} } \end{split}\]
\[ \bbox[5px, border:2px solid black]{ \begin{array} {cccc} \bbox[Plum, 5px, border:2px solid Plum]{\ \text{1: Antracita}\ } & \bbox[SkyBlue, 5px, border:2px solid SkyBlue]{\ \text{2: Hulla}\ } & \bbox[PaleGreen, 5px, border:2px solid PaleGreen]{\ \text{3: Turba}\ } & \bbox[LightGoldenRodYellow, 5px, border:2px solid LightGoldenRodYellow]{\ \text{4: Lignito}\ } \end{array} } \]

Cada tonelada de carbón extraída de la zona \(i\in M\) le cuesta a la empresa \(c_i\) pesos. Adicionalmente, cada zona tiene una capacidad máxima de extracción de carbón de \(n_i\) toneladas. La Tabla 1 muestra la información de capacidades máximas y costos para cada zona.

Tabla 1. Información de extracción de cada zona

Zona

Capacidad máxima de extracción de carbón (ton)

Costo por tonelada de carbón extraído ($)

1

189

16

2

196

6

3

143

11

4

136

8

5

106

5

6

151

25

7

170

16

8

129

17

9

184

25

10

122

8

11

146

15

12

190

8

13

160

10

14

109

20

15

133

6

16

198

17

17

138

6

18

107

20

19

117

5

20

150

8

21

171

25

22

103

11

23

157

8

24

143

7

25

170

7

26

130

28

27

140

26

28

126

27

29

180

9

30

153

24

31

108

15

32

132

14

33

105

22

34

145

20

35

145

19

36

114

8

Por último. Supon que las directivas desean explotar un mínimo de \(b_j\) toneladas de cada tipo de carbón \(j\in K\). La Tabla 2 presenta dichas cantidades.

Tabla 2. Requerimiento mínimo de cada tipo de carbón

Tipo de carbón

Mínimo de toneladas a explotar

Antracita

862

Huila

898

Turba

562

Lignito

742

Formula un programa lineal que le permita a la empresa decidir cuánto deben extraer en cada zona de manera que se cumplan los requerimientos a un mínimo costo.

Formulación#

a. Formula matemáticamente un modelo de optimización de forma general que represente la situación anterior. Defina clara y rigurosamente:

  • Conjuntos

  • Parámetros

  • Variables de decisión

  • Restricciones

  • Naturaleza de las variables

  • Función objetivo

Conjuntos#

  • \(M\): Zonas

  • \(K\): Tipos de carbón

Parámetros#

  • \(a_{ij}:\begin{cases}1\text{,}&\text{ si en la zona }i\in M\text{ se puede extraer carbón del tipo }j\in K \\ 0\text{,}& \text{ d.l.c} \end{cases}\)

  • \(c_i\): costo por cada tonelada de carbón extraída en la zona \(i\in M\)

  • \(n_i\): capacidad máxima de extracción de carbón (en toneladas) en la zona \(i\in M\)

  • \(b_j\): mínimo de toneladas a explotar del tipo de carbón \(j\in K\)

Variables de decisión#

  • \(x_i\): toneladas de carbón extraído en la zona \(i\in M\)

Restricciones#

Se garantiza que la cantidad de carbón extraído del tipo \(j\in K\) sea como mínimo \(b_j\):

\[ \sum_{i\in M\vert a_{ij} = 1}x_i \ge b_j,\ \forall j\in K; \]

Se garantiza que no se puedan extraer más de \(n_i\) toneladas de carbón en la zona \(i\in M\):

\[ x_i \le n_i,\ \forall i\in M; \]

Naturaleza de las Variables#

Solo pueden extraerse cantidades positivas de carbón de cada zona \(i\in M\):

\[ x_i \ge 0,\ \forall i\in M; \]

Función objetivo#

Debe minimizarse el costo total de extracción:

\[ \operatorname{mín}\ \sum_{i\in M}c_{i}x_{i} \]

Implementación#

b. Resuelve el modelo planteado utilizando la librería de PuLP en Python. ¿Cuál es la solución óptima del problema?

Librerías#

Importa la librería pulp para crear y resolver el modelo.

import pulp as lp

Conjuntos#

Define los conjuntos M y K que representan las zonas y los tipos de carbón respectivamente.

Recuerda que por conveniencia de preservar el orden de los elementos de los conjuntos, no siempre deberás definirlos con el tipo set.

# Zonas
M = list(range(1, 37))

# Tipos de carbón
K = ["Antracita", "Hulla", "Turba", "Lignito"]

Parámetros#

Define o importa los parámetros del modelo.

# Si en la zona i en M se puede extraer carbón tipo j en K
a = {
    (1, "Hulla"): 1,
    (2, "Hulla"): 1,
    (3, "Turba"): 1,
    (4, "Turba"): 1,
    (5, "Hulla"): 1,
    (6, "Turba"): 1,
    (7, "Hulla"): 1,
    (8, "Hulla"): 1,
    (9, "Antracita"): 1,
    (10, "Antracita"): 1,
    (11, "Hulla"): 1,
    (12, "Turba"): 1,
    (13, "Antracita"): 1,
    (14, "Turba"): 1,
    (15, "Turba"): 1,
    (16, "Hulla"): 1,
    (17, "Hulla"): 1,
    (18, "Turba"): 1,
    (19, "Antracita"): 1,
    (20, "Lignito"): 1,
    (21, "Turba"): 1,
    (22, "Antracita"): 1,
    (23, "Antracita"): 1,
    (24, "Antracita"): 1,
    (25, "Antracita"): 1,
    (26, "Turba"): 1,
    (27, "Turba"): 1,
    (28, "Lignito"): 1,
    (29, "Lignito"): 1,
    (30, "Lignito"): 1,
    (31, "Hulla"): 1,
    (32, "Turba"): 1,
    (33, "Hulla"): 1,
    (34, "Antracita"): 1,
    (35, "Lignito"): 1,
    (36, "Antracita"): 1,
}

for i in M:
    for j in K:
        if (i, j) not in a:
            a[i, j] = 0

# Costo por tonelada de carbón en la zona i en M
c = {
    1: 16,
    2: 6,
    3: 11,
    4: 8,
    5: 5,
    6: 25,
    7: 16,
    8: 17,
    9: 25,
    10: 8,
    11: 15,
    12: 8,
    13: 10,
    14: 20,
    15: 6,
    16: 17,
    17: 6,
    18: 20,
    19: 5,
    20: 8,
    21: 25,
    22: 11,
    23: 8,
    24: 7,
    25: 7,
    26: 28,
    27: 26,
    28: 27,
    29: 9,
    30: 24,
    31: 15,
    32: 14,
    33: 22,
    34: 20,
    35: 19,
    36: 8,
}

# Máximas toneladas a extraer de la zona i en M
n = {
    1: 189,
    2: 196,
    3: 143,
    4: 136,
    5: 106,
    6: 151,
    7: 170,
    8: 129,
    9: 184,
    10: 122,
    11: 146,
    12: 190,
    13: 160,
    14: 109,
    15: 133,
    16: 198,
    17: 138,
    18: 107,
    19: 117,
    20: 150,
    21: 171,
    22: 103,
    23: 157,
    24: 143,
    25: 170,
    26: 130,
    27: 140,
    28: 126,
    29: 180,
    30: 153,
    31: 108,
    32: 132,
    33: 105,
    34: 145,
    35: 145,
    36: 114,
}

# Mínimo de toneladas a explotar del tipo de carbón j en K
b = {
    "Antracita": 862,
    "Hulla": 898,
    "Turba": 562,
    "Lignito": 742,
}

Objeto del modelo#

Construye un problema al que luego agregarás las restricciones y la función objetivo.

problema = lp.LpProblem(name="extraccion_minera", sense=lp.LpMinimize)

Variables de decisión#

Define las variables del problema de manera que estén contenidas en diccionarios indexados en los conjuntos de sus variables respectivas.

# Toneladas a extraer de la zona i en M
x = {
    i: lp.LpVariable(
        name=f"extraccion_zona_{i}", lowBound=0, upBound=None, cat=lp.LpContinuous
    )
    for i in M
}

Función objetivo#

Agrega al problema la función objetivo. Recuerda que al definir el problema, ya definiste si este es de maximización o minimización.

problema += sum(c[i] * x[i] for i in M), "costo_total"

Restricciones#

Agrega al problema las restricciones del modelo.

# Garantiza el mínimo de toneladas extraidas del tipo de carbón j en K
for j in K:
    problema += (
        sum(x[i] for i in M if a[i, j] == 1) >= b[j],
        f"minimo_toneladas_extracción_carbon_{j}",
    )

# Garantiza el máximo de toneladas de carbón que se pueden extraer de la zona i en M
for i in M:
    problema += x[i] <= n[i], f"maximo_toneladas_extraccion_zona_{i}"

Resolver el problema#

Invoca el optimizador. Este paso le asigna un valor a las variables incluidas en las restricciones o función objetivo del modelo.

problema.solve()

Welcome to the CBC MILP Solver 
Version: 2.10.8 
Build Date: May  6 2022 

command line - cbc /tmp/0a4b410a84b6478793fadfda6c327b61-pulp.mps timeMode elapsed branch printingOptions all solution /tmp/0a4b410a84b6478793fadfda6c327b61-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 45 COLUMNS
At line 154 RHS
At line 195 BOUNDS
At line 196 ENDATA
Problem MODEL has 40 rows, 36 columns and 72 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Presolve 4 (-36) rows, 33 (-3) columns and 33 (-39) elements
0  Obj 11596.299 Primal inf 2370.5 (4)
5  Obj 32782
Optimal - objective value 32782
After Postsolve, objective 32782, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 32782 - 5 iterations time 0.002, Presolve 0.00
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

1

Imprimir resultados#

Antes de estudiar el óptimo del modelo, identifica en el estado del optimizador si pudo resolver el problema.

f"Estado del optimizador: {lp.LpStatus[problema.status]}"

'Estado del optimizador: Optimal'

Identifica también el valor de la función objetivo.

f"Costo total: {lp.value(problema.objective)}"

'Costo total: 32782.0'

Por último, imprime de manera estructurada el valor de las variables de decisión y otras expresiones de interés.

for i in M:
    if lp.value(x[i]) >= 0.1:
        for j in K:
            if a[i, j] == 1:
                print(f"{i}: {lp.value(x[i])} toneladas de {j}")

1: 189.0 toneladas de Hulla
2: 196.0 toneladas de Hulla
3: 103.0 toneladas de Turba
4: 136.0 toneladas de Turba
5: 106.0 toneladas de Hulla
7: 15.0 toneladas de Hulla
10: 122.0 toneladas de Antracita
11: 146.0 toneladas de Hulla
12: 190.0 toneladas de Turba
13: 39.0 toneladas de Antracita
15: 133.0 toneladas de Turba
17: 138.0 toneladas de Hulla
19: 117.0 toneladas de Antracita
20: 150.0 toneladas de Lignito
23: 157.0 toneladas de Antracita
24: 143.0 toneladas de Antracita
25: 170.0 toneladas de Antracita
28: 114.0 toneladas de Lignito
29: 180.0 toneladas de Lignito
30: 153.0 toneladas de Lignito
31: 108.0 toneladas de Hulla
35: 145.0 toneladas de Lignito
36: 114.0 toneladas de Antracita

Enunciado - Punto 2#

Ahora considera el escenario en que la empresa incurre en un costo fijo de \(q_i\) pesos cuando decide extraer carbón de la zona \(i\in M\). Por lo tanto, si se decide explotar la mina \(i\in M\), no se pueden extraer más de (\(n_i\)) toneladas de carbón. Pero si no se decide explotar, la extracción de carbón en esa zona debe ser igual a cero (0). La Tabla 3 presenta la información sobre los costos fijos de extraer cada zona.

Tabla 3. Costo fijo por extracción para cada zona

Zona

Costo fijo de extracción en la zona

1

240

2

155

3

240

4

125

5

177

6

342

7

157

8

457

9

396

10

411

11

341

12

469

13

402

14

186

15

404

16

344

17

290

18

340

19

482

20

472

21

394

22

102

23

330

24

433

25

205

26

394

27

156

28

298

29

134

30

462

31

432

32

362

33

127

34

203

35

417

36

215

Formulación#

a. Agrega a la formulación matemática los componentes de un modelo de optimización de forma general que represente la nueva situación. Defina clara y rigurosamente:

  • Conjuntos

  • Parámetros

  • Variables de decisión

  • Restricciones

  • Naturaleza de las variables

  • Función objetivo

Conjuntos#

  • \(M\): Zonas

  • \(K\): Tipos de carbón

Parámetros#

  • \(a_{ij}:\begin{cases}1\text{,}&\text{ si en la zona }i\in M\text{ se puede extraer carbón del tipo }j\in K \\ 0\text{,}& \text{ d.l.c} \end{cases}\)

  • \(c_i\): costo por cada tonelada de carbón extraída en la zona \(i\in M\)

  • \(n_i\): capacidad máxima de extracción de carbón (en toneladas) en la zona \(i\in M\)

  • \(b_j\): mínimo de toneladas a explotar del tipo de carbón \(j\in K\)

  • \(q_i\): costo fijo si se decide extraer carbón de la zona \(i\in M\)

Variables de decisión#

  • \(x_i\): toneladas de carbón extraído en la zona \(i\in M\)

  • \(y_i:\begin{cases} 1\text{,}&\text{si se decide explotar la zona }i\in M \\ 0\text{,}& \text{d.l.c} \end{cases}\)

Restricciones#

Se garantiza que la cantidad de carbón extraído del tipo \(j\in K\) sea como mínimo \(b_j\):

\[ \sum_{i\in M\vert a_{ij} = 1}x_i \ge b_j,\ \forall j\in K; \]

Se garantiza que no se puedan extraer más de \(n_i\) toneladas de carbón en la zona \(i\in M\), 0 si no se decide explotar la zona:

\[ x_i \le n_{i}y_{i},\ \forall i\in M; \]

Naturaleza de las Variables#

Solo pueden extraerse cantidades positivas de carbón de cada zona \(i\in M\):

\[ x_i \ge 0,\ \forall i\in M; \]

Cada zona \(i\in M\) puede o no ser explotada:

\[ y_{i} \in \{0, 1\},\ \forall i \in M; \]

Función objetivo#

Debe minimizarse el costo total de extracción y explotación:

\[ \operatorname{mín}\ \sum_{i\in M}c_{i}x_{i} + q_{i}y_{i}. \]

Implementación#

b. Resuelve el nuevo modelo planteado utilizando la librería de PuLP en Python. ¿Cuál es la solución óptima del problema?

Librerías#

Importa la librería pulp para crear y resolver el modelo.

import pulp as lp

Conjuntos#

Define los conjuntos M y K que representan las zonas y los tipos de carbón respectivamente.

Recuerda que por conveniencia de preservar el orden de los elementos de los conjuntos, no siempre deberás definirlos con el tipo set.

# Zonas
M = list(range(1, 37))

# Tipos de carbón
K = ["Antracita", "Hulla", "Turba", "Lignito"]

Parámetros#

Define o importa los parámetros del modelo.

# Si en la zona i en M se puede extraer carbón tipo j en K
a = {
    (1, "Hulla"): 1,
    (2, "Hulla"): 1,
    (3, "Turba"): 1,
    (4, "Turba"): 1,
    (5, "Hulla"): 1,
    (6, "Turba"): 1,
    (7, "Hulla"): 1,
    (8, "Hulla"): 1,
    (9, "Antracita"): 1,
    (10, "Antracita"): 1,
    (11, "Hulla"): 1,
    (12, "Turba"): 1,
    (13, "Antracita"): 1,
    (14, "Turba"): 1,
    (15, "Turba"): 1,
    (16, "Hulla"): 1,
    (17, "Hulla"): 1,
    (18, "Turba"): 1,
    (19, "Antracita"): 1,
    (20, "Lignito"): 1,
    (21, "Turba"): 1,
    (22, "Antracita"): 1,
    (23, "Antracita"): 1,
    (24, "Antracita"): 1,
    (25, "Antracita"): 1,
    (26, "Turba"): 1,
    (27, "Turba"): 1,
    (28, "Lignito"): 1,
    (29, "Lignito"): 1,
    (30, "Lignito"): 1,
    (31, "Hulla"): 1,
    (32, "Turba"): 1,
    (33, "Hulla"): 1,
    (34, "Antracita"): 1,
    (35, "Lignito"): 1,
    (36, "Antracita"): 1,
}

for i in M:
    for j in K:
        if (i, j) not in a:
            a[i, j] = 0

# Costo por tonelada de carbón en la zona i en M
c = {
    1: 16,
    2: 6,
    3: 11,
    4: 8,
    5: 5,
    6: 25,
    7: 16,
    8: 17,
    9: 25,
    10: 8,
    11: 15,
    12: 8,
    13: 10,
    14: 20,
    15: 6,
    16: 17,
    17: 6,
    18: 20,
    19: 5,
    20: 8,
    21: 25,
    22: 11,
    23: 8,
    24: 7,
    25: 7,
    26: 28,
    27: 26,
    28: 27,
    29: 9,
    30: 24,
    31: 15,
    32: 14,
    33: 22,
    34: 20,
    35: 19,
    36: 8,
}

# Máximas toneladas a extraer de la zona i en M
n = {
    1: 189,
    2: 196,
    3: 143,
    4: 136,
    5: 106,
    6: 151,
    7: 170,
    8: 129,
    9: 184,
    10: 122,
    11: 146,
    12: 190,
    13: 160,
    14: 109,
    15: 133,
    16: 198,
    17: 138,
    18: 107,
    19: 117,
    20: 150,
    21: 171,
    22: 103,
    23: 157,
    24: 143,
    25: 170,
    26: 130,
    27: 140,
    28: 126,
    29: 180,
    30: 153,
    31: 108,
    32: 132,
    33: 105,
    34: 145,
    35: 145,
    36: 114,
}

# Mínimo de toneladas a explotar del tipo de carbón j en K
b = {
    "Antracita": 862,
    "Hulla": 898,
    "Turba": 562,
    "Lignito": 742,
}

# Costo fijo de explotación de la zona i en M
q = {
    1: 240,
    2: 155,
    3: 240,
    4: 125,
    5: 177,
    6: 342,
    7: 157,
    8: 457,
    9: 396,
    10: 411,
    11: 341,
    12: 469,
    13: 402,
    14: 186,
    15: 404,
    16: 344,
    17: 290,
    18: 340,
    19: 482,
    20: 472,
    21: 394,
    22: 102,
    23: 330,
    24: 433,
    25: 205,
    26: 394,
    27: 156,
    28: 298,
    29: 134,
    30: 462,
    31: 432,
    32: 362,
    33: 127,
    34: 203,
    35: 417,
    36: 215,
}

Objeto del modelo#

Construye un problema al que luego agregarás las restricciones y la función objetivo.

problema = lp.LpProblem(name="extraccion_minera", sense=lp.LpMinimize)

Variables de decisión#

Define las variables del problema de manera que estén contenidas en diccionarios indexados en los conjuntos de sus variables respectivas.

# Toneladas a extraer de la zona i en M
x = {
    i: lp.LpVariable(
        name=f"extraccion_zona_{i}", lowBound=0, upBound=None, cat=lp.LpContinuous
    )
    for i in M
}

# Si se explota la zona i en M
y = {
    i: lp.LpVariable(
        name=f"binaria_zona_{i}", lowBound=0, upBound=None, cat=lp.LpBinary
    )
    for i in M
}

Función objetivo#

Agrega al problema la función objetivo. Recuerda que al definir el problema, ya definiste si este es de maximización o minimización.

problema += sum(c[i] * x[i] + q[i] * y[i] for i in M), "costo_total"

Restricciones#

Agrega al problema las restricciones del modelo.

# Garantiza el mínimo de toneladas extraidas del tipo de carbón j en K
for j in K:
    problema += (
        sum(x[i] for i in M if a[i, j] == 1) >= b[j],
        f"minimo_toneladas_extracción_carbon_{j}",
    )

# Garantiza el máximo de toneladas de carbón que se pueden extraer de la zona i en M
for i in M:
    problema += x[i] <= n[i] * y[i], f"maximo_toneladas_extraccion_zona_{i}"

Resolver el problema#

Invoca el optimizador. Este paso le asigna un valor a las variables incluidas en las restricciones o función objetivo del modelo.

problema.solve()

Welcome to the CBC MILP Solver 
Version: 2.10.8 
Build Date: May  6 2022 

command line - cbc /tmp/c9afc750df044a21ad207b14cb03a845-pulp.mps timeMode elapsed branch printingOptions all solution /tmp/c9afc750df044a21ad207b14cb03a845-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 45 COLUMNS
At line 298 RHS
At line 339 BOUNDS
At line 376 ENDATA
Problem MODEL has 40 rows, 72 columns and 108 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 39266.3 - 0.00 seconds
Cgl0004I processed model has 35 rows, 67 columns (31 integer (31 of which binary)) and 98 elements
Cbc0038I Initial state - 3 integers unsatisfied sum - 0.980279
Cbc0038I Pass   1: suminf.    0.21196 (1) obj. 40605.9 iterations 3
Cbc0038I Solution found of 40918
Cbc0038I Relaxing continuous gives 40328
Cbc0038I Before mini branch and bound, 27 integers at bound fixed and 17 continuous
Cbc0038I Full problem 35 rows 67 columns, reduced to 5 rows 16 columns
Cbc0038I Mini branch and bound improved solution from 40328 to 39488 (0.01 seconds)
Cbc0038I Round again with cutoff of 39468.7
Cbc0038I Reduced cost fixing fixed 11 variables on major pass 2
Cbc0038I Pass   2: suminf.    0.24243 (2) obj. 39468.7 iterations 3
Cbc0038I Pass   3: suminf.    0.30267 (1) obj. 39468.7 iterations 1
Cbc0038I Pass   4: suminf.    1.24478 (4) obj. 39468.7 iterations 4
Cbc0038I Pass   5: suminf.    0.81088 (3) obj. 39468.7 iterations 2
Cbc0038I Pass   6: suminf.    0.31260 (2) obj. 39468.7 iterations 2
Cbc0038I Pass   7: suminf.    0.80858 (4) obj. 39468.7 iterations 5
Cbc0038I Pass   8: suminf.    0.80858 (4) obj. 39468.7 iterations 0
Cbc0038I Pass   9: suminf.    0.83938 (3) obj. 39468.7 iterations 3
Cbc0038I Pass  10: suminf.    0.83938 (3) obj. 39468.7 iterations 1
Cbc0038I Pass  11: suminf.    0.62277 (3) obj. 39468.7 iterations 3
Cbc0038I Pass  12: suminf.    0.08054 (1) obj. 39468.7 iterations 3
Cbc0038I Pass  13: suminf.    0.53305 (2) obj. 39468.7 iterations 4
Cbc0038I Pass  14: suminf.    0.30386 (2) obj. 39468.7 iterations 4
Cbc0038I Pass  15: suminf.    0.24243 (2) obj. 39468.7 iterations 1
Cbc0038I Pass  16: suminf.    0.30267 (1) obj. 39468.7 iterations 1
Cbc0038I Pass  17: suminf.    0.31300 (2) obj. 39468.7 iterations 1
Cbc0038I Pass  18: suminf.    0.74938 (2) obj. 39468.7 iterations 2
Cbc0038I Pass  19: suminf.    0.80858 (4) obj. 39468.7 iterations 3
Cbc0038I Pass  20: suminf.    0.80858 (4) obj. 39468.7 iterations 0
Cbc0038I Pass  21: suminf.    0.83938 (3) obj. 39468.7 iterations 2
Cbc0038I Pass  22: suminf.    0.83938 (3) obj. 39468.7 iterations 0
Cbc0038I Pass  23: suminf.    0.62277 (3) obj. 39468.7 iterations 4
Cbc0038I Pass  24: suminf.    0.08054 (1) obj. 39468.7 iterations 4
Cbc0038I Pass  25: suminf.    0.53305 (2) obj. 39468.7 iterations 3
Cbc0038I Pass  26: suminf.    1.13452 (4) obj. 39468.7 iterations 6
Cbc0038I Pass  27: suminf.    1.13452 (4) obj. 39468.7 iterations 0
Cbc0038I Pass  28: suminf.    0.32169 (2) obj. 39468.7 iterations 6
Cbc0038I Pass  29: suminf.    0.83938 (3) obj. 39468.7 iterations 4
Cbc0038I Pass  30: suminf.    0.83938 (3) obj. 39468.7 iterations 0
Cbc0038I Pass  31: suminf.    0.62277 (3) obj. 39468.7 iterations 4
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 17 integers at bound fixed and 12 continuous
Cbc0038I Full problem 35 rows 67 columns, reduced to 14 rows 31 columns
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 39488 - took 0.02 seconds
Cbc0012I Integer solution of 39488 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0038I Full problem 35 rows 67 columns, reduced to 4 rows 18 columns
Cbc0031I 4 added rows had average density of 6
Cbc0013I At root node, 10 cuts changed objective from 39294.715 to 39488 in 1 passes
Cbc0014I Cut generator 0 (Probing) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 1 (Gomory) - 3 row cuts average 6.7 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 4 row cuts average 6.0 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is 1
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 3 row cuts average 6.7 elements, 0 column cuts (0 active)  in 0.000 seconds - new frequency is -100
Cbc0001I Search completed - best objective 39488, took 6 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 11 variables fixed on reduced cost
Cuts at root node changed objective from 39294.7 to 39488
Probing was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 1 times and created 4 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                39488.00000000
Enumerated nodes:               0
Total iterations:               6
Time (CPU seconds):             0.03
Time (Wallclock seconds):       0.02

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.03   (Wallclock seconds):       0.02

1

Imprimir resultados#

Antes de estudiar el óptimo del modelo, identifica en el estado del optimizador si pudo resolver el problema.

f"Estado del optimizador: {lp.LpStatus[problema.status]}"

'Estado del optimizador: Optimal'

Identifica también el valor de la función objetivo.

f"Costo total: {lp.value(problema.objective)}"

'Costo total: 39488.0'

Por último, imprime de manera estructurada el valor de las variables de decisión y otras expresiones de interés.

for i in M:
    if lp.value(x[i]) >= 0.1:
        for j in K:
            if a[i, j] == 1:
                print(f"{i}: {lp.value(x[i])} toneladas de {j}")

1: 189.0 toneladas de Hulla
2: 196.0 toneladas de Hulla
3: 103.0 toneladas de Turba
4: 136.0 toneladas de Turba
5: 106.0 toneladas de Hulla
7: 123.0 toneladas de Hulla
10: 122.0 toneladas de Antracita
11: 146.0 toneladas de Hulla
12: 190.0 toneladas de Turba
15: 133.0 toneladas de Turba
17: 138.0 toneladas de Hulla
19: 117.0 toneladas de Antracita
20: 150.0 toneladas de Lignito
22: 39.0 toneladas de Antracita
23: 157.0 toneladas de Antracita
24: 143.0 toneladas de Antracita
25: 170.0 toneladas de Antracita
28: 114.0 toneladas de Lignito
29: 180.0 toneladas de Lignito
30: 153.0 toneladas de Lignito
35: 145.0 toneladas de Lignito
36: 114.0 toneladas de Antracita

Créditos#

Equipo Principios de Optimización
Autores: Alfaima Solano, Alejandro Mantilla
Desarrollo:Alfaima Solano, Alejandro Mantilla, Juan Felipe Rengifo
Última fecha de modificación: 08/04/2023