POLITECHNIKA WROCŁAWSKA	SPRAWOZDANIE Laboratorium Podstaw Optymalizacji
Imię i nazwisko: Maciej Zydlewicz (nr ind: 163051)	Wydział: EKA Kierunek: Automatyka przemysłowa i Robotyka
Temat ćwiczenia: Zrewidowana metoda simpleks	Data wykonania: 6.11.2010

Zadanie 1. Przedsiębiorstwo produkuje dwa wyroby W1 i W2. Ograniczeniem w produkcji są zapasy dwóch surowców A i B. W tablicy 1 podano jednostkowe nakłady surowców na produkcję wyrobów, zapasy surowców oraz ceny wyrobów.

Surowce	Zużycie surowca (w j.p.) na 1 szt. wyrobu	Zapas surowca (w j.p.)
	W1	W2
A	3	2
B	4	5
Cena (w zł)	11	12

Ustalić rozmiary produkcji wyrobów W1 i W2, które gwarantują maksymalny przychód z ich sprzedaży przy istniejących zapasach surowców.

Rozwiązanie. W modelu matematycznym opisującym przedstawioną sytuację decyzyjną występują dwie zmienne decyzyjne: x₁ oznacza wielkość produkcji wyrobu W1, a x₂ to wielkość produkcji wyrobu W2. Model jest następujący:

1. (14 − 3)x₁ + (15 − 3)x₂max,

2. 3x₁ + 2x₂ ≤ 12

3. 4x₁ + 5x₂ ≤ 23

4. x₁, x₂ ≥ 0

Ponieważ w modelu występują tylko dwie zmienne decyzyjne, można go rozwiązać metodą graficzną co czyniliśmy na pierwszym laboratorium. Podany wyżej prosty model, zostanie rozwiązany obecnie za pomocą zrewidowanego algorytmu simpleks. Zaczynamy od sprowadzenia zadania do postaci kanonicznej, dodając do lewych stron nierówności (które są mniejsze) zmienne swobodne: x₃ i x₄. A zatem:

11x₁ + 12x₂ + 0x₃ + 0x₄max,

3x₁ + 2x₂ + x₃                  = 12,

4x₁ + 5x₂ +                  x₄ = 23,

Zauważmy iż zmienną swobodną x₃ można interpretować jako niewykorzystany zasób surowca A, zmienną x₄ - jako nie wykorzystany zasób surowca B.

Tablica simpleksowa dla początkowego rozwiązania bazowego (pierwsza tablica simpleksowa) ma zatem postać tabl. 2 (jak wiemy z teorii, do pierwszej bazy wchodzą zmienne swobodne, natomiast zmienne decyzyjne przyjmują wartość zero).

I tak, zgodnie z oznaczeniami przyjętymi w moim programie: A - jest macierzą współczynników występujących po lewej stronie warunków ograniczających, b – wektorem wyrazów wolnych warunków ograniczających, c – wektorem wierszowym współczynników funkcji celu, x_b – wektorem zmiennych bazowych, a c_b – wektorem kolumnowym współczynników funkcji celu dla zmiennych bazowych. I – jest macierzą jednostkową (macierzą współczynników przy zmiennych występujących w pierwszej bazie). Wartość z_j dla poszczególnych zmiennych (kolumn tablicy) oblicza się jako sumę iloczynów współczynników odpowiadających poszczególnym zmiennym (a_ij) i współczynników funkcji celu dla zmiennych bazowych (c_bi), czyli $z_{j} = \sum_{i = 1}^{m}{a_{\text{ij}}c_{\text{bi}}}$. Ostatni wiersz tablicy simpleksowej c_j − z_j, zwany jest wierszem zerowym lub (kryterium simpleks), jej elementy informują nas o ile zmieni się aktualna w danej iteracji wartość funkcji celu, jeżeli jedną jednostkę tej zmiennej wprowadzimy do nowej (kolejnej) bazy.

I tak, macierzowa postać tablicy simpleksowej dla l-tej iteracji ma postać:

	Zmienne bazowe	c	Rozwiązanie
