Algoritmo Simplex

// Universidade Federal do Rio Grande do Norte
// Centro de Tecnologia
// Programa de Pós-graduação em Engenharia Elétrica
// Aluno: Evandson Claude Seabra Dantas
// Matrícula: 20191003111

// Apaga as variáveis do Scilab e limpa a tela
clear; clc; funcprot(0);

// Tipo de problema
// 1 = Maximização
// 0 = Minimização
tipo = 1;

// Função objetivo
// z = x1 x2 x3 ...
z = [2 1];

// Restrições
xij = [1 0;
       0 1;
       1 3];

// Sinal das desilgualdades
// -1 representa <=
// 0 representa =
// +1 representa >=
s = [ -1;
      -1;
      +1];

// Sinal das desilgualdades
b = [+3;
     +4;
     +13];


// ==================================================================
// A PARTIR DAQUI NÃO PRECISA ALTERAR MAIS NADA
// ==================================================================
// ALGORITMO SIMPLEX (SEMPRE MAXIMIZAR)
// s_objetivo (caracter) : Letra correspondente a função objetivo
// s_M = Matriz expandida do Simplex
// s_imprime = 1 -> Imprime o passo-a-passo
// s_imprime = 0 -> Imprime apenas o resultado final
// H -> Matriz resolvida
// Bases -> Variáveis básicas encontradas
function [H, bases] = simplex(s_objetivo, s_M, bases, s_imprime)
    m = size(s_M, 1) - 1; // Exclui a linha da função objetivo
    n = size(s_M, 2) - 2; // Exclui a coluna z e a coluna b.
    p_col = -1;
    passo = 0;
    while(p_col<0)

        // Política gulosa (maior peso na função objetivo)
        [p_col, col] = min(s_M(1,2:n+1));

        // Vetor de divisão
        vdiv = s_M(2:m+1,col+1);

        // Troca os zeros por NaN (Not a Number)
        vdiv(vdiv==0) = %nan;

        // Calcula os fatores limitantes
        q = s_M(2:m+1,n+2)./vdiv;

        // Copia os fatores limitantes
        qnn = q;

        // Troca os NaN por Inf
        qnn(isnan(qnn)) = %inf;

        // Substitui fatores limitantes negativos por valores infinitos positivos
        qnn(q<0) = %inf;
        // Procura pelo menor valor de q

        [q_min, lin] = min(qnn);
        // Corrige a linha em função da primeira linha (função objetivo)

        // Separa o elemento pivo
        pivo = s_M(lin+1,col+1);

        // Imprime a matriz M
       // ===========================================
       if (s_imprime > 0)
            passo = passo + 1;
            mprintf('\n\n Passo %d:\n ',passo);
            for v = 1:(33+n*10), mprintf('-');end;
            mprintf('\n | Fator Limitante | Base | %c |',s_objetivo);
            for v = 1:n, mprintf('   x%d   |',v);end;
            mprintf('   b    |\n |                 |      |');
            for v = 1:(7+n*10), mprintf('-');end;
            mprintf('\n |                 |      | 1 |');
            for v = 1:(n+1), mprintf(' %6.2f |',s_M(1,1+v));end;
            mprintf('\n ');
            for v = 1:(33+n*10), mprintf('-');end;
            mprintf('\n ');
            // ===========================================
            for v = 1:length(q)
                if(isnan(q(v))) then
                  mprintf('|        Infinito |  x%d  | 0 |',bases(v));
                else
                  mprintf('| %15.2f |  x%d  | 0 |',q(v),bases(v));
                end
              for k = 1:(n+1)
                mprintf(' %6.2f |',s_M(1+v,1+k));
              end
              mprintf(' \n ');
            end
            for v = 1:(33+n*10), mprintf('-');end;
            mprintf('\n');
            // ===========================================
            mprintf(' Coluna Pivotal: %d\n',col);
            mprintf(' Linha Pivotal: %d\n',lin);
            mprintf(' Pivô: %.2f\n\n',pivo);
            mprintf(' Operações:\n');
        end

        if (isinf(q_min)) then
            disp(' Não há convergência para estes parâmetros! ');
            abort;
        end

        // Lista as linhas a serem atualizadas
        atualizar = 1:m+1;   // Para todas as linhas
        atualizar(lin+1) = []; // Excerto a linha pivotal

        // Para cada linha da tabela...
        for i = atualizar
           // Calcula o escalar
           alpha = s_M(i,col+1)./pivo;
           s_M(i,:) = s_M(i,:) - alpha*s_M(lin+1,:);
           if (s_imprime > 0)
            mprintf(' L%d = L%d - %.2f*L%d\n',i-1,i-1,alpha,lin);
           end
        end

        // Atualiza a linha pivotal por ultimo
        s_M(lin+1,:) = s_M(lin+1,:)./pivo;
        if (s_imprime > 0)
            mprintf(' L%d = %.2f*L%d\n\n\n ',lin,1/pivo,lin);
        end
        // Troca as bases
        bases(lin) = col;

    end
    // Copia a matriz resolvida
    H = s_M;
