Algoritmo Simplex aplicado ao problema do transporte

// 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 = 0;

// Tabela de custos
Cij = [20 15 10;
       12 08 16];

// Ofertas
oferta = [50;
          60];

// Demandas
demanda = [20;
           40;
           60];

// Método
// 1 = Noroeste
// 2 = Custo mínimo
metodo = 2;

// ============================================================================
// Não alterar mais a partir daqui
disp('=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-');
disp('Simplex para o problema do transporte');
disp('=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-');

// Transforma o problema de maximização em minimização
if tipo == 1 then
    Cij = max(Cij) - Cij;
end

// Verifica se a oferta é equivalente a demanda
if sum(oferta) > sum(demanda) then
    // Adiciona uma demanda virtual
    demanda(length(demanda)+1) = sum(oferta) - sum(demanda);
    // Adiciona mais uma coluna de custo zero
    Cij(:,length(demanda)) = zeros(length(oferta),1);
    // Mensagem
    disp('Detectado uma OFERTA > DEMANDA');
elseif sum(oferta) < sum(demanda) then
    // Adiciona uma oferta virtual
    oferta(length(oferta)+1) = sum(demanda) - sum(oferta);
    // Adiciona mais uma coluna de custo zero
    Cij(length(oferta),:) = zeros(1,length(demanda));
    disp('Detectado uma OFERTA < DEMANDA');
end

// Imprime a tabela
for k=1:30; mprintf('='); end
mprintf('\n Tabela de Custos:\n');
for k=1:(28+15*length(demanda)); mprintf('-'); end
mprintf('\n Destino \\ Origem |');
for k=1:length(demanda)
    mprintf(' Orig. (%3d)  |',k);
end
mprintf(' OFERTA |\n');
for k=1:(28+15*length(demanda)); mprintf('-'); end
for g = 1:length(oferta);
    mprintf('\n| Destino (%4d)  |',g);
    for k=1:length(demanda);
        mprintf(' %11.2f  |',Cij(g,k));
    end
    mprintf(' %6.2f |\n',oferta(g));
    for k=1:(28+15*length(demanda)); mprintf('-'); end
end
mprintf('\n| DEMANDA         |');
for k=1:length(demanda); mprintf(' %11.2f  |',demanda(k)); end
mprintf(' %d \\ %d |\n',sum(demanda),sum(oferta));
for k=1:(28+15*length(demanda)); mprintf('-'); end
mprintf('\n');

// Descobre o tamanho da matriz
[m,n] = size(Cij);

// Cria a matriz de solução
Xij = zeros(m,n);

// Inicialização do Transporte
if (metodo == 1) then
    disp('Inicialização pelo método do canto Noroeste: ');
    // Copia as variáveis temporariamente
    oferta_temp = oferta;
    demanda_temp = demanda;
    // Percorre todos os elementos
    for coluna = 1:n
        for linha = 1:m
            // Pecorre todas as linhas até não poder mais
            Xij(linha,coluna) = min(demanda_temp(coluna),oferta_temp(linha));
            // Atualiza as demandas e ofertas
            demanda_temp(coluna) = demanda_temp(coluna) - Xij(linha,coluna);
            oferta_temp(linha) = oferta_temp(linha) - Xij(linha,coluna);
        end
    end
elseif (metodo == 2) then
    disp('Inicialização pelo método do custo mínimo: ');
    // Copia as variáveis temporariamente
    Cij_temp = Cij;
    oferta_temp = oferta;
    demanda_temp = demanda;
    // Começa a otimização
    while(~isinf(min(Cij_temp)))
        [i,j] = min(Cij_temp);
        linha = j(1);
        coluna = j(2);
        // Atualiza o valor de Xij e Cij_temp
        Xij(linha,coluna) = min(demanda_temp(coluna),oferta_temp(linha));
        Cij_temp(linha,coluna) = %inf;
        // Atualiza as demandas e ofertas
        demanda_temp(coluna) = demanda_temp(coluna) - Xij(linha,coluna);
        oferta_temp(linha) = oferta_temp(linha) - Xij(linha,coluna);
    end
    if rank(Xij) < min(m,n) then
        error('Ciclo fechado');
    end
end

// Solução inicial
disp('Solução Inicial');
for k=1:30; mprintf('='); end
mprintf('\n Tabela de Custos:\n');
for k=1:(28+15*length(demanda)); mprintf('-'); end
mprintf('\n Destino \\ Origem |');
for k=1:length(demanda)
    mprintf(' Orig. (%3d)  |',k);
end
mprintf(' OFERTA |\n');
for k=1:(28+15*length(demanda)); mprintf('~'); end
for g = 1:length(oferta);
    mprintf('\n| Destino (%4d)  |',g);
    for k=1:length(demanda);
        mprintf(' %11.2f  |',Xij(g,k));
    end
    mprintf(' %6.2f |\n',oferta(g));
    for k=1:(28+15*length(demanda)); mprintf('-'); end
end
mprintf('\n| DEMANDA         |');
for k=1:length(demanda); mprintf(' %11.2f  |',demanda(k)); end
mprintf(' %d \\ %d |\n',sum(demanda),sum(oferta));
for k=1:(28+15*length(demanda)); mprintf('~'); end
mprintf('\n\n\n');


