R Notebook

BIBLIOTECAS

library(magrittr)
library(combinat)

CONFIGURAÇÕES INICIAIS

‘A’ é a matriz de restrições, de tamanho m x n (m = nº de restrições e n = número de variáveis, na forma padrão). ‘b’ é o vetor right hand side, conjunto de valores numéricos do lado direito da igualdade. ‘f_obj’ é o vetor de coeficientes da função objetivo, chamado de c na teoria. Neste exemplo, o sistema tem duas incógnitas, descontadas as incógnitas de folga, e três restrições (questão 4a).

# nomes das variáveis
var_names <- c("x1","x2","xf1","xf2","xf3")

# montagem de A
A <- rbind(c(1,0,1,0,0),
           c(0,1,0,1,0),
           c(4,-2,0,0,1))
colnames(A) <- var_names

# montagem de b
b <- matrix(c(6,8,10),ncol=1)

# montagem de f_obj
f_obj <- c(-4,-7,0,0,0)

# m: nº de restrições
# n: nº de variáveis
m <- nrow(A)
n <- ncol(A)

cat("Matriz A: \n");print(A);cat("\n")

Matriz A: 
     x1 x2 xf1 xf2 xf3
[1,]  1  0   1   0   0
[2,]  0  1   0   1   0
[3,]  4 -2   0   0   1

cat("Vetor b: \n");print(b);cat("\n")

Vetor b: 
     [,1]
[1,]    6
[2,]    8
[3,]   10

cat("Vetor c: \n");print(f_obj);cat("\n")

Vetor c: 
[1] -4 -7  0  0  0

DEFINIÇÃO DA MATRIZ ‘B’ E MULTIPLICADOR SIMPLEX

A matriz ‘B’ representa uma matriz LI de um sistema do tipo B.x = b tal que x é uma solução factível (i.e., um vetor de elementos não-negativos). Portanto, a escolha de ‘B’ deve satisfazer: a) det(B) != 0 e b) (B^-1).b não tem elementos negativos. O código a seguir testa as permutações do conjunto {1,2,3,…,n} das colunas de ‘A’ até encontrar ‘B’ adequada:

# todas as combinações possíveis de m colunas de A 
comb = combn(1:n,m)
# busca de uma matriz básica factível
  # processo dura, no máximo, até o nº de permutações
  for (i in 1:ncol(comb)){
    # estimativa inicial da matriz básica
    B <- A[,comb[,i]]
    
    # condição 1: det(B) != 0
    if (det(B) != 0){
      xB = inv(B) %*% b
      
      # condição 2: não-negatividade da solução básica inicial
      if (sum(xB<0)==0) {
        i_sel <- i
        # encerramento da busca
        break
      }
    }
  }
cat("Combinações: \n");print(comb);cat("\n")

Combinações: 
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    1    1    1    1    2    2    2     3
[2,]    2    2    2    3    3    4    3    3    4     4
[3,]    3    4    5    4    5    5    4    5    5     5

cat("Nº de iterações para achar 'B': ");cat(i_sel);cat("\n\n")

Nº de iterações para achar 'B': 2

cat("Matriz básica inicial: \n")

Matriz básica inicial: 

print(B);cat("\n")

     x1 x2 xf2
[1,]  1  0   0
[2,]  0  1   1
[3,]  4 -2   0

cat("Solução básica inicial: \n")

Solução básica inicial: 

print(xB)

    [,1]
x1     6
x2     7
xf2    1

Como foi achada na segunda iteração, a matriz ‘B’ corresponde ao índices da segunda coluna da matriz de permutações (colunas 1, 2 e 4 de A). De posse da matriz básica, é hora de calcular o multiplicador simplex, simplesmente multiplicando os coeficientes da função objetivo correspondentes às variáveis básicas (aquelas de índice 1, 2 e 4, correspondentes aos índices das colunas que construíram ‘B’. No nosso exemplo, corresponde às variáveis x1, x2 e xf2):

# coeficientes de custo relativo das variáveis básicas  
  cB = f_obj[comb[,i_sel]]
  
  # multiplicador simplex
  lambdaT = matrix(cB,nrow = 1) %*% inv(B)
  
  # exibição do multiplicador da iteração
  cat("Multiplicador simplex:\n");print(lambdaT);cat("\n")

Multiplicador simplex:
     [,1] [,2] [,3]
[1,]  -18    0  3.5

CUSTOS RELATIVOS