cbl	xbl	Bl^−1A	Bl^−1	Bl^−1b
	zj	cb^TBl^−1A	cb^TBl^−1	cb^TBl^−1b
	cj − zj	c−cb^TBl^−1A	−cb^TBl^−1

Wracając do naszego zadania, tablica 2 ma postać:

cb	cj	11	12	0	0	Rozwiązanie (bi)
	Zmienne bazowe	x1	x2	x3	x4
0	x3	3	2	1	0	12
0	x4	4	5	0	1	23
	zj	0	0	0	0	0
	cj − zj	11	12	0	0

gdzie,

A = $\begin{bmatrix} 3 & 2 \\ 4 & 5 \\ \end{bmatrix},\ \ I\ = \ \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ \end{bmatrix}\ $, $b = \ \begin{bmatrix} 12 \\ 23 \\ \end{bmatrix},\ \ c = \left\lbrack \begin{matrix} 11 & 12 \\ \end{matrix}\text{\ \ \ \ }\begin{matrix} 0 & 0 \\ \end{matrix} \right\rbrack,\ {\ x}_{b} = \ \begin{bmatrix} x_{3} \\ x_{4} \\ \end{bmatrix},\ \ c_{b} = \ \begin{bmatrix} 0 \\ 0 \\ \end{bmatrix}$

Wartości wiersza z_j dla wszystkich zmiennych w tej tablicy są równe zeru, bo współczynniki przy wszystkich zmiennych bazowych są równe zeru. Np. dla x₁: z₁=3*0+4*0=0.

Ponieważ funkcja celu jest maksymalizowana, do kolejnego rozwiązania bazowego wejdzie – zmienna o największej wartości kryterium simpleks. (czyli największej wartości w wierszu zerowym c_j − z_j). Jest to zmienna x₂ (max(c_j − z_j)=c₂ − z₂=12). Należy teraz ustalić w miejsce której z dotychczasowych zmiennych zostanie wprowadzona. W tym celu wartości zmiennych (b_i) należy podzielić przez współczynniki stojące przy wprowadzanej zmiennej w aktualnej tablicy simpleksowej (tablica 2) i wybrać zmienną, dla której ten iloraz jest najmniejszy. (Uwaga! Jeżeli w kolumnie zmiennej wprowadzanej do bazy występują współczynniki ujemne lub równe zeru, pomijamy te zmienne w dalszych obliczeniach). Spośród dwóch ilorazów: 12/2=6, 23/5=4.6. Najmniejszy odpowiada zmiennej x₄. Co można wytłumaczyć tym, że zapas surowca B ogranicza produkcję wyrobu W2 do 4.6 szt. Zatem w kolejnej tablicy simpleksowej (tablica 3) zmiennymi bazowymi są x₃ i x₂. A poszczególne elementy tej tablicy można obliczyć, stosując wzory macierzowe podane w tablicy simpleksowej. Wobec tego tablica 3 ma postać:

A macierz B₁ = $\begin{bmatrix} 1 & 2 \\ 0 & 5 \\ \end{bmatrix}$

cb	cj	11	12	0	0	Rozwiązanie (bi)
	Zmienne bazowe	x1	x2	x3	x4
0	x3	1.4	0	1	-0.4	2.8
12	x2	0.8	1	0	0.2	4.6
	zj	9.6	12	0	2.4	55.2
	cj − zj	1.4	0	0	-2.4

Tym razem x₁ wchodzi do bazy i podmienia zmienną x₃.

Kolejna iteracja doprowadzi nas do rozwiązania opisanego tablicą 4:

cb	cj	11	12	0	0	Rozwiązanie (bi)
	Zmienne bazowe	x1	x2	x3	x4
11	x1	1	0	0.7143	-0.2857	2
12	x2	0	1	-0.5714	0.4286	3
	zj	11	12	1	2	58
	cj − zj	0	0	-1	-2

