Znaleźć minimum
\( z =\vec{c}^{T} \cdot \vec{x} \ \ (1) \)
przy warunkach :
\( A\cdot \vec{x} \leq \vec{b}, \ \ (2)\)
\( \vec{x}\geq 0 \ \ (3)\)
gdzie
\( A \) jest daną macierzą o wymiarach \( m\times n, \ \ m\leq n. \)
\( \vec{c}^{T} \) - wektorem kosztów
\( \vec{b} \) - wektorem ograniczeń.
\( \vec{x} \) - wektorem niewiadomych.
Wiele praktycznych zagadnień w ekonomii, technice , finansach w badaniach operacyjnych w teorii podejmowania decyzji może być sformułowanych jako modele LP.
To sprawia, że programowanie liniowe jest popularne i że w ostatnich czterdziestu latach dokonał się fenomenalny jego rozwój.
Do tego przyczyniły się: rozwój algebry opartej na równaniach i nierównościach liniowych i szybki rozwój komputerów jako narzędzia PL.
George Dantzing, Leonid Kantorowicz, Charles Koopmans, Allan Tucker, Cleve Moller - należą do tych wielkich matematyków, którzy rozwinęli teorię, metody obliczeniowe i zastosowania PL.
Na czym polega zrewidowana metoda symplex?
Rozważamy zagadnienie PL w postaci standardowej ze zmiennymi dopełniającymi
\( A\cdot \vec{x} = \vec{b}\)
Zakładając, że macierz \( A \) jest rzędu \( m \) - wydzielamy z niej podmacierz \( B \) stopnia \( m \) i pozostałą podmacierz \( N \) stopnia \( n- m.\)
Zagadnienie PL możemy zapisać w postaci
\( \begin{bmatrix} B & N \end{bmatrix}\cdot \begin{bmatrix} x_{B} \\ x_{N} \end{bmatrix} = \vec{b} \ \ (4)\)
gdzie: \( x_{B}, \ \ x_{N} \) mają odpowiednio \( m \) i \( n-m \) składowych.
Jeśli \( x_{N} \) równa się zero to rozwiązanie \( (4) \) nazywa się rozwiązaniem bazowym \( \vec{x}_{B} = B^{-1}\cdot \vec{b}, \) a nieosobliwą macierz \( B \) nazywa się bazą.
Jeśli ponadto \( x \geq 0 \) to \( \vec{x} \) nazywamy bazowym rozwiązaniem dopuszczalnym.
Następujące twierdzenie pochodzące od twórcy metody simplex George Dantzinga [1] określa ważną rolę bazowych rozwiązań dopuszczalnych.
Twierdzenie
Jeśli istnieje rozwiązanie dopuszczalne spełniające \( (1) - (3) \) to istnieje bazowe rozwiązanie dopuszczalne. Co więcej, jeśli istnieje rozwiązanie optymalne minimalizujące \( z \), to istnieje optymalne rozwiązanie bazowe.
Dowód tego twierdzenia można znaleźć na przykład w [2].
To podstawowe twierdzenie stanowi klucz do algorytmu simplex, najpotężniejszej metody rozwiązywania programów liniowych.
Istnieje kilka różnych wersji tej metody i wiele implementacji numerycznych.
Opiszemy algorytm prymalnej metody symplex w wersji zrewidowanej.
Krok 1
Wybieramy taką bazę, że \( \vec{x}_{B} = = B^{-1}\cdot \vec{b}\geq 0\)
Krok 2
Riozwiązanie równania \( B^{T}\cdot \lambda = c_{B} \) względem wektora mnożników sympleksowych \( \lambda\)
Krok 3
Wybranie z \( N \) takiej kolumny \( a_{k} \), że \( \vec{p}_{k} = = \vec{c}_{N_{k}} -\lambda^{T}\cdot a_{k} < 0.\)
Jeżeli \( p_{T} = \vec{c}_{N}^{T} - \lambda^{T}\cdot N \geq 0, \) to bieżące rozwiązanie jest optymalne.
Krok 4
Rozwiązać układ równań \( B\cdot \vec{y} = a_{k}, \) wzlędem \( \vec{y}.\)
Krok 5
Znaleźć \( \theta = \frac{x_{0l}}{y_{l}} = \min\left [ \frac{x_{0i}}{y_{i}}: y_{i}> 0, \ \ 1\leq i \leq m \right], \)
Jeśli żadne \( y_{i} \) nie jest dodatnie, to zbiór rozwiązań układu \( A\cdot \vec{x} = \vec{b}, \vec{x} \geq 0 \) jest nieograniczony i \( z \) może przyjąć dowolnie dużą wartość ujemną, Zakończyć obliczenia.
Krok 6
Uaktualnić rozwiązanie bazowe, korzystając z równania \( \vec{x}_{B} = \vec{x}_{0} - B^{-1}\cdot \vec{a}_{k}\cdot \vec{x}_{N_{k}} \)
\( \overline{x}_{i} = x_{i} - \theta\cdot , \ \ i\neq k , \)
\( \overline{x}_{k} \) jest nową zmienną bazową.
Krok 7
Uaktualnić bazę \( B \) i wrócić do kroku 2.
Implementacja zrewidowanej metody sympleks w programie MATLAB.
Kod: Zaznacz cały
function[Binv ,B, xB]=rsm(A, c ,B, xB, Binv , v )
% wejście:macierz rozszerzona A, wektor kosztów c ,
%baza B, vector bazowy xB,
% macierz dentycznościowa Binv , znienna niebazowa v
%wyjście: baza B, vektor bazowy xB
[ ~ , n]= size (A) ;
format rat
Cb=c(B);
Cdd=setdiff(1 : n ,B) ;
D=A( : ,Cdd );
Cd=c(1 ,Cdd );
L=Cb*Binv ;
disp('L=');
disp(L );
r=Cd -(L*D);
disp('r= ');
disp(r) ;
[~ , aa]= pivotcolumn(r,v);
% sprawdzenie warunku optymalności
if aa==0
disp('siągnięto rozwiązanie optymalne');
return
else
a=Cdd(aa);
disp('q=');
disp (a) ;
y=Binv*A(: , a);
Y=sprintf('%s ',' y= ');
disp(Y);
disp(y);
revise =[Binv xB y];
[m,n]= size(revise);
min=Inf;
arg_min=0;
end
% wybór elementu głównego
for k=1:m
if revise(k , n)>0
row=revise(k,n-1)/revise(k,n);
if row < min
min=row ;
arg_min=k ;
end
end
end
% rozwiązanie nieograniczone
if arg_min==0
x=sprintf('%c','unbounded');
disp(x);
return
else
disp('p=');
disp(arg_min) ;
% poprawiona tablca simplex
revise=simplex(revise, arg_min, n);
disp('revise');
disp(revise) ;
Binv=revise(:, 1:n-2);
B(arg_min)=a;
xB=revise(:, n-1);
end
end
Kod: Zaznacz cały
> A=[2 -1 2 1 0; 1 0 4 0 1]
A =
2 -1 2 1 0
1 0 4 0 1
>> c= [-6 2 -3 0 0]
c =
-6 2 -3 0 0
>> Binv=eye(2)
Binv =
1 0
0 1
>> v=3
v =
3
>> B=[4 5]
B =
4 5
>> xB=[2;4]
xB =
2
4
>> [Binv,B,xB]=rsm(A,c,B,xB,Binv,v)
L=
0 0
r=
-6 2 -3
q=
1
y=
2
1
p=
1
revise
1/2 0 1 1
-1/2 1 3 0
Binv =
1/2 0
-1/2 1
B =
1 5
xB =
1
3
Kod: Zaznacz cały
>> A=[1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0;
78 71 36 46 45 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
0 0 0 16 98 66 57 54 88 33 41 36 0 0 0 0 0 0 0 1 0 0;
4125 4021 2018 2898 4335 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
0 0 0 638 7548 4638 3812 3905 12497 3990 3841 3277 0 0 0 0 0 0 0 0 0 1]
A =
Columns 1 through 15
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 1 1 1 1 1 1 1 1 1 0 0 0
78 71 36 46 45 0 0 0 0 0 0 0 0 0 0
0 0 0 16 98 66 57 54 88 33 41 36 0 0 0
4125 4021 2018 2898 4335 0 0 0 0 0 0 0 0 0 0
0 0 0 638 7548 4638 3812 3905 12497 3990 3841 3277 0 0 0
Columns 16 through 22
0 0 1 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 1
>> c = [-5334 -4670 -3701 -68305 -16387 -5024 -4193 -5572 -18907 -12165 -7176 -8732 40 40 0.14 0.14 0 0 0 0 0 0]
c =
1.0e+04 *
Columns 1 through 18
-0.5334 -0.4670 -0.3701 -6.8305 -1.6387 -0.5024 -0.4193 -0.5572 -1.8907 -1.2165 -0.7176 -0.8732 0.0040 0.0040 0.0000 0.0000 0 0
Columns 19 through 22
0 0 0 0
>> Binv = eye(6)
Binv =
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
>> v = 16
v =
16
>> xB=[3.44;3.44;661;5800;5800]
xB =
1.0e+03 *
0.0034
0.0034
0.6610
5.8000
5.8000
>> B=[17 18 19 20 21 22]
B =
17 18 19 20 21 22
>> [Binv,b,xB]=rsm(A,c,B,xB,Binv,v)
L=
0 0 0 0 0 0
r=
Columns 1 through 13
-5334 -4670 -3701 -68305 -16387 -5024 -4193 -5572 -18907 -12165 -7176 -8732 40
Columns 14 through 16
40 7/50 7/50
q=
4
y=
1
1
46
16
2898
638
>> xB=[3.44;3.44;661;661;5800;5800]
xB =
86/25
86/25
661
661
5800
5800
>> [Binv,b,xB]=rsm(A,c,B,xB,Binv,v)
L=
0 0 0 0 0 0
r=
Columns 1 through 13
-5334 -4670 -3701 -68305 -16387 -5024 -4193 -5572 -18907 -12165 -7176 -8732 40
Columns 14 through 16
40 7/50 7/50
q=
4
y=
1
1
46
16
2898
638
p=
5
revise
1 0 0 0 -1/2898 0 1125/782 0
0 1 0 0 -1/2898 0 1125/782 0
0 0 1 0 -1/63 0 35843/63 0
0 0 0 1 -8/1449 0 28304/45 0
0 0 0 0 1/2898 0 1451/725 1
0 0 0 0 -319/1449 1 113078/25 0
Binv =
1 0 0 0 -1/2898 0
0 1 0 0 -1/2898 0
0 0 1 0 -1/63 0
0 0 0 1 -8/1449 0
0 0 0 0 1/2898 0
0 0 0 0 -319/1449 1
b =
17 18 19 20 4 22
xB =
1125/782
1125/782
35843/63
28304/45
1451/725
113078/25
[1] George Bernard Dantzing. Linear Programming and Extension. Princeton University Press 1963.
[2] David Luenberger. INTRODUCTION to DYNAMIC SYSTEMS. Theory and models. Stanford University Press 1979.
[3] Maciej M. Sysło Narshing Deo Janusz Kowalik. ALGORYTMY OPTYMALIZACJI DYSKRETNEJ Z PROGRAMAMI W JĘZYKU PASCAL.
Wydawnictwo Naukowe PWN Warszawa 1993,