// Descobre as bases iniciais
nao_bases = (Xij' == 0);
nao_bases = find(nao_bases(:));

// Monta o sistema de equações
// A primeira coluna corresponde ao eixo z
z = Cij';
z = z(:)';
M = [1 -z 0];
// Adiciona o sistema de equações
// Equações de oferta
for linha = 1:m
    // Adiciona o primeiro elemento (z=0)
    M(1+linha,1) = 0;
    // Adiciona os demais coeficientes
    for coluna = 1:n
        M(1+linha,1+coluna + (linha-1)*n) = 1;
    end
    // Adiciona o valores de oferta
    M(1+linha,m*n+2) = oferta(linha);
end
// Equações de demanda
for coluna = 1:n
    // Adiciona o primeiro elemento (z=0)
    M(1+m+coluna,1) = 0;
    // Adiciona os demais coeficientes
    for linha = 1:m
        M(1+m+coluna,1+coluna + (linha-1)*n) = 1;
    end
    // Adicioa os valores de demanda
    M(1+m+coluna,m*n+2) = demanda(coluna);
end

// Elimia a ultima linha
M(m+n+1,:) = [];

// Obtém o novo tamanho da matriz
[m2, n2] = size(M);

// Cria o de novas bases com NaN
bases = zeros(1,m2-1);

// Não bases
nao_bases = Xij' == 0;
nao_bases = find(nao_bases(:));

// Repete até encontrar a solução ótima
// 1. Enquanto houver valores não-nulos nos coeficientes bases
// 2. Enquanto houver valores positivos nas variáveis básicas
// 3. A base ainda não está completa
while (sum([(sum(M(1,1+bases)) ~= 0) (max(M(1,1+nao_bases)) > 0) (sum(bases==0) > 0)])>0)

        // Preferencialmente, resolve-se as variáveis básicas (até descobrir todas)
        if (sum(bases==0) > 0) then
            v = 1:m*n;
            v([bases(bases>0) nao_bases]) = [];
            col = v;
        elseif (sum(M(1,1+bases)) ~= 0) then
            col = bases(M(1,1+bases) ~= 0);
        else
            col = nao_bases(M(1,1+nao_bases) ~= 0);
        end
        col = col(1);

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

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

        // Calcula os fatores limitantes
        q = M(2:m2,n2)./vdiv;

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

        // Copia fatores limitantes
        qnn = q;

        // Substitui fatores limitantes negativos por valores infinitos positivos
        qnn(q<=0) = %inf;

        // Procura pelo menor valor de q
        [q_min, lin] = min(qnn);

        // Se empatar no rank q, escolhe de forma aleatória
        if sum(qnn==q_min) > 1 then
              disp('Solução degenerada encontrada na próxima iteração.');
              v2 = find(qnn==q_min);
              lin = v2(round(rand()*(length(v2)-1))+1);
        end

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

       // Verifica se o elemento que vai sair da base existe
        if (bases(lin) > 0) then
            // Diz qual elemento vai saiu da base
            nao_bases(nao_bases==col) = lin;
        end

         // Preenche o elemento que vai entrar na base
        bases(lin) = col;

        // Operações
        for v = 1:(13+n2*10), mprintf('-');end;
            mprintf('\n | Fator Limitante | Base | %c |','z');
            for v = 1:n2-2, mprintf('   x%d   |',v);end;
            mprintf('   b    |\n |                 |      |');
            for v = 1:(-13+n2*10), mprintf('-');end;
            mprintf('\n |                 |      | 1 |');
            for v = 1:(n2-1), mprintf(' %6.2f |',M(1,1+v));end;
            mprintf('\n ');
            for v = 1:(13+n2*10), mprintf('-');end;
            mprintf('\n ');
            // ===========================================
            for v = 1:length(q)
                if(isinf(q(v))) then
                  mprintf('|        Infinito ');
                else
                    mprintf('| %15.2f ',q(v));
                end
                if isnan(bases(v)) then
                    mprintf('|  %2s  | 0 |','-');
                else
                    mprintf('|  x%d  | 0 |',bases(v));
                end

              for k = 1:(n2-1)
                mprintf(' %6.2f |',M(1+v,1+k));
              end
              mprintf(' \n ');
            end
            for v = 1:(13+n2*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');

            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:m2;      // Para todas as linhas
            atualizar(lin+1) = []; // Excerto a linha pivotal

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

            // Atualiza a linha pivotal por ultimo
            M(lin+1,:) = M(lin+1,:)./pivo;
            mprintf(' L%d = %.2f*L%d\n\n\n ',lin,1/pivo,lin);
end


// Operações
for v = 1:(13+n2*10), mprintf('-');end;
mprintf('\n | Fator Limitante | Base | %c |','z');
for v = 1:n2-2, mprintf('   x%d   |',v);end;
mprintf('   b    |\n |                 |      |');
for v = 1:(-13+n2*10), mprintf('-');end;
mprintf('\n |                 |      | 1 |');
for v = 1:(n2-1), mprintf(' %6.2f |',M(1,1+v));end;
mprintf('\n ');
for v = 1:(13+n2*10), mprintf('-');end;
mprintf('\n ');
// ===========================================
for v = 1:length(q)
    if(isinf(q(v))) then
        mprintf('|        Infinito ');
    else
        mprintf('| %15.2f ',q(v));
    end
    if isnan(bases(v)) then
        mprintf('|  %2s  | 0 |','-');
    else
        mprintf('|  x%d  | 0 |',bases(v));
    end
    for k = 1:(n2-1)
        mprintf(' %6.2f |',M(1+v,1+k));
    end
    mprintf(' \n ');
end
for v = 1:(13+n2*10), mprintf('-');end;
mprintf('\n');
// ===========================================