endfunction

// Verifica se existem erros de digitação
if (2*size(xij,1)~=length(s)+length(b)) then
    error('O número de linhas, sinais e desigualdades não coincidem');
end

// ==================================================================
// Imprime o sistema de equações
mprintf('Função Objetivo: \n ');
if (tipo>0) then
    mprintf('Maximizar ');
else
    mprintf('Minimizar ');
end
mprintf(' z = ');
for i = 1:length(z)
    if(z(i)~=0) then
      mprintf(' + (%.2f)x%d',z(i),i);
    end
end
mprintf('\n\n');
for lin = 1:size(xij,1)
    for col = 1:size(xij,2)
        if(xij(lin,col)~=0) then
            mprintf(' + (%.2f)x%d',xij(lin,col),col);
        end
    end
    if (s(lin)<0) then
        mprintf(' ≤ ');
    elseif (s(lin)>0) then
        mprintf(' ≥ ');
    else
        mprintf(' = ');
    end
     mprintf('%.2f\n',b(lin));
end

// Fim da impressão do sistema de equações
// ==================================================================

// Tira qualquer negatividade do segundo membro da equação
for lin = 1:length(s)
    if (b(lin) < 0) then
       b(lin) = -b(lin);
       s(lin) = -s(lin); // troca o sinal da desigualdade
       xij(lin,:) = -xij(lin,:);
    end
end

// Copia a matriz
M = xij;

// Trata as restrições de sinais (expandindo a matriz)
// Armazena as posições das variáveis básicas
bases = [];
// Armazena as posições das variáveis artificiais
v_artificiais = [];
for i = 1:length(s)
    if (s(i) < 0) then
        // Adiciona uma variável de folga
        M(i,size(M,2)+1) = 1;
        // Adiciona a posição da base
        bases = [bases size(M,2)];
    elseif (s(i) == 0)
        // Adiciona uma variável artificial
        M(i,size(M,2)+1) = 1;
        // Salva a posição da variável artificial
        v_artificiais = [v_artificiais size(M,2)];
        // Adiciona a posição da base
        bases = [bases size(M,2)];
    elseif (s(i) > 0)
        // Adiciona uma variável de folta e uma artificial
        M(i,size(M,2)+[1 2]) = [-1 1];
        // Salva a posição da variável artificial
        v_artificiais = [v_artificiais size(M,2)];
        // Adiciona a posição da base
        bases = [bases size(M,2)];
    end
end

// Obtém m e n
m = size(M, 1);
n = size(M, 2);

// Completa a função objetivo com zeros (variáveis básicas)
z = [z zeros(1,n-length(z))];

// Após adicionar as variáveis de folga e artificiais a super matriz
// Adiciona-se a igualdade
M = [M b];
// Adiciona a coluna referente ao eixo Z
M = [zeros(m,1) M];
// Adiciona a linha referente ao eixo Z trantado o tipo (maximizar ou minimizar)
M = [1 ((-1)^(tipo>0)).*z 0;M];

// Relatório Inicial
mprintf(' =======================\n');
mprintf('  Sistema de Equações\n');
mprintf(' =======================\n');
for v = 1:(m+1)
    mprintf(' (%2.f)z',M(v,1));
    for k = 2:(n+1)
        mprintf(' + (%2.f)x%d',M(v,k),k-1);
    end
    mprintf(' = (%2.f)\n',M(v,n+2));
