Estudos de otimização: Programação linear

Um dos tópicos de revisão de otimização que estou estudando é o de Programação Linear. Neste tipo de problema de otimização tanto a função objetivo quanto as restrições são funções lineares das variáveis.

Apesar de problemas deste tipo poderem ser expressos de diversas maneiras, a maneira abaixo, chamada de forma canônica, é frequentemente utilizada por materiais da área (e também por pacotes computacionais).

$$\begin{cases} \textrm{minimize} f(x) = c^Tx \textrm{ s.a.}\\ Ax = b\\ x \geq 0 \end{cases}$$

Conversão para a forma canônica

Seguem abaixo alguns ajustes que podem ser feitos para converter um problema linear qualquer para sua forma canônica:

1. Se problema for $maximize c^Tx$, converta-o para $minimize -c^Tx$;
2. Se houverem restrições de $\geq$ ou $\leq$, adicione uma variável $s_i$ para cada restrição e converta-as em restrições de igualdade da seguinte maneira: se $A_ix \geq b$, então $A_ix -s_i = b$; se $A_ix \leq b$, então $A_ix + s_i = b$.
3. Se houverem variáveis que podem ser negativas (sem restrição $x_i \geq 0$), faça a substituição $x_i = x_i^{(+)} - x_i^{(-)}, x_i^{(+)} \geq 0, x_i^{(-)} \geq 0$.

As transformações (1) e (3) são relativamente fáceis de serem entendidas. No caso da transformação (2), é possível atribuir a $s_i$ qualquer valor (pois s_i não faz parte da função objetivo). Desta maneira, se durante o proceso de otimização $A_ix = b$, então basta que $s_i = 0$.

O problema Dual

Todo problema de minimização (primal) em PL possui um problema de maximização (dual) relacionado de modo que ambos possuem o mesmo valor de função objetivo e é trivial encontrar a resposta de um a partir do outro.

Dado o problema linear na forma canônica acima, o problema dual consiste em:

$$\begin{cases} \textrm{maximize} f(y) = b^Ty \textrm{ s.a.}\\ A^Ty \leq c\\ \end{cases}$$

A partir da solução ótima $y^*$ do problema dual, é possível obter a solução ótima $x^*$ do problema primal.

Seja $u = c - A^Ty$. O vetor $u$ contem a "distância" que cada restrição dual está da igualdade. Uma restrição somente restringe a função objetivo se está "saturada", ou seja, $u_i = 0$. Desta maneira, podemos ignorar toda coluna de $A$ tal que $u_j \neq 0$, pois esta restrição não está ativa na solução do dual.Se o problem dual tiver solução única o número de componentes não zero de $u$ será igual ao número de linhas de $A$ (número de restrições primais). Logo, podemos obter $x$ resolvendo o seguinte sistema linear:

$$A[:,u_i = 0]x^* = b$$ .

As posições de x^* tal que $u_i \neq 0$ são completadas com $x^*_i = 0$, pois todas as restrições de $A$ já são satisfeitas com igualdade pelas outras componentes de $x^*_i$.

Exercício exemplo

Como exemplo, resolveremos um exercício simples de programação linear que ilustra os conceitos revistos.

Enunciado:

Uma startup busca fabricar máquinas de lavar falantes ao menor custo possível. Existem três maneiras de produzí-las: manualmente, semi automaticamente e automaticamente. A produção manual de uma máquina de lavar requer 1 minuto de trabalho qualificado, 40 minutos de trabalho não qualificado e 3 minutos de montagem. Para a produção semi automáticas os valores são, respectivamente, 4, 30 e 2 minutos e para a produção automática os valores são 8, 20 e 4 minutos. A startup dispões de 4500 minutos de trabalho qualificado, 36000 de trabalho não qualificado e 2700 minutos de montagem. Os custos de produção manual, semi automático e automático são, respectivamente, 70, 80 e 85 euros.

1. Formule o problema como um PL e dê uma solução.

No problema acima, as variáveis são o número $x_1, x_2 e x_3$ de máquinas feitas por cada tipo de produção (manual, semi-automática e automática). O custo a ser minimizado é, portanto $70x_1 + 80x_2 +85x_3$. Existem restrições quanto ao número de máquinas a serem produzidas e à capacidade de produção da fábrica (número de minutos para cada tipo de trabalho).

