Simplex Analyse

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ausarbeitung zur Simplex Analyse in Matlab 7.1 mit ausführlichen Kommentar
% Vorlesung: Optimierung
% Copyright 2007 Thomas Ludwig
% This file is distributed under the terms of the GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameter c ist für den Beispiel Switcher in der Funktion
% simplex_beispiele(c) - c= Integer (1 - 5)
function simplex_analyse(c);
 
    % Funktion mit mehreren Beispielen aufrufen
    [f,f_wert,A,b,Aeq,beq,lb,ub,x0,options]= simplex_beispiele(c);

    % Linprog Parameter Erklaerung
    % [x,fval,exitflag,output,lambda]=
    % Ausgangs - Argumente:
    %        x: Koeffizienten Vektor des optimalen Punktes
    %        fval: Funktionswert am optimalen Punkt
    %        exitflag: Integer Wert gibt den Grund an, weshalb der Algorithmus
    %                  abgebrochen wurde
    %                   1 : Erfolgreich. Bei der Lösung x
    %                   0 : due to maximum iteration
    %                  -2 : kein durchfuehrbarer Punkt
    %                  -3 : Problem ist unbegrenzt (kein Convexes Problem)
    %        output: Information ueber die Optimierung (Iteration, Algorithmus)
    %        lambda: Lagrange Vervielfacher bei der x'ten Loesung
    %
    % linprog(-f,A,b,Aeq,beq,lb,ub,x0,options);
    % Eingangs - Argumente:
    %        f: Koeffizientenvektor der Zielfunktion
    %        A: Koeffizientmatrix der Nebenbedingungen
    %        b: Restriktionen der Nebenbedingungen
    %        Aeq: Koeffizientenmatrix der Gleichheitsbegrenzungen
    %        beq: Restriktionenvektor der Gleichheitsbegrenzungen
    %        lb: Koeffizientenvektor der untereren Grenze (Lower Bound)
    %        ub: Koeffizientenvektor der oberen Grenze (Upper Bound)
    %        x0: Koeffizientenvektor der Ausgangsvermutung
    %        options: Definiert die Optimierungsoptionen
    %                 algorithm, tolerance, Jacobian
    %                 display, die Anzahl der Wiederholungen, usw.

    % Options Erklaerung - optimset(Parameter, value)
    % Largescale: off  -> Die Simplex Methode unterstützt keine Largescale Probleme
    % Simplex: on      -> Benutze das Simplex Verfahren
    % Display: iter    -> Jeden Durchgang anzeigen (iteration = Wiederholung)

    % Optionen pruefen die bereits in der Funktion simplex_beispiele
    % aufgerufen wurden
    if isempty(options);
        options = optimset('Largescale','off', 'Simplex','on', 'Display','iter');
    end
   
    [x,fval,exitflag,output,lambda]= linprog(-f,A,b,Aeq,beq,lb,ub,x0,options);
    fval= -(fval - f_wert); % Vorherigen Zielfunktionswert mit einbeziehen
   
    % Fehlermeldungen Auswerten
    switch exitflag
        case 0
          disp('Fehler: Abgebrochen wegen der maximalen Wiederholung');
        case -2
          disp('Fehler: Kein durchführbarer Punkt');
        case -3
          disp('Fehler: Das Problem hat keine Grenzen (Kein Konvexes Problem)');

        otherwise
            % Ausgabe
            output
            lambda 
            fprintf(1,'f  = %g\n', fval);
            for (i=1:length(x));
                fprintf(1,'x%g = %g \n', i, x(i));
            end       

            % Grafik Funktion laden. Diese entscheidet selbst über 2 oder
            % mehrdimensionale Darstellung der Optimalen Lösung
            grafik(f,A,b,x,lb,ub,fval,f_wert);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotten
function grafik(f,A,b,x,lb,ub,fval,f_wert);

     % Grenzen setzen
     if isempty(ub);
        ub= max(b);
     end
     if isempty(lb);
        lb= min(b);
     end

     x_range= min(lb):max(ub); % Skalenlänge der x Achse
     [kx, ky]= size(x_range);    
    
    %
    % Ordnung der Zielzeile
    %
    if (length(f) == 2)
       
        subplot(2,1,1);

            % Gleichungen berechnen
            for (i=1:length(b));
                y(i,kx:ky)= max( ( b(i)- ( A(i,1)*x_range ) ) / A(i,2), 0 );
            end

            area(x_range,min([y])), hold on;
            colormap summer;
            [u, v]= meshgrid(x_range, y);

            % Funktion erstellen
            z= f(1)*u+f(2)*v+f_wert;

            contour(u,v+v,z);

            % Optimalen Lösungspunkt Plotten
            for i=1:x(1);
                for j=1:x(2);
                    plot(i,x(2),'-r',x(1),j,'-r');
                end
            end

            % Zeichnet alle Funktionen inkl. des Optimalen Punktes ein
            plot(x_range,y,'-r','LineWidth',2);
            plot(x(1),x(2),'--bs');

            % Immer lesbare Beschriftung der Optimalen Lösung
            if ( (length(x_range) - x(1)) <= (length(x_range) / 3) )
                text(x(1),x(2),'Optimale Lösung \rightarrow    ', 'HorizontalAlignment','right');
            else
                text(x(1),x(2),'    \leftarrow Optimale Lösung', 'HorizontalAlignment','left');
            end

            % Beschriftung der Grafik
            title('Geometrische Interpretation für die Optimale Lösung');
            xlabel('x'), ylabel('y'), hold off;
 
        subplot(2,1,2);
            mesh(u,v,z), xlabel('x1'), ylabel('x2'), hold on;
           
            % Optimierten Bereich kennzeichnen
            [a, b]= meshgrid(0:x(1), 0:x(2));
            z3= f(1)*a + f(2)*b + f_wert;
            meshc(a, b, z3);
           
            [a1, b1]= meshgrid(0:x(1), x(2));
            z4= f(1)*a1 + f(2)*b1 + f_wert;
            plot3(a1, b1, z4,'-r', 'LineWidth',3);

            [a2, b2]= meshgrid(x(1), 0:x(2));
            z5= f(1)*a2 + f(2)*b2 + f_wert;
            plot3(a2, b2, z5,'-r', 'LineWidth',3);
           
            plot3(x(1),x(2),fval,'--bs'), hidden off; % Optimale Loesung markieren

            plot3(x(1),x(2),0:fval/100:fval, '-r', 'LineWidth',5); % Linie Z-Achse
            text(x(1),x(2),fval,'\leftarrow Optimale Lösung', 'HorizontalAlignment','right', 'VerticalAlignment','bottom' )
            hold off;
      
    %
    % Mehrdimensionale Darstellung
    %
    elseif (length(f) == 3)
 
        [u, v]= meshgrid(x_range, x_range);
      
        % 1. Grafik initialisieren
        z= f(1)*u + f(2)*v + f(3)*0.00001 + f_wert;
        mesh(u,v,z), hold on;

        for k = 0:(x(3)/50):x(3)
            z= f(1)*u + f(2)*v + f(3)*k + f_wert;
            mesh(u,v,z), xlabel('x1'), ylabel('x2'), zlabel('x3');
            title('Geometrische Interpretation für die Optimale Lösung');
           
            drawnow;
            pause(.1);
        end

        hold off;
    end
   
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Beispiele
function [f,f_wert,A,b,Aeq,beq,lb,ub,x0,options]= simplex_beispiele(c);
   
    % Beispiel Aufgaben Switcher
    switch c
       
        case 1
            % Aufgabe
            % -> 2 Produkte (x1, x2)
            % -> 3 Maschienen
            % -> Maximiere f
            %
            % Zielfunktion
            % f= 300x1 + 500x2 - 36000
            %
            % Nebenbedingungen
            % 1*x1 + 2*x2 <= 170 (Maschine A)
            % 1*x1 + 1*x2 <= 150 (Maschine B)
            % 0*x1 + 3*x2 <= 180 (Maschine C)
            %
            % Grenzwerte
            % x1, x2 >= 0

            f= [300; 500];
            f_wert= -36000;
            A= [1 2;
                1 1;
                0 3];
            b= [170; 150; 180];
            Aeq= [];
            beq= [];           
            lb= [0];
            ub= [];
            x0= [];
            options= [];
           
           
        case 2
            % Aufgabe
            % -> 3 Produkte (x1, x2, x3)
            % -> 2 Maschienen
            % -> Maximiere f
            %
            % Zielfunktion
            % f= 1*x1 + 2*x2 + 3*x3
            %
            % Nebenbedingungen
            % -1*x1 +   x2 + x3 <= 20 (Maschine A)
            %    x1 - 3*x2 + x3 <= 30 (Maschine B)
            %
            % Grenzwerte
            % 0 <= x1 <= 40
            % x2,x3 >= 0
        
            f= [1; 2; 3];
            f_wert= 0;
            A= [-1  1 1;
                 1 -3 1];
            b= [20; 30];
            lb= [0];
            ub= [40];
            Aeq= [];
            beq= [];
            x0=  [];
            options= [];
           
           
        case 3
            % Aufgabe
            % -> 3 Produkte (x1, x2, x3)
            % -> 3 Maschienen
            % -> Sonderfall = in Nebenbedingung
            % -> Maximiere f
            %
            % Zielfunktion
            % f= 4*x1 + 2*x2 + 1*x3
            %
            % Nebenbedingungen
            % 2*x1 + 1*x2 + 0*x3 <= 1
            % 1*x1 + 0*x2 + 2*x3 <= 2
            % 1*x1 + 1*x2 + 1*x3  = 1
            %
            % Grenzwerte
            % 0 <= x1 <= 1
            % 0 <= x2 <= 1
            % 0 <= x3 <= 2
         
            f= [4; 2; 1];
            f_wert= 0;
            A= [2 1 0;
                1 0 2];
            b= [1; 2];
            lb= [0; 0; 0];
            ub= [1; 1; 2];
            Aeq= [1 1 1];
            beq= [1];
            x0=  [];
            options= [];  

        case 4
            % Aufgabe
            % -> 2 Produkte (x1, x2)
            % -> 3 Maschienen
            % -> Maximiere f
            %
            % Zielfunktion
            % f= 143x1 + 60x2
            %
            % Nebenbedingungen
            % 120*x1 + 210*x2 <= 15000 (Maschine A)
            % 110*x1 +  30*x2 <=  4000 (Maschine B)
            %   1*x1 +   1*x2 <=    75 (Maschine C)
            %
            % Grenzwerte
            % x1, x2 >= 0
           
            f = [143 60];
            f_wert= 0;
            A = [120 210;
                 110  30;
                   1   1];
            b = [15000; 4000; 75];        
            Aeq= [];
            beq= [];           
            lb= [0];
            ub= [100; 80;];
            x0=  [];
            options= [];
           
        case 5
            % Aufgabe Script Optimieren (S.24)
            % -> x1= Zangen, x2= Schraubenschlüssel
            % -> 3 Maschinen
            % -> Maximiere f
            %
            % Zielfunktion
            % f= 300x1 + 400x2
            %
            % Nebenbedingungen
            % 1  *x1 +  1,5*x2 <= 15
            % 1  *x1 +  1  *x2 <= 12
            % 0,5*x1 +  0,4*x2 <=  5
            % 1  *x1 +  0  *x2 <= 10
            % 0  *x1 +  1  *x2 <=  8
            %
            % Grenzwerte
            % x1, x2 >= 0           
           
            f = [300 400];
            f_wert= 0;
            A = [1 1.5; 1 1; 0.5 0.4];
            b = [15; 12; 5;];        
            Aeq= [];
            beq= [];           
            lb= [ 0; 0];
            ub= [10; 8];
            x0=  [];
            options= [];           
           
    end

end