przy macierzy B₂ = $\begin{bmatrix} 3 & 2 \\ 4 & 5 \\ \end{bmatrix}$

Jak nie trudno zauważyć w wierszu zerowym tablicy 4 nie występują już liczby dodatnie, zatem rozwiązanie jest optymalne.

Zadanie 2. Określenie zagadnienia oraz równań PL sprowadzonych do postaci kanonicznej odpowiednio dla metody simpleks.

1. x₁₁ + 2x₁₂ + 3x₁₃ + 3x₂₁ + 2x₂₂ + 2x₂₃max,

2. x₁₁ + x₂₁ = 2,

3. x₁₂ + x₂₂ = 3,

4. x₁₃ + x₂₃ = 4,

5. x₁₂ + x₁₂ + x₁₃ +  z₁ ≤ 5,

6. x₂₁ + x₂₂ + x₂₃ +  z₂ ≤ 6,

7. x₁₁, x₁₂, x₁₃, x₂₁, x₂₂, x₂₃ ≥ 0.

Postać kanoniczna:

1. x₁₁ + 2x₁₂ + 3x₁₃ + 3x₂₁ + 2x₂₂ + 2x₂₃max,

2. 1x₁₁ + 0x₁₂ + 0x₁₃ + 1x₂₁ + 0x₂₂ + 0x₂₃ = 2,

3. 0x₁₁ + 1x₁₂ + 0x₁₃ + 0x₂₁ + 1x₂₂ + 0x₂₃ = 3,

4. 0x₁₁ + 0x₁₂ + 1x₁₃ + 0x₂₁ + 0x₂₂ + 1x₂₃ = 4,

5. 1x₁₁ + 1x₁₂ + 1x₁₃ + 0x₂₁ + 0x₂₂ + 0x₂₃ = 5,

6. 0x₁₁ + 0x₁₂ + 0x₁₃ + 1x₂₁ + 1x₂₂ + 1x₂₃ = 6.

Implementacja algorytmu w matlabie:

function solution = rsm (c, A, b)

%

% Zadanie: fc-->max przy Ax <= b & x >= 0

%

% m liczba wierszy w macierzy A

% n liczba kolumn w macierzy A

%

[m n] = size(A);

%

% Określenie warunków początkowych

%

% I macierz jednostkowa o wymiarach [mxm]

% cb wektor kolumnowy współczynnikó funkcji celu dla zmiennych

% bazowych o wymiarach [mx1]

% cc wektor wierszowych wspołczynników funkcji celu

% (jak w tablicy simpleks)

% X wektor z numerami zmiennych bazowych

X=[n+1:m+n]';

B = [A,eye(m)];

D=B;

cb=zeros(m,1);