Abaixo o problema completo:

$$\begin{cases} \textrm{minimize} f(x) = 70x_1 + 80x_2 +85x_3 \textrm{ s.a.}\\ x_1 + x_2 + x_3 = 999\\ x_1 + 4x_2 + 8x_3 \leq 4500\\ 40x_1 + 30x_2 + 20x_3 \leq 36000\\ 3x_1 + 2x_2 + 4x_3 \leq 2700\\ x \geq 0 \end{cases}$$

Solução usando scipy:

import numpy as np
from scipy.optimize import linprog
from numpy.linalg import solve

A_eq = np.array([[1,1,1]])
b_eq = np.array([999])

A_ub = np.array([
[1, 4, 8],
[40,30,20],
[3,2,4]])

b_ub = np.array([4500, 36000,2700])

c = np.array([70, 80, 85])

res = linprog(c, A_eq=A_eq, b_eq=b_eq, A_ub=A_ub, b_ub=b_ub,
bounds=(0, None))
print('Valor otimo:', res.values()[3], '\nX:', res.values()[4])

Valor otimo: 73725.0
X: [ 636.  330.   33.]

2. Converta o problema para sua forma canônica

Para converter o problema acima para a forma canônica é necessário adicionar variáveis de folga para as restrições de desigualdade. Desta forma, o problema se torna:

$$\begin{cases} \textrm{minimize} f(x) = 70x_1 + 80x_2 +85x_3 \textrm{ s.a.}\\ x_1 + x_2 + x_3 = 999\\ x_1 + 4x_2 + 8x_3 + s_1 = 4500\\ 40x_1 + 30x_2 + 20x_3 + s_2 = 36000\\ 3x_1 + 2x_2 + 4x_3 + s_3 = 2700\\ x \geq 0, s \geq 0 \end{cases}$$

Solução usando scipy:

A = np.array([
[1, 1, 1, 0, 0, 0],
[1, 4, 8, 1, 0, 0],
[40, 30, 20, 0, 1, 0],
[3, 2, 4, 0, 0, 1]])

b = np.array([999, 4500, 36000, 2700])
c = np.array([70, 80, 85, 0, 0, 0])

res = linprog(c, A_eq=A, b_eq=b, bounds=(0, None))
print('Valor otimo:', res.values()[3], '\nX:', res.values()[4])

Valor otimo: 73725.0
X: [  636.   330.    33.  2280.     0.     0.]

Note que $s_1 = 2280$ implica que $x_1 + 4x_2 + 8x_3 \leq 4500$.

3. Apresente o dual do problema e mostre que ambas formulações obtem os mesmos resultados.

O problema dual correspondente é dado por:

$$\begin{cases} \textrm{maximize} f(y) = 999y_1 + 4500y_2 + 360000y_3 + 2700y_4 \textrm{ s.a.} \\ y_1 + y_2 + 40y_3 + 3y_4 \leq 70 \\ y_1 + 4y_2 + 30y_3 + 2y_4 \leq 80 \\ y_1 + 8y_2 + 20y_3 + 4y_4 \leq 85 \\ y_2 \leq 0 \\ y_3 \leq 0 \\ y_4 \leq 0 \end{cases}$$

Solução usando scipy: Para a solução usando $scipy$ podemos utilizar as variáveis declaradas anteriormente. Note que como linprog resolve problemas de minimização, o sinal da função objetivo está trocado.

res = linprog(-b, A_ub=A.T, b_ub=c, bounds=[(None,None), (None,None),
(None,None), (None,None)])
y = res.values()[4]
print('Valor otimo:', -res.values()[3], '\nY:', y)

u = c - A.T.dot(y)
Ar = A[:, np.abs(u)< 1e-10]
x_1 = solve(Ar, b)
print(x_1)
x = np.array([0.0] * len(c))
x[np.abs(u)< 1e-10] = x_1
x[np.abs(u)> 1e-10] = 0
print('Solucao X a partir de Y:', x)

Valor otimo: 73725.0
Y: [ 108.33333333    0.           -0.83333333   -1.66666667]
[  636.   330.    33.  2280.]
Solucao X a partir de Y: [  636.   330.    33.  2280.     0.     0.]