Nos interessam os custos associados às variáveis nulas (como ‘xB’ foi encontrado considerando que as variáveis remanescentes - 3 e 5, xf1 e xf2 no nosso exemplo - são nulas, tudo relativo a elas recebe o adjetivo “nulo”). Esse custo é simplesmente os coeficientes da função objetivo associados às variáveis não-básicas (nulas), de índices 3 e 5, menos lambdaT.N, onde N é formado pelas colunas 3 e 5 de A. lambdaT é o multiplicador simplex do passo anterior:

# índice das variáveis nulas
  sel = 1:ncol(A)
  sel = sel[-comb[,i_sel]]
  
  # coeficientes de custo relativo das variáveis nulas  
  cN <- f_obj[sel] - lambdaT %*% A[,sel]
  
  # exibição dos coeficientes de custo relativo das variáveis nulas  
  cat("Índices das variáveis nulas:\n");print(sel);cat("\n")

Índices das variáveis nulas:
[1] 3 5

  cat("Custos relativos cN:\n");print(cN);cat("\n")

Custos relativos cN:
     xf1  xf3
[1,]  18 -3.5

CONDIÇÕES DE PARADA E OTIMIZAÇÃO DA SOLUÇÃO INICIAL

O algoritmo encerra quando todos os custos da variáveis nulas são maiores que 0 (para o caso de minimização). No nosso exemplo, isto não se configura, e já sabemos que haverá outra iteração. Caso tivéssemos cN >= 0 em todos os elementos, ‘xB’ já seria ótima.

A otimização começa pela determinação da variável de mínimo custo, a de índice 5 (xf3) no nosso exemplo. a direção simplex é o produto da inversa de B pela coluna 5, equivalente à variável de menor custo:

# índice do mínimo custo relativo de cN
  k = which(cN == min(cN,na.rm = T))
# direção simplex
    y = inv(B) %*% matrix(A[,sel[k]],ncol=1)
    cat("Direção simplex y:\n");print(y);cat("\n")

Direção simplex y:
    [,1]
x1   0.0
x2  -0.5
xf2  0.5

Se não houvesse elemento de ‘y’ estritamente positivo (>0), o problema seria ilimitado; no nosso exemplo, temos apenas o terceiro elemento de y satisfazendo essa condição. Com todos os elementos de ‘y’ maiores que zero, construímos o vetor ‘eps’ composto pelos quocientes entre ‘xB’ e ‘y’ (novamente, somente os elementos de ‘xB’ que correspondem aos elementos de ‘y’ maiores que 0):

# seleção de y>0
      sel_y <- y>0
      # passo
      eps = xB[sel_y]/y[sel_y]
cat("Passo mínimo: ");cat(eps)

Passo mínimo: 2

Agora selecionamos o menor passo mínimo, trivial no nosso exemplo, posto que só há um elemento. Esse elemento corresponde à variável básica dessa iteração que será substituída na próxima:

# índice do passo mínimo
      l = which(eps == min(eps))
      
      # seleção da variável substituída
      sel2 <- (((1:n)[-sel])[sel_y])[l]
cat("Variável substituída: ");cat(sel2);cat(", ");cat(var_names[sel2])

Variável substituída: 4, xf2

Para que xf2 possa sair, alguma terá que entrar, e será aquela com o menor custo mínimo dentre as variáveis nulas:

# coluna que vai sair da base
      L = A[,sel2]
      nL <- var_names[sel2]
# coluna que vai entrar da base
      K = A[,sel[k]]
      nK <- var_names[sel[k]]
      cat("Variável de substituição: ");cat(sel[k]);cat(", ");cat(nK)

Variável de substituição: 5, xf3

Agora, resta apenas proceder à substituição e recomeçar o processo, até que não haja mais custos de variáveis nulas negativos:

# substituição das colunas
      A[,sel[k]] = L
      A[,sel2] = K
      colnames(A)[sel[k]] <- nL
      colnames(A)[sel2] <- nK
      
      cat("Nova matriz A:\n");print(A);cat("\n---------------------\n\n")

Nova matriz A:
     x1 x2 xf1 xf3 xf2
[1,]  1  0   1   0   0
[2,]  0  1   0   0   1
[3,]  4 -2   0   1   0

---------------------

Note que houve a inversão na prioridade de busca entre ‘xf2’ e ‘xf3’. Como o sistema implementado aqui é sequencial, as colunas 1, 2 e 3 de ‘A’ serão descartadas, como na iteração anterior, e, de novo, as colunas 1, 2 e 4 serão selecionadas, mas agora com ‘xf3’ ocupando a 4ª posição. Num cálculo manual, bastaria substituir a última coluna de ‘B’, (0 1 0)‘, por (0 0 1)’ e executar os passos. A partir daqui, volta-se à seção # DEFINIÇÃO DA MATRIZ ‘B’ E MULTIPLICADOR SIMPLEX.