cc=[c, cb'];

%

% iters numer l-tej iteracji

%

iters=0;

while 1==1

iters=iters+1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 1

% wyznaczanie macierzy odwrotnej (B^-1)

%

% Binv odwrócenie macierzy bazowej

%

Binv = inv(B(:,n+1:m+n));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 2

% wyznaczanie bierzącego wektora rozwiązań (wartości zmiennych bazowych)

%

% xb bieżący wektor rozwiązań

%

xb = Binv * b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 3

% wyznaczanie wartosci bieżącej funkcji celu fc = cb * B^-1 * b

%

% fc modyfikowany wektor kosztów

%

fc = cb'*Binv*b

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 4

% wyszukiwanie największej wartosci kryterium simpleks

[cs j]=max(cc-[(cb'*Binv*A),cb'*Binv]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 5

% if cs <= 0 , wówczas znaleziono rozwiązanie optymalne

%

if (cs <= 0.0005) %uwzględniam dokładność przy odejmowani (dla funkcji max)

rozwiazanie = zeros(n,1);

rozwiazanie = xb

return;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 6

% rozwiąż w = B^-1 * a[j]

%

% w względna waga (wektor) kolumny wchodzącej do bazy

w = Binv*D(:,j);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 7

% rozwiązanie w[i]>0 i wyszukanie najmniejszego dodatniego ilorazu b[i]/w[i]

% wymiana kolumn "j" z bazą i wymiana kolumn "i" wychodzących z bazy

% mn minimalna wartość ilorazu b[i]/w[i] gdzie w[i] > 0

% i wiersze korespondujące z mn -- określające wychodzącą kolumne

% C tymczasowy wektor

mn = inf;

i=0;

zz = find (w > 0)';

[yy, ii] = min(xb(zz) ./ w(zz));

i = zz(ii(1));

if (i == 0)

error ('Rozwiązanie jest nieokreślone.')

end

cb(i,1) = cc(1,j);

B(:,X(i))=D(:,j);

X(i)=j;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Krok 8

% Powtórz procedurę

end; % while

WYWOŁANIA FUNKCJI Z LISTY 1:

% c=[11 12]; A=[3 2;4 5]; b=[12 23];rsm(c,A,b') %wywołanie funkcji z zad 1

% c=[1 2 3 3 2 2]; A=[1 0 0 1 0 0;0 1 0 0 1 0;0 0 1 0 0 1;1 1 1 0 0 0;0 0 0 1 1 1]; b=[2 3 4 5 6]; rsm(c,A,b') %wywołanie funkcji z zad 2

ROZWIĄZANIA ZADAŃ:

Zad.1(lista 1):

>> c=[11 12]; A=[3 2;4 5]; b=[12 23];rsm(c,A,b')

iters = 1

fc = 0

iters = 2

fc = 55.2000

iters = 3

fc = 58.0000

rozwiazanie =

2.0000

3.0000

Zad.2 (lista 1):

>> c=[1 2 3 3 2 2]; A=[1 0 0 1 0 0;0 1 0 0 1 0;0 0 1 0 0 1;1 1 1 0 0 0;0 0 0 1 1 1]; b=[2 3 4 5 6]; rsm(c,A,b')

iters = 1

fc = 0

iters = 2

fc = 12

iters = 3

fc = 18

iters = 4

fc = 20

iters = 5

fc = 24

rozwiazanie =

2

2

4

1

2

Analiza wrażliwości. Rozwiązanie optymalne programu liniowego umożliwia podjęcie trafnej decyzji. Równocześnie jednak może stanowić punkt wyjścia do analizy, jak i zmiany warunków działania przedsiębiorstwa, mające odzwierciedlenie w zmianach parametrów modelu, mogą wpłynąć na rozwiązanie optymalne. Badanie wpływu zmian wartości parametrów na rozwiązanie optymalne programu liniowego nosi nazwę analizy wrażliwości. Jakkolwiek programy liniowe należą do modeli deterministycznych, a więc z założenia parametry są znane i stałe, jednakże w dłuższych okresach czasu parametry mogą ulegać zmianom. Zmieniać się mogą w szczególności ceny (produktów sprzedawanych, lub kosztów produkcji), ale także zasobów środków produkcji, zdolności produkcyjne, itd. Przedmiotem analizy wrażliwości są więc:

1. Współczynników funkcji celu – analiza wrażliwości pozwala na określenie w jakich granicach mogą się zmieniać te parametry, aby dotychczasowe rozwiązanie pozostało optymalne,

2. Wyrazów wolnych w warunkach ograniczających – co pozwala określić w jakich przedziałach liczbowych mogą się zmieniać wyrazy wolne (prawostronne ograniczenia), aby w rozwiązaniu optymalnym pozostały dotychczasowe zmienne bazowe, oraz wyznaczyć nowe optymalne wartości tych zmiennych,

3. Współczynników występujących po lewej stronie układu warunków ograniczających, a więc pewnych norm technicznych,

4. Dodanie nowych warunków ograniczających.

W praktyce najczęściej ogranicza się do badania p.pkt. a) i b).

Analizie poddamy tylko przykład z zadania 1. Tym razem jednak wykorzystamy uzyskane na końcu rozwiązanie. Oceny wrażliwości dokonamy tylko w oparciu o p.pkt. a). zastępując dotychczasowe ceny c_j wyrażeniami c_j^′ = c_j + δ_j, gdzie δ_j jest dopuszczalną zmianą ceny c_j, przy jakiej dotychczasowe rozwiązanie optymalne nie ulegnie zmianie. Wielkości δ_j określa się pamiętając, że w zagadnieniach maksymalizacji wszystkie elementy wiersza zerowego końcowej tablicy simpleksowej muszą być niedodatnie, natomiast w zagadnieniach minimalizacji elementy te powinny być nieujemne.

Zatem dla wyrobu W1 przyjmujemy, że c₁^′ = c₁ + δ₁. Wstawiając tę wartość w miejsce c₁ do ostatniej tablicy simpleksowej, otrzymujemy:

cb	cj	11+δ1	12	0	0	Rozwiązanie (bi)
	Zmienne bazowe	x1	x2	x3	x4
11+δ1	x1	1	0	0.7143	-0.2857	2
12	x2	0	1	-0.5714	0.4286	3
	zj	11+δ1	12	1+0.7143δ1	2−0.2857δ1	58
	cj − zj	0	0	-1−0.7143δ1	-2+0.2857δ1

Aby znaleźć δ₁ należy rozwiązać układ dwóch nierówności liniowych (tylko te dwa elementy wiersza zerowego zależą od δ₁):

1. -1−0.7143δ₁ ≤ 0

2. -2+0.2857δ₁ ≤ 0

Pierwsza nierówność jest spełniona dla δ₁≥ $- \frac{10000}{7143}$, druga dla δ₁≤ $\frac{20000}{2857}$. Zatem $\delta_{1} \in < - \frac{10000}{7143},\frac{20000}{2857} >$, a $c_{1} \in < 11 - \frac{10000}{7143},11 + \frac{20000}{2857} >$. Czyli obecne rozwiązanie optymalne nie ulegnie zmianie, jeżeli cena wyrobu W1 będzie przyjmować wartości z przedziału <9.6 , 18>. Należy jednak pamiętać, że nie zmienią się oczywiście tylko optymalne wielkości x₁^* i x₂^*, natomiast zmieni się (wzrośnie lub spadnie) wartość funkcji celu (o δ₁x₁^* ). Analogicznie liczy się dopuszczalne granice zmian ceny wyrobu W2, które wynoszą < 7.33, 13.75>.

Wnioski. Metoda simpleks jest uniwersalną metodą pozwalającą na rozwiązywanie zagadnień programowania liniowego. Istnieją co prawda dwie metody rozwiązywania programów liniowych: metoda graficzna (geometryczna, z którą mieliśmy okazję zaznajomić się w sposób praktyczny na pierwszym laboratorium) oraz metoda Simpleks. Jednakże rozwiązywanie układów programowania linowego w praktycznych problemach składających się z większej liczby zmiennych decyzyjnych i ilości ograniczeń rzędu kilkudziesięciu (co miało miejsce w zadaniu 2) metodą graficzną staje się bardzo utrudnione. Rozwiązań programu liniowego należy bowiem poszukiwać wśród wierzchołków wielościanu określonego przez ograniczenia i warunki brzegowe. Istota metody Simpleks polega na badaniu rozwiązań bazowych programu kanonicznego w taki sposób, że znajdujemy wyjściowe rozwiązanie bazowe programu liniowego. Mając rozwiązanie bazowe sprawdzamy czy jest ono optymalne czy nie. Jeśli dane rozwiązanie bazowe nie jest optymalne budujemy następne rozwiązanie bazowe lepsze lub przynajmniej nie gorsze od poprzedniego, z którym to rozwiązaniem postępujemy tak samo jak z rozwiązaniem wyjściowym. Aż do momentu znalezienia rozwiązania optymalnego. Ponadto co przedstawiłem w swoim sprawozdaniu, metoda simpleks umożliwia łatwe przeprowadzenie analizy wrażliwości zagadnień programowania liniowego, co może być pomocne przy podejmowaniu decyzji.