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:
- Se problema for \(maximize c^Tx\), converta-o para \(minimize -c^Tx\);
- 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\).
- 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.]
