primalsimplex_phase_2.m

primalsimplex_phase_2.m — Objective-C source code, 3 kB (3.738 bytes)

Dateiinhalt

%----------------------------------------------------------------------
%
% Phase II (Matrixversion) des primalen Simplexverfahrens zur Loesung 
% des linearen Optimierungsproblems 
% 
%     Minimiere c'x u.d.N. Ax = b, x>=0
%
% Die Pivotelemente werden nach der Bland'schen Regel oder 
% nach der Regel von Dantzig gewaehlt. 
%
% Eingabe:
% ========
% A     : Matrix A 
% b     : Vektor b
% c     : Vektor c
% J     : Basisindexmenge einer zulaessigen Startbasisloesung
% pivot : 0 (=Dantzig) oder <>0(=Bland'sche Regel)
%
% Ausgabe:
% ========
% x     : Loesungsvektor
% y     : dualer Loesungsvektor
% f     : Zielfunktionswert
% Jopt  : Basisindexmenge fuer x
%
%----------------------------------------------------------------------
   function [x,y,f,Jopt] = primalsimplex_phase_2(A,b,c,J,pivot)
   [m,n] = size(A);
   if (m ~= size(b,1)) error('Falsche Dimension von A und b!'); end;
   if (n ~= size(c,1)) error('Falsche Dimension von A und c!'); end;
   if (m ~= size(J,2)) error('Basisindexmenge hat falsche Dimension!'); end;
   x = zeros(n,1);
   zaehler = 0;
   %
   % Toleranz fuer Abfragen
   %
   tol = 1.0e-13;
   %
   % Bestimme Nichtbasisindexmenge Jc
   % 
   Jc = [];
   for i1=1:1:n
      flag = 0;
      for i2=1:1:m
         if (J(i2) == i1) flag=1; end;
      end;
      if (flag == 0) Jc = [Jc,i1]; end;
   end;
   % 
   % Schleife
   % 
   while (1 > 0)
     %
     % Berechne Hilfsgroessen
     % 
     XJJc   = A(:,J)\A(:,Jc);
     pJ     = A(:,J)\b;
     zetaJc = XJJc'*c(J) - c(Jc);
     x(J)   = pJ;
     x(Jc)  = 0;
     f      = c(J)'*pJ;
     Jopt   = J;
     y      = A(:,J)'\c(J);
     %
     % Tableauausgabe 
     %
     fprintf(' Tableau %d \n',zaehler);
     fprintf('      |'); 
     for i=1:1:n-m
	     fprintf('%14d',Jc(i));
     end;
     fprintf(' |\n');
     for i=1:1:m
	     fprintf('%5d |',J(i));
	     for i2=1:1:n-m
            fprintf('%14.5e',XJJc(i,i2));
	     end;
         fprintf(' | %14.5e \n',pJ(i));
     end;
     fprintf('      |'); 
     for i=1:1:n-m
        fprintf('%14.5e',zetaJc(i));
     end;
     fprintf(' | %14.5e \n',f);
     %
     % Optimalitaetskriterium
     % 
     if all(zetaJc <= tol) 
        display('Tableau ist optimal!')
        return
     end;
     %
     % Waehle Pivotspalte (Bland'sche Regel oder Dantzig)
     % 
     if (pivot == 0) 
        %
        % Dantzig
        %
        [zetak,ik] = max(zetaJc);
        k = Jc(ik);
     else
        % 
        % Bland
        %
        k = n+1;
        for i=1:1:n-m
          if (zetaJc(i) >= tol) 
             if (Jc(i) < k) 
                k     = Jc(i);
                zetak = zetaJc(i);
                ik    = i;
             end;
          end;
        end;
     end;
     if (k < 1) || (k > n), 
        disp('Fehler bei Bestimmung der Pivotspalte!');
        return;
     end;
     %
     % Unbeschraenkheitskriterium
     %
     if all(XJJc(:,ik) <=0) 
        display('Problem ist unbeschraenkt!')
        return;
     end;
     %
     % Wahl der Pivotzeile (Bland'sche Regel)
     %
     l    = 0;
     qmin = 1.0e21;
     for i=1:1:m
        if (XJJc(i,ik) > tol) 
           tmp = pJ(i)/XJJc(i,ik);
           if (tmp <= qmin) 
              if (tmp == qmin && J(i) < l) 
                 l  = J(i);
                 il = i;
              elseif (tmp < qmin) 
                 qmin = tmp;
                 l    = J(i);
                 il   = i;
              end;
           end;
        end;
     end;
     if (l < 1) || (l > n), 
        disp('Fehler bei Bestimmung der Pivotzeile!');
        return;
     end;
     fprintf(' Pivotzeile  = %d \n Pivotspalte = %d \n',l,k);
     %
     % Basiswechsel
     %
     J(il)  = k;
     Jc(ik) = l;
     zaehler = zaehler + 1;
   end;