end
mprintf('\n ');

// Verifica se vai logo para o Simplex ou tem variáveis auxiliares
// para resolver antes
if (v_artificiais~=[]) then
    disp('Método da Duas Etapas primeiro!');
    disp('Objetivo: Zerar as variáveis artificiais');
    disp('Objetivo: Minimizar w: ');

    // Cria a nova função objetivo
    w = zeros(1,n);
    w(v_artificiais) = 1;

   // Imprime a nova função objetivo
    mprintf(' w = ');
    for i = 1:length(w)
        if(w(i)~=0) then
          mprintf(' + (%.2f)x%d',w(i),i);
        end
    end
    mprintf('\n ');

    // Copia a mega-matriz
    M2 = M;

    // Substitui a função objetivo original pela nova função objetivo
    M2(1,:) = [1 w 0];

    // Em função da variável não básica
    M2(1,:) = M2(1,:) - M2(m+1,:);

    // Roda o SIMPLEX para minimizar w
    [M2, bases] = simplex('w', M2, bases, 1);

    // Substitui a função objetivo pela original
    M2(1,:) = M(1,:);

    // Exclui as colunas pivotais
    M2(:,v_artificiais+1) = [];

    // Atualiza a super matriz original (sem as variáveis artificiais)
    M = M2;

    // Atualiza n (removendo as variáveis artificiais)
    n = size(M, 2) - 2;

end
disp(' ');
disp('Simplex!');
disp('===================================================');
[M, bases] = simplex('z',M, bases, 1);
disp('===================================================');
disp('Resultados ótimos encontrados: ');

// Imprime a melhor solução
// Verifica se existem múltiplas soluções
 if (tipo>0) then
     mprintf('\n\tMax z = %.2f\n\n',M(1,n+2));
 else
     mprintf('\n\tMin z = %.2f\n\n',M(1,n+2));
 end
 mprintf(' Com variáveis básicas: \n');
 for i = 1:length(bases)
     mprintf('\tx%d = %.2f\n',bases(i),M(i+1,n+2));
 end

// Procura pelas variáveis NÃO-BÁSICAS
nao_bases = 1:n;
// Devido a um bug do Scilab temos que usar bases' se não a variável
// bases muda o seu valor.
nao_bases(bases') = [];
if sum(M(1,1+nao_bases)==0) then
    disp('Opa! Existem outras soluções para este problema.');
    disp('-------------------------------');
    var_nao_bases_nulas = nao_bases(M(1,1+nao_bases)==0);
    for outra_base = var_nao_bases_nulas
        // Subtrai a linha 0 pela linha que contém maior elemento na coluna da variável não-basica
        [p_linha, linha] = max(M(2:m+1,outra_base));
        // Trocando a Linha X pela Coluna Y
        mprintf(' Coluna Pivotal: %d\n ',outra_base);
        mprintf('Linha Pivotal: %d\n ',linha);
        mprintf('Pivô: %d\n ',M(1+linha,1+outra_base));
        mprintf('Nova solução: ');
        atualizar = 1:m+1;   // Para todas as linhas
        atualizar(linha+1) = []; // Excerto a linha pivotal
        for i = atualizar
           // Calcula o escalar
           alpha = M(i,outra_base+1)./M(1+linha,1+outra_base);
           M(i,:) = M(i,:) - alpha*M(linha+1,:);
        end
        // atualiza a linha pivotal
        M(linha+1,:) = M(linha+1,:)./M(1+linha,1+outra_base);
        // Atualiza as bases
        bases(linha) = outra_base;

        if (tipo>0) then
         mprintf('\n\tMax z = %.2f\n\n',M(1,n+2));
         else
             mprintf('\n\tMin z = %.2f\n\n',M(1,n+2));
         end
         mprintf(' Com variáveis básicas: \n');
         for i = 1:length(bases)
             mprintf('\tx%d = %.2f\n',bases(i),M(i+1,n+2));
         end
    end

end