SHARE
TWEET

Algoritmo Simplex aplicado ao problema do transporte

evnds Apr 4th, 2019 103 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // Universidade Federal do Rio Grande do Norte
  2. // Centro de Tecnologia
  3. // Programa de Pós-graduação em Engenharia Elétrica
  4. // Aluno: Evandson Claude Seabra Dantas
  5. // Matrícula: 20191003111
  6.  
  7. // Apaga as variáveis do Scilab e limpa a tela
  8. clear; clc; funcprot(0);
  9.  
  10. // Tipo de problema
  11. // 1 = Maximização
  12. // 0 = Minimização
  13. tipo = 0;
  14.  
  15. // Tabela de custos
  16. Cij = [20 15 10;
  17.        12 08 16];
  18.  
  19. // Ofertas
  20. oferta = [50;
  21.           60];
  22.  
  23. // Demandas
  24. demanda = [20;
  25.            40;
  26.            60];
  27.  
  28. // Método
  29. // 1 = Noroeste
  30. // 2 = Custo mínimo
  31. metodo = 2;
  32.  
  33. // ============================================================================
  34. // Não alterar mais a partir daqui
  35. disp('=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-');
  36. disp('Simplex para o problema do transporte');
  37. disp('=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-');
  38.  
  39. // Transforma o problema de maximização em minimização
  40. if tipo == 1 then
  41.     Cij = max(Cij) - Cij;
  42. end
  43.  
  44. // Verifica se a oferta é equivalente a demanda
  45. if sum(oferta) > sum(demanda) then
  46.     // Adiciona uma demanda virtual
  47.     demanda(length(demanda)+1) = sum(oferta) - sum(demanda);
  48.     // Adiciona mais uma coluna de custo zero
  49.     Cij(:,length(demanda)) = zeros(length(oferta),1);
  50.     // Mensagem
  51.     disp('Detectado uma OFERTA > DEMANDA');
  52. elseif sum(oferta) < sum(demanda) then
  53.     // Adiciona uma oferta virtual
  54.     oferta(length(oferta)+1) = sum(demanda) - sum(oferta);
  55.     // Adiciona mais uma coluna de custo zero
  56.     Cij(length(oferta),:) = zeros(1,length(demanda));
  57.     disp('Detectado uma OFERTA < DEMANDA');
  58. end
  59.  
  60. // Imprime a tabela
  61. for k=1:30; mprintf('='); end
  62. mprintf('\n Tabela de Custos:\n');
  63. for k=1:(28+15*length(demanda)); mprintf('-'); end
  64. mprintf('\n Destino \\ Origem |');
  65. for k=1:length(demanda)
  66.     mprintf(' Orig. (%3d)  |',k);
  67. end
  68. mprintf(' OFERTA |\n');
  69. for k=1:(28+15*length(demanda)); mprintf('-'); end
  70. for g = 1:length(oferta);
  71.     mprintf('\n| Destino (%4d)  |',g);
  72.     for k=1:length(demanda);
  73.         mprintf(' %11.2f  |',Cij(g,k));
  74.     end
  75.     mprintf(' %6.2f |\n',oferta(g));
  76.     for k=1:(28+15*length(demanda)); mprintf('-'); end
  77. end
  78. mprintf('\n| DEMANDA         |');
  79. for k=1:length(demanda); mprintf(' %11.2f  |',demanda(k)); end
  80. mprintf(' %d \\ %d |\n',sum(demanda),sum(oferta));
  81. for k=1:(28+15*length(demanda)); mprintf('-'); end
  82. mprintf('\n');
  83.  
  84. // Descobre o tamanho da matriz
  85. [m,n] = size(Cij);
  86.  
  87. // Cria a matriz de solução
  88. Xij = zeros(m,n);
  89.  
  90. // Inicialização do Transporte
  91. if (metodo == 1) then
  92.     disp('Inicialização pelo método do canto Noroeste: ');
  93.     // Copia as variáveis temporariamente
  94.     oferta_temp = oferta;
  95.     demanda_temp = demanda;
  96.     // Percorre todos os elementos
  97.     for coluna = 1:n
  98.         for linha = 1:m
  99.             // Pecorre todas as linhas até não poder mais
  100.             Xij(linha,coluna) = min(demanda_temp(coluna),oferta_temp(linha));
  101.             // Atualiza as demandas e ofertas
  102.             demanda_temp(coluna) = demanda_temp(coluna) - Xij(linha,coluna);
  103.             oferta_temp(linha) = oferta_temp(linha) - Xij(linha,coluna);
  104.         end
  105.     end
  106. elseif (metodo == 2) then
  107.     disp('Inicialização pelo método do custo mínimo: ');
  108.     // Copia as variáveis temporariamente
  109.     Cij_temp = Cij;
  110.     oferta_temp = oferta;
  111.     demanda_temp = demanda;
  112.     // Começa a otimização
  113.     while(~isinf(min(Cij_temp)))
  114.         [i,j] = min(Cij_temp);
  115.         linha = j(1);
  116.         coluna = j(2);
  117.         // Atualiza o valor de Xij e Cij_temp
  118.         Xij(linha,coluna) = min(demanda_temp(coluna),oferta_temp(linha));
  119.         Cij_temp(linha,coluna) = %inf;
  120.         // Atualiza as demandas e ofertas
  121.         demanda_temp(coluna) = demanda_temp(coluna) - Xij(linha,coluna);
  122.         oferta_temp(linha) = oferta_temp(linha) - Xij(linha,coluna);
  123.     end
  124.     if rank(Xij) < min(m,n) then
  125.         error('Ciclo fechado');
  126.     end
  127. end
  128.  
  129. // Solução inicial
  130. disp('Solução Inicial');
  131. for k=1:30; mprintf('='); end
  132. mprintf('\n Tabela de Custos:\n');
  133. for k=1:(28+15*length(demanda)); mprintf('-'); end
  134. mprintf('\n Destino \\ Origem |');
  135. for k=1:length(demanda)
  136.     mprintf(' Orig. (%3d)  |',k);
  137. end
  138. mprintf(' OFERTA |\n');
  139. for k=1:(28+15*length(demanda)); mprintf('~'); end
  140. for g = 1:length(oferta);
  141.     mprintf('\n| Destino (%4d)  |',g);
  142.     for k=1:length(demanda);
  143.         mprintf(' %11.2f  |',Xij(g,k));
  144.     end
  145.     mprintf(' %6.2f |\n',oferta(g));
  146.     for k=1:(28+15*length(demanda)); mprintf('-'); end
  147. end
  148. mprintf('\n| DEMANDA         |');
  149. for k=1:length(demanda); mprintf(' %11.2f  |',demanda(k)); end
  150. mprintf(' %d \\ %d |\n',sum(demanda),sum(oferta));
  151. for k=1:(28+15*length(demanda)); mprintf('~'); end
  152. mprintf('\n\n\n');
  153.  
  154.  
  155. // Descobre as bases iniciais
  156. nao_bases = (Xij' == 0);
  157. nao_bases = find(nao_bases(:));
  158.  
  159. // Monta o sistema de equações
  160. // A primeira coluna corresponde ao eixo z
  161. z = Cij';
  162. z = z(:)';
  163. M = [1 -z 0];
  164. // Adiciona o sistema de equações
  165. // Equações de oferta
  166. for linha = 1:m
  167.     // Adiciona o primeiro elemento (z=0)
  168.     M(1+linha,1) = 0;
  169.     // Adiciona os demais coeficientes
  170.     for coluna = 1:n
  171.         M(1+linha,1+coluna + (linha-1)*n) = 1;
  172.     end
  173.     // Adiciona o valores de oferta
  174.     M(1+linha,m*n+2) = oferta(linha);
  175. end
  176. // Equações de demanda
  177. for coluna = 1:n
  178.     // Adiciona o primeiro elemento (z=0)
  179.     M(1+m+coluna,1) = 0;
  180.     // Adiciona os demais coeficientes
  181.     for linha = 1:m
  182.         M(1+m+coluna,1+coluna + (linha-1)*n) = 1;
  183.     end
  184.     // Adicioa os valores de demanda
  185.     M(1+m+coluna,m*n+2) = demanda(coluna);
  186. end
  187.  
  188. // Elimia a ultima linha
  189. M(m+n+1,:) = [];
  190.  
  191. // Obtém o novo tamanho da matriz
  192. [m2, n2] = size(M);
  193.  
  194. // Cria o de novas bases com NaN
  195. bases = zeros(1,m2-1);
  196.  
  197. // Não bases
  198. nao_bases = Xij' == 0;
  199. nao_bases = find(nao_bases(:));
  200.  
  201. // Repete até encontrar a solução ótima
  202. // 1. Enquanto houver valores não-nulos nos coeficientes bases
  203. // 2. Enquanto houver valores positivos nas variáveis básicas
  204. // 3. A base ainda não está completa
  205. while (sum([(sum(M(1,1+bases)) ~= 0) (max(M(1,1+nao_bases)) > 0) (sum(bases==0) > 0)])>0)
  206.      
  207.         // Preferencialmente, resolve-se as variáveis básicas (até descobrir todas)
  208.         if (sum(bases==0) > 0) then
  209.             v = 1:m*n;
  210.             v([bases(bases>0) nao_bases]) = [];
  211.             col = v;
  212.         elseif (sum(M(1,1+bases)) ~= 0) then
  213.             col = bases(M(1,1+bases) ~= 0);
  214.         else
  215.             col = nao_bases(M(1,1+nao_bases) ~= 0);
  216.         end
  217.         col = col(1);
  218.                
  219.         // Vetor de divisão
  220.         vdiv = M(2:m2,col+1);
  221.        
  222.         // Troca os zeros por NaN (Not a Number)
  223.         vdiv(vdiv==0) = %nan;
  224.        
  225.         // Calcula os fatores limitantes
  226.         q = M(2:m2,n2)./vdiv;
  227.        
  228.         // Troca os NaN por Inf
  229.         q(isnan(q)) = %inf;
  230.        
  231.         // Copia fatores limitantes
  232.         qnn = q;
  233.        
  234.         // Substitui fatores limitantes negativos por valores infinitos positivos
  235.         qnn(q<=0) = %inf;
  236.        
  237.         // Procura pelo menor valor de q
  238.         [q_min, lin] = min(qnn);
  239.        
  240.         // Se empatar no rank q, escolhe de forma aleatória
  241.         if sum(qnn==q_min) > 1 then
  242.               disp('Solução degenerada encontrada na próxima iteração.');
  243.               v2 = find(qnn==q_min);
  244.               lin = v2(round(rand()*(length(v2)-1))+1);
  245.         end
  246.            
  247.         // Separa o elemento pivo
  248.         pivo = M(lin+1,col+1);
  249.        
  250.        // Verifica se o elemento que vai sair da base existe
  251.         if (bases(lin) > 0) then
  252.             // Diz qual elemento vai saiu da base
  253.             nao_bases(nao_bases==col) = lin;
  254.         end
  255.        
  256.          // Preenche o elemento que vai entrar na base
  257.         bases(lin) = col;
  258.                
  259.         // Operações
  260.         for v = 1:(13+n2*10), mprintf('-');end;
  261.             mprintf('\n | Fator Limitante | Base | %c |','z');
  262.             for v = 1:n2-2, mprintf('   x%d   |',v);end;    
  263.             mprintf('   b    |\n |                 |      |');
  264.             for v = 1:(-13+n2*10), mprintf('-');end;
  265.             mprintf('\n |                 |      | 1 |');
  266.             for v = 1:(n2-1), mprintf(' %6.2f |',M(1,1+v));end;  
  267.             mprintf('\n ');
  268.             for v = 1:(13+n2*10), mprintf('-');end;
  269.             mprintf('\n ');
  270.             // ===========================================
  271.             for v = 1:length(q)
  272.                 if(isinf(q(v))) then
  273.                   mprintf('|        Infinito ');
  274.                 else
  275.                     mprintf('| %15.2f ',q(v));
  276.                 end
  277.                 if isnan(bases(v)) then
  278.                     mprintf('|  %2s  | 0 |','-');  
  279.                 else
  280.                     mprintf('|  x%d  | 0 |',bases(v));
  281.                 end
  282.  
  283.               for k = 1:(n2-1)
  284.                 mprintf(' %6.2f |',M(1+v,1+k));
  285.               end
  286.               mprintf(' \n ');
  287.             end
  288.             for v = 1:(13+n2*10), mprintf('-');end;
  289.             mprintf('\n');    
  290.             // ===========================================
  291.             mprintf(' Coluna Pivotal: %d\n',col);
  292.             mprintf(' Linha Pivotal: %d\n',lin);
  293.             mprintf(' Pivô: %.2f\n\n',pivo);
  294.             mprintf(' Operações:\n');
  295.            
  296.             if (isinf(q_min)) then
  297.                 disp(' Não há convergência para estes parâmetros! ');
  298.                 abort;
  299.             end
  300.            
  301.             // Lista as linhas a serem atualizadas
  302.             atualizar = 1:m2;      // Para todas as linhas
  303.             atualizar(lin+1) = []; // Excerto a linha pivotal
  304.    
  305.             // Para cada linha da tabela...
  306.             for i = atualizar
  307.                // Calcula o escalar
  308.                alpha = M(i,col+1)./pivo;
  309.                M(i,:) = M(i,:) - alpha*M(lin+1,:);
  310.                mprintf(' L%d = L%d - %.2f*L%d\n',i-1,i-1,alpha,lin);
  311.             end
  312.            
  313.             // Atualiza a linha pivotal por ultimo
  314.             M(lin+1,:) = M(lin+1,:)./pivo;
  315.             mprintf(' L%d = %.2f*L%d\n\n\n ',lin,1/pivo,lin);
  316. end
  317.  
  318.  
  319. // Operações
  320. for v = 1:(13+n2*10), mprintf('-');end;
  321. mprintf('\n | Fator Limitante | Base | %c |','z');
  322. for v = 1:n2-2, mprintf('   x%d   |',v);end;    
  323. mprintf('   b    |\n |                 |      |');
  324. for v = 1:(-13+n2*10), mprintf('-');end;
  325. mprintf('\n |                 |      | 1 |');
  326. for v = 1:(n2-1), mprintf(' %6.2f |',M(1,1+v));end;  
  327. mprintf('\n ');
  328. for v = 1:(13+n2*10), mprintf('-');end;
  329. mprintf('\n ');
  330. // ===========================================
  331. for v = 1:length(q)
  332.     if(isinf(q(v))) then
  333.         mprintf('|        Infinito ');
  334.     else
  335.         mprintf('| %15.2f ',q(v));
  336.     end
  337.     if isnan(bases(v)) then
  338.         mprintf('|  %2s  | 0 |','-');  
  339.     else
  340.         mprintf('|  x%d  | 0 |',bases(v));
  341.     end
  342.     for k = 1:(n2-1)
  343.         mprintf(' %6.2f |',M(1+v,1+k));
  344.     end
  345.     mprintf(' \n ');
  346. end
  347. for v = 1:(13+n2*10), mprintf('-');end;
  348. mprintf('\n');    
  349. // ===========================================
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top