Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 1
Poszukiwanie minimów funkcji wielu zmiennych
Określenie zagadnienia
Problem polega na poszukiwaniu minimum funkcji wielu zmiennych. Funkcja może być dowolna (tzn. dowolnego rzędu,
także nieliniowa). Do badania w przykładach wybierzemy następującą funkcję (jest to obwód jakiejś figury geometrycznej
uzależniony od długości dwóch jej boków przy ustalonej wartości powierzchni całej figury, oznaczymy ją przez P(x,y)):
2
4 2 2 4
1 x 2 x y y 4
2
P(x,y) x y x y x y
x y x y
Funkcja badana faktycznie ma minimum wynoszące 2.8284271507 w punkcie x=y=0.7070110163 (są to wartości
obliczone numerycznie). Podczas ilustrowania działania omawianych dalej metod szukania minimum będę zakładać jako
punkt wyjściowy punkt o współrzędnych (2.25,3.00).
Wykres badanej funkcji
10
8
6
4
0.5
1
1.5
2
0.5
2.5
1
1.5
y
2
2.5
3
3
x
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 2
(czerwona gwiazdka na rysunku oznacza faktyczne położenie minimum lokalnego).
Znamy analityczne wyrażenia na składowe jej gradientu:
2 4 2 2 4 4 3 3 4
P(x,y) x y x 2 x y y 4 x 2 x y 2 x y y 4
x 2 4 2 2 4
x y x 2 x y y 4
2 4 2 2 4 4 3 3 4
P(x,y) x y x 2 x y y 4 x 2 x y 2 x y y 4
y 2 4 2 2 4
x y x 2 x y y 4
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 3
Jest to dość  perfidna funkcja, widać to dopiero po obejrzeniu jej wykresu dla większego zakresu argumentów  w
zasadzie z faktu, że ma to być obwód figury w funkcji długości boków wynika, że x i y musza być dodatnie, jednak podczas
obliczeń numerycznych wartości pośrednie zmiennych mogą być bardzo różne.
Obok przedstawiam wykres tej właśnie funkcji dla większego zakresu zmienności argumentów  widać na nim, że
funkcja ma zarówno minimum lokalne jak i maksimum lokalne, ale globalnie minimum dąży do  , a maksimum do + .
Metoda połowienia kroku
Metoda jest najprostszą adaptacją metod poznanych dla funkcji jednej zmiennej. Przebieg obliczeń jest następujący:
1. Wybieramy punkt startowy (xo,yo), obliczamy wartość P(xo,yo) i gradient grad(P(x,y)) badanej funkcji w tym punkcie,
2. W kierunku przeciwnym do obliczonego gradientu  wykonujemy krok o wybranej wcześniej długości  obliczamy
nowe współrzędne punktu (x1,y1) i nową wartość funkcji P(x1,y1),
3. Jeśli nowa wartość funkcji jest mniejsza niż wartość z poprzedniego punktu obliczamy jej gradient w tym punkcie
grad(P(x1,y1)) i wykonujemy następny krok (czyli powtarzamy etap 2 powyżej), jeśli jednak wartość funkcji nie zmalała
to wracamy do poprzedniego punktu, zmniejszamy krok (np. o połowę), i powracamy do etapu 2  wykonujemy ten
skrócony krok i ponownie sprawdzamy wartość funkcji.
Już na podstawie tego opisu możemy zauważyć, że jeśli trafimy wreszcie (być może przypadkiem) do minimum to
kolejny krok przeniesie nas dalej  do punktu o większej wartości funkcji i możemy powtarzać etap 2 w nieskończoność
(albo do chwili gdy osiągniemy zadaną minimalną długość kroku).
Oto przykładowy przebieg obliczeń w metodzie  połowienia kroku :
cykl 1, wartość funkcji: 6.0912 błąd wzgl.: 1.15357
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 4
cykl 2, wartość funkcji: 4.72901 błąd wzgl.: 0.671957
cykl 3, wartość funkcji: 3.69826 błąd wzgl.: 0.307534
w cyklu 4 zmniejszam krok z 1 na 0.5
cykl 4, wartość funkcji: 2.93152 błąd wzgl.: 0.036449
w cyklu 5 zmniejszam krok z 0.5 na 0.25
cykl 5, wartość funkcji: 2.8836 błąd wzgl.: 0.0195064
w cyklu 6 zmniejszam krok z 0.25 na 0.125
cykl 6, wartość funkcji: 2.83676 błąd wzgl.: 0.00294435
w cyklu 7 zmniejszam krok z 0.125 na 0.0625
cykl 7, wartość funkcji: 2.83149 błąd wzgl.: 0.00108422
w cyklu 8 zmniejszam krok z 0.0625 na 0.03125
cykl 8, wartość funkcji: 2.82883 błąd wzgl.: 0.000142345
w cyklu 9 zmniejszam krok z 0.03125 na 0.015625
cykl 9, wartość funkcji: 2.82872 błąd wzgl.: 0.000102758
w cyklu 10 zmniejszam krok z 0.015625 na 0.0078125
w cyklu 10 zmniejszam krok z 0.0078125 na 0.00390625
w cyklu 10 zmniejszam krok z 0.00390625 na 0.00195313
cykl 10, wartość funkcji: 2.82843 błąd wzgl.: 9.56914e  007
Wyniki metody  połowienia kroku :
współrzędne minimum: x_min=0.708081, y_min=0.708099, błędy wzgl. x i y: 0.00151403, 0.00153917
wartość funkcji: 2.82843 błąd wzgl.wartości: 9.56914e-007, wartość prawdziwa: 2.82843
wartości składowych gradientu: x-owa: 0.00276355, y-owa: 0.00278872
Ogólnie taka metoda postępowania działa dość dobrze, jednak jeżeli mamy pecha i wybierzemy zbyt dużą wartość
pierwszego kroku  możemy  przeskoczyć szukane minimum i potem już tylko oddalać się od niego do nieskończoności
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 5
(wybrana funkcja przykładowa ma właśnie taką własność  istnieje dobrze określone minimum lokalne , ale w połowie
swej dziedziny funkcja maleje do  ).
Drugą wadą jest stosunkowo powolna zbieżność  jeśli wybraliśmy zbyt krótki krok początkowy poszukiwania
minimum mogą się bardzo dłużyć. Rysunek powyżej pokazuje przebieg poszukiwań minimum tą właśnie metodą.
Metoda najszybszego spadku (the steepest descent method)
Jest to metoda nieco bardziej skomplikowana  najpierw, jak zwykle, wybieramy punkt startu i obliczamy gradient
badanej funkcji w tym punkcie. Wartości składowych gradientu określają nam linię (mającą kierunek gradientu). Jeśli
założymy, że x i y będą leżeć na tej właśnie linii to możemy badaną funkcję przekształcić do postaci parametrycznej:
P(x,y)
x x0
t
x
P(x,y)
y y0
t
y
P(x,y) P(x,y)
P(x,y) P x0 t , y0 t
x y
Badana funkcja będzie wtedy funkcją jednej zmiennej i będziemy mogli zastosować do niej znane już nam metody
szukania minimum funkcji jednej zmiennej (jest oczywiście pewien kłopot  nie będziemy znać właściwych granic
przedziału zmiennej t w jakim powinniśmy poszukiwać minimum  musimy przyjąć coś  z kapelusza ). Gdy już
znajdziemy wartość parametru t minimalizującą badaną funkcję to na jej podstawie obliczymy nowe współrzędne x i y.
P(x,y)
x1 x0
t
min
x
P(x,y)
y1 y0
t
min
y
Operację tę powtarzamy tak długo aż zadowolą nas wartości składowych gradientu  w minimum muszą oczywiście
wynosić zero, zatem to na ile różnią się od zera może stanowić pewną  miarę odległości od minimum. W metodzie tej nie
zwracamy zbytniej uwagi na zmiany wartości badanej funkcji  w każdym cyklu przesuwamy się na inną, trudną do
przewidzenia, odległość.
Przykładowy przebieg obliczeń w metodzie  steepest descent :
cykl: 1, wartość funkcji: 4.51207 błąd wzgl.: 0.595258
cykl: 2, wartość funkcji: 2.97947 błąd wzgl.: 0.0534011
cykl: 3, wartość funkcji: 2.83593 błąd wzgl.: 0.00265193
cykl: 4, wartość funkcji: 2.82944 błąd wzgl.: 0.000357943
cykl: 5, wartość funkcji: 2.82853 błąd wzgl.: 3.72817e  005
cykl: 6, wartość funkcji: 2.82844 błąd wzgl.: 4.22643e  006
cykl: 7, wartość funkcji: 2.82843 błąd wzgl.: 4.57242e  007
cykl: 8, wartość funkcji: 2.82843 błąd wzgl.: 4.26757e  008
cykl: 9, wartość funkcji: 2.82843 błąd wzgl.:  3.42853e  009
cykl: 10, wartość funkcji: 2.82843 błąd wzgl.:  8.53522e  009
Wyniki metody  steepest descent :
współrzędne: x_min=0.706976, y_min=0.707148, błędy wzgl. x i y:  5.01624e  005, 0.000193142,
wartość funkcji: 2.82843, błąd wzgl. wartości:  3.42853e  009, wartość prawdziwa: 2.82843,
wartości składowych gradientu: x-owa:  0.000249537, y-owa:  6.28107e  006.
Metoda ta może być łatwo zaadaptowana do przypadku funkcji o dowolnej liczbie zmiennych  niezależnie od liczby
zmiennych problem zawsze sprowadzamy do poszukiwania minimum funkcji jednej zmiennej (funkcji parametru t).
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 6
Metoda  rank-one
Jest to metoda zbliżona do metody najszybszego spadku, jednak oprócz gradientu (czyli pierwszych pochodnych)
wykorzystuje się w niej również drugie pochodne. Podobnie jak w metodzie najszybszego spadku wybieramy jakiś punkt
początkowy i obliczamy w nim wartości funkcji i jej gradientu. Zakładamy, że istnieje jakiś  właściwy kierunek , który
oznaczymy jako (s1,s2) i wzdłuż niego będziemy zmieniać wartości x i y (w metodzie najszybszego spadku kierunek ten
jednoznacznie wyznaczał nam gradient funkcji):
x a t s1
y b t s2
Zakładamy, że właśnie takie równania parametryczne można podstawić do wyrażenia na badaną funkcję P(x,y) (pomimo
faktu, że kierunków tych na razie nie znamy) po czym rozwijamy ją w szereg Taylora (dwuwymiarowy) wokół punktu o
współrzędnych (a,b):
f (x,y) f (a s1t,b s2t)
2 2 2 2
t 2 2
f (a,b) t s1 f s2 f s1 f 2s1s2 f s2 f
2 2
x y 2 x y
x y
Wszystkie pochodne w powyższym wzorze obliczane są w punkcie (a,b). Jeśli za s1 i s2 przyjmiemy wartości spełniające
następujący układ równań:
2 2
f
s1 f s2 f
2
x y x
x
2 2
f
s1 f s2 f
2
x y y
y
to dostaniemy:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 7
f (x,y) f (a,b)
2 2 2
1 2
2
t 1 1 s1 f 2s1s2 f s2 f
2 2
2 x y
x y
Będzie zatem zachodzić:
2 2 2
2 2
J e li |t| < 1 i s1 f 2s1s2 f s2 f > 0 to f(x,y) < f(a,b)
2 2
x y
x y
Spełnienie tych warunków jest dość łatwe jeżeli tylko szukając kolejnego nie oddalimy się zbytnio od poprzedniego
przybliżenia.
Przykładowy przebieg obliczeń w metodzie  rank-one :
cykl: 1, wartości: funkcji: -32.0051 błąd wzgl.:  12.3155, t_min: 0.499999,
cykl: 2, wartości: funkcji: -3.21376 błąd wzgl.:  1.13624, t_min: 0.235521
cykl: 3, wartości: funkcji: 2.82844 błąd wzgl.: 6.21862e  006, t_min: 0.196386
cykl: 4, wartości: funkcji: 2.82843 błąd wzgl.: 3.20727e  007, t_min: 0.230203
cykl: 5, wartości: funkcji: 2.82843 błąd wzgl.: 1.00811e  009, t_min: 0.280505
cykl: 6, wartości: funkcji: 2.82843 błąd wzgl.:  9.05838e  009, t_min: 0.362769
Wyniki metody  rank-one :
współrzędne: x_min=0.714811, y_min=0.714811, błędy wzgl. x i y: 0.0110318, 0.0110318,
wartość funkcji: 2.82859, błąd wzgl. wartości: 5.87003e  005, wartość prawdziwa: 2.82843,
wartości składowych gradientu: x-owa: 0.0214388, y-owa: 0.0214388.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 8
Adaptacja metody  rank-one do większej liczby wymiarów.
Zmiana polega na zastąpieniu układu dwóch równań służących w przypadku dwuwymiarowym do wyznaczenia
składowych kierunku poszukiwań (s1,s2) układem odpowiedniej (większej) liczby równań:
G s f
w przypadku dwóch wymiarów układ ten wygląda tak:
2 2
f f
f
2
s1
x y
x x
2 2
s2 f
f f
2
y
x y
y
dla trzech i więcej:
2 2 2
f f f
f
2
x y x z
x
s1 x
2 2 2
f f f
s2 f
2
x y y z y
y
s3 f
2 2 2
f f f
z
2
x z y z
z
(może to być całkiem spory układ równań).
Wszystkie pochodne byłyby oczywiście obliczane w aktualnym punkcie. Rozwiązanie tego układu powinno dać nam
wektor składowych kierunku poszukiwania następnego punktu (i tak właśnie postępowaliśmy w przypadku zagadnienia
dwuwymiarowego). Jednak w przypadku większej liczby wymiarów zamiast rozwiązywać dokładnie duży układ równań
posłużymy się przybliżeniem:
s H f
w którym macierz H ma zbiegać się do odwrotności G w miarę zbliżania się do lokalnego minimum, którego poszukujemy.
Na początek przyjmiemy, że jest to macierz jednostkowa (zawierająca jedynki na głównej przekątnej), w miarę postępu
obliczeń będziemy ją modyfikować  modyfikacja będzie polegać na dodawaniu poprawki E:
T
E Ä…u u
Poprawkę tę będziemy wyznaczać następująco. Najpierw zauważamy, że (jeśli drugie pochodne obliczamy numerycznie)
możliwy jest zapis:
G s f G "x "f
Wówczas, jeśli dodatkowo założymy (o czym była mowa już wcześniej), że kolejne poprawki E zbliżają macierz H do
odwrotności macierzy G otrzymamy:
T
H E "f "x H Ä…uu "f "x
które to równanie będzie spełnione jeśli zajdzie równocześnie:
T
u "x H"f i Ä…u "f 1
Tryb postępowania podczas poszukiwania minimum tą metodą będzie następujący:
1 Wybieramy punkt poczÄ…tkowy a,
2 Obliczamy wektor pochodnych cząstkowych f a (jeśli trzeba numerycznie).
3 Poszukujemy minimum funkcji jednej zmiennej t (parametru): f t f a tH f a
4 Zastępujemy a przez a t f a
min
5 Koniec obliczeń następuje gdy uznamy, że zmiany a są wystarczająco małe
6 Wyznaczamy nowÄ… macierz H:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 9
6.1 wektor u "x H"f a1 a0 H f a1 f a0
T
6.2 współczynnik ą wyznaczamy z równania: ą u "f 1
T
6.3 cała poprawka: E ą uu
T
6.4 i nowa macierz: H1 H0 Ä… uu
Programowanie liniowe
Określenie zagadnienia
Zagadnienia programowania liniowego polegajÄ… na znalezieniu maksimum zadanej liniowej funkcji wielu zmiennych
przy uwzględnieniu wielu dodatkowych warunków narzuconych na te zmienne.
Metoda znana jako  programowanie liniowe zaczęła być powszechnie stosowana podczas II Wojny Światowej do
rozwiązywania zagadnień logistycznych (w armiach angielskiej i amerykńskiej). Nazwa ta jest uzasadniona co najwyżej
historycznie, bo w swej oryginalnej postaci metoda ta nie miała nic wspólnego z programowaniem komputerów.
Przykład sformułowania problemu
Załóżmy, że mamy fabrykę produkującą dwa typy płaszczy: zwykłe przeciwdeszczowe i luksusowe. Cały proces
produkcji obejmuje trzy procesy cząstkowe: krojenie materiałów, szycie i pakowanie. Płaszcze produkuje się partiami (po
ileś tam sztuk)  jednak ze względu na uwarunkowania lokalne (dostępność wykwalifikowanego personelu itp.) liczba
godzin, które można przeznaczyć na każdy z procesów jest ograniczona. Konkretnie: na cięcie można przeznaczyć 42
godziny (roboczogodziny), na szycie  60 godzin i na pakowanie  32 godziny. Ilość czasu jaką trzeba przeznaczyć na
każdą z czynności jest inna dla każdego typu płaszcza. Wycinanie materiału na płaszcz przeciwdeszczowy zajmuje 8 godzin,
podczas gdy wycinanie płaszczy luksusowych zajmuje 5 godzin. Szycie płaszcza zwykłego zajmuje 3 godziny, a
luksusowego 16. Pakowanie: 5 godzin dla partii płaszczy zwykłych i 6 godzin dla luksusowych. Na każdej partii płaszczy
zwykłych zarabiamy 500 zł, a na każdej partii płaszczy luksusowych  400 zł.
Dodatkowo zakładamy, że wszyscy pracownicy otrzymują takie samo wynagrodzenie (niezależnie od tego czy pracują
czy nie, że materiały ( substraty ) i (ewentualne) narzuty kosztują tyle samo (niezależnie od procesu) oraz że fabryka może
sprzedać wszystko co wyprodukuje. Chcemy wiedzieć w jaki sposób zmaksymalizować zysk.
Dane:
zwykły luksusowy dostępny
przeciwdeszczowy czas pracy
cięcie 8 5 42
szycie 3 16 60
pakowanie 5 6 32
zysk 500 400
Wprowadz
my oznaczenia: x1  liczba partii zwykłych płaszczy przeciwdeszczowych, jakie fabryka może
wyprodukować w ciągu dnia, x2  analogiczna liczba dla płaszczy luksusowych.
Tzw.  funkcją celu (którą chcemy zmaksymalizować) będzie zysk: z = 500x1 + 400x2. Mamy też ograniczenia
(dotyczÄ…ce czasu pracy):
8x1 + 5x2 42 (cięcie)
3x1 + 16x2 60 (szycie)
5x1 + 6x2 32 (pakowanie)
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 10
Ddatkowo dodamy jeszcze ograniczenia oczywiste (zwane trywialnymi):
x1 0 i x2 0
W przypadku tak prostych problemów (pod względem liczby zmiennych) można pozwolić sobie na zastosowanie
metody graficznej. Posłuży ona jako ilustracja i umożliwi póz
niejsze wprowadzenie metody numerycznej.
Metoda graficzna
Najpierw wykreślamy nasz  obszar zainteresowania  tzn. obszar, w którym wartości obu zmiennych spełniają
narzucone na nie warunki.
Na tym samym wykresie możemy też zaznaczyć linie odpowiadające kilku wybranym wartościom zysku  w ten sposób
możemy szybko zorientować się uzyskanie jakich zysków w ogóle wchodzi w grę, a jakie wartości są niemożliwe do
uzyskania. Na koniec, również z rysunku, łatwo jest zorientować się, którędy będzie przebiegać linia zysku odpowiadająca
wartości maksymalnej  jeśli już określimy przez który  narożnik obszaru ma przechodzić zauważymy, przez punkt
wspólny których dwóch warunków musi przechodzić linia zysku i wyznaczymy wartości x1 i x2 odpowiadające
maksymalnemy zyskowi, a potem i sam zysk. W naszym przykładzie maksymalizują zysk wartości x1=4 i x2=2, a odpowiada
im zysk wynoszÄ…cy 2800.
Przypadek ogólniejszy
Jak łatwo zauważyć parametry zostały  ręcznie tak dobrane aby otrzymać wynik wyrażający się liczbami całkowitymi
 gdyby wynikiem nie były liczby całkowite moglibyśmy mieć spory problem  np. ze względów technicznych nie byłoby
możliwe wykonywanie tylko części partii towaru (przykładem może być choćby i huta i (niemożliwa do uzyskania)
 ułamkowa część wytopu z wielkiego pieca).
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 11
Wracając do naszego przykładu  był on bardzo prymitywny: różni pracownicy otrzymują różne wynagrodzenia, do
wytworzenia jednego wyrobu potrzeba wielu różnych materiałów (o różnych cenach), może dojść jeszcze zużycie energii,
koszty utylizacji odpadków itd., itp.. W rzeczywistości problem nigdy nie będzie tylko dwuwymiarowy  liczba wymiarów
będzie o wiele większa  zatem metody graficznej nie da się zastosować. Jednak prawdziwy pozostanie fakt, że
rozwiązaniem będzie wierzchołek wielowymiarowej bryły (będącej naszym obszarem zainteresowania), który będzie
styczny do  linii zysku (albo ścieślej do czegoś  o jeden mniej wymiarowego niż nasz obszar zainteresowania).
 Na oko wydaje się, że rozwiązanie można by znalez
ć badając wartości zysku odpowiadające wszystkim  narożnikom
 jednak okazuje się, że w przypadku np. 16 zmiennych i 20 warunków (tych  nietrywialnych , tzn. pomijając oczywiste
warunki głoszące, że dana zmienna ma być większa od zera) otrzymujemy ok. 7 milionów narożników (a jest to problem,
wg standardów programowania liniowego, bardzo mały). Potrzebne jest więc inne podejście.
Nosi ono nazwę metody  simplex , a jej zasadę działania można zilustrować za pomocą adaptacji do przypadku
wielowymiarowego metody, którą posłużyliśmy się w (dwuwymiarowym) przypadku rozwiązywanym graficznie.
W największym skrócie i uproszczeniu metoda ta działa następująco: wybieramy jakiś wierzchołek (choćby i ten dla
wszystkich zmiennych wynoszących zero), sprawdzamy jaka jest w nim wartość naszej funkcji celu (zysku), następnie
spośród sąsiednich wierzchołków wybierzemy ten, dla przejścia do którego przyrost funkcji celu będzie maksymalny,
przejdziemy do niego i z kolei wśród jego sąsiadów będziemy szukać wierzchołka dającego maksymalny przyrost funkcji
celu. Operacje takie będziemy powtarzać aż do znalezienia się w wierzchołku, którego wszyscy sąsiedzi dają ujemny
przyrost funkcji celu  uznamy wówczas problem za rozwiązany.
Oczywiście aby operację tę mógł samodzielnie przeprowdzić komputer musimy odpowiednio przygotować dane i
sformułować algorytm rozwiązywania takich problemów dla dowolnej liczby zmiennych i dowolnej liczby ograniczeń.
Postać kanoniczna problemu
Najpierw określmy postać standardową zagadnienia  odpowiada ona uporządkowaniu wszystkich nierówności
opisujÄ…cych zagadnienie:
z c1x1 c2x2 c x
n n
c11x1 c12x2 c1x b1
n
c21x1 c22x2 c2x b2
n
c x1 c x2 c x b
m1 m2 mn n m
x1,x2, ,x 0
n
Warto zauważyć, że w postaci standardowej wszystkie warunki  nietrywialne zapisane są w postaci  mniejszy lub
równy  zatem podczas przepisywania zagadnienia do postaci standardowej może być konieczne odwrócenie niektórych
nierówności (pomnożenie przez  1, z czym będzie się wiązać zmiana znaków współczynników na przeciwne). Natomiast
wszystkie warunki trywialne mają postać  większy lub równy (to zazwyczaj nie będzie ulegać zmianie).
W postaci kanonicznej zawsze mówimy o maksymalizacji funkcji, a warunki nietrywialne muszą mieć postać
równości  co można uzyskać przez wprowadzenie zmiennych dodatkowych.
Rozważmy to na następującym przykładzie  mamy zminimalizować funkcję z:
z x1 3x2 4x3
2x1 4x3 4
x1 x2 5x3 4
2x2 x3 2
x1,x2,x3 0
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 12
W postaci standardowej zagadnienie to przybierze postać:
z x1 3x2 4x3
2x1 4x3 4
x1 x2 5x3 4
2x2 x3 2
x1,x2,x3 0
Znak funkcji z został zmieniony abyśmy mogli mówić o maksymalizacji. Aby zamienić tę postać na kanoniczną musimy
sprowadzić nierówności będące warunkami (ograniczeniami) nietrywialnymi do równości  zrobimy to wprowadzając
dodatkowe zmienne:
2x1 4x3 4
2x1 4x3 x4 4
x4 0
Jeśli przeprowadzimy taką operację dla wszystkich warunków otrzymamy:
z x1 3x2 4x3
2x1 4x3 x4 4
x1 x2 5x3 x5 4
2x2 x3 x6 2
x1,x2,x3,x4,x5,x6 0
co jest właśnie postacią kanoniczną naszego zagadnienia.
Przykład realizacji metody simplex
Rozpatrzmy działanie takiej metody na poprzednim przykładzie (produkcja płaszczy). Problem maksymalizacji zysku
z produkcji płaszczy ma, w postaci kanonicznej, następującą postać:
z 500x1 400x2
8x1 5x2 x3 42
3x1 16x2 x4 60
5x1 6x2 x5 32
x1,x2,x3,x4,x5 0
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 13
Na rysunku powyżej widzimy nasz  obszar zainteresowania przedstawiony jako obszar ograniczony prostymi o
równaniach postaci xi = 0. Istotne dla nas wierzchołki, w których będziemy badać zmiany funkcji celu oznaczone są na
rysunku wielkimi literami A...E, natomiast wierzchołki, które wprawdzie również stanowią punkty przecięcia par prostych
xi = 0, ale jako nie należące do krawędzi interesującego nas obszaru są z naszego punktu widzenia nieistotne, oznaczone
są jako X, Y, Z (są jeszcze dwa takie wierzchołki, ale znalazły się poza wykresem).
Można zauważyć, że w każdym z wierzchołków pokazanych na wykresie powyżej dwie zmienne mają wartość zero,
a trzy pozostałe nie.
W ogólnym przypadku, jeśli problem był określony przez n zmiennych, do których dodaliśmy jeszcze m zmiennych
dodatkowych (mamy ich więc w sumie n+m) to w każdym z wierzchołków n spośród zmiennych przyjmuje wartość zero.
Zmienne, które w danym wierzchołku mają wartości niezerowe określamy jako podstawowe (w danym wierzchołku), a
pozostałe jako niepodstawowe (w tymże wierzchołku). Zbiór wartości wszystkich zmiennych (podstawowych i
niepodstawowych) w danym wierzchołku określamy jako rozwiązanie odpowiadające temu wierzchołkowi. W tabeli
pokazano, które zmienne są podstawowe, a które nie, w poszczególnych wierzchołkach w naszym przykładzie:
wierzchołek zmienne zmienne
podstawowe niepodstawowe
A x1, x2, x4 x3, x5
B x1, x4, x5 x2, x3
C x3, x4, x5 x1, x2
D x2, x3, x5 x1, x4
E x1, x2, x3 x4, x5
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 14
Można zauważyć, że przemieszczając się od wierzchołka do wierzchołka po obwodzie interesującego nas obszaru za
każdym razem zmieniamy jedną ze zmiennych podstawowych na niepodstawową i odwrotnie: A B x2 x5, B C
x1 x3, C D x2 x4, D E x1 x5, E A x3 x4.
Rozwiązanie metodą simplex rozpoczynamy od nadania wszystkim zmiennym występującym w oryginalnym problemie
wartości zerowych. Jeśli wierzchołek ten należy do krawędzi obszaru zainteresowania to wszystkie wartości zerowe są
dopuszczalnym rozwiązaniem naszego problemu (zazwyczaj tak jest) i możemy rozpocząć obliczenia  jeśli nie 
będziemy musieli znalez
ć inny wierzchołek, od którego rozpoczniemy (o tym jak to zrobić  póz
niej).
W przykładowym problemie rozwiązaniem początkowym będzie x1=0 i x2=0 (jak widać na rysunku takie rozwiązanie
jest dopuszczalne  jest to punkt C). W przypadku gdy nie mamy rysunku (np. problem jest wielowymiarowy)
podstawiamy założone wartości zerowe do warunków problemu (sformułowanego w postaci kanonicznej  a więc już z
wprowadzonymi zmiennymi dodatkowymi) i sprawdzamy czy warunki te nie są sprzeczne. Jeśli nie są  rozwiązanie jest
dopuszczalne. A
atwo zauważyć (patrząc na funkcję celu z), że zwiększanie x1 (od wartości zerowej) spowoduje większy
przyrost funkcji celu (zysku  z) niż analogiczne zwiększanie wartości x2  poszukamy więc zmiennej niepodstawowej,
która ma największy współczynnik w (maksymalizowanej przez nas) funkcji celu (w wierzchołku od którego rozpoczęliśmy
zmiennymi niepodstawowymi były właśnie x1 i x2). Zdecydowawszy się na x1 pozwalamy jej przyjąć wartość niezerową (a
więc stać się zmienną podstawową)  musimy określić, która z dotychczasowych zmiennych podstawowych ma stać się
niepodstawową. W tym celu równania określające warunki przepiszemy w następującej postaci:
x3 42 8x1 5x2
x4 60 3x1 16x2
x5 32 5x1 6x2
Sprawdzamy teraz, która z pozostałych zmiennych pierwsza osiągnie zero gdy wartość zmiennej x1 będzie wzrastać
ponad zero  w tym celu wystarczy podzielić wolny wyraz w każdym z równań przez współczynnik stojący przy x1 
otrzymamy: 42/8 (5.25) dla x3, 60/3 (20) dla x4 i 32/5 (6.4) dla x5  najszybciej spadnie do zera zmienna x3 i tą też zmienną
wybieramy jako  zamiennik dla x1. Z naszej tabeli wynika, że przeniesiemy się do wierzchołka B.
W wierzchołku B powtarzamy całą procedurę  sprawdzamy, która ze zmiennych niepodstawowych stając się
podstawową będzie mieć największy wpływ na wzrost funkcji celu z. W wierzchołku tym (B) zmiennymi niepodstawowymi
są x2 i x3  funkcję celu przepiszemy więc z ich wykorzystaniem:
42 5x2 x3
x1
8
42 5x2 x3
z 500 400x2
8
2625 87.5x2 62.5x3
Wynika stąd, że jeśli zmienna x2 stanie się podstawową to będzie to mieć większy wpływ na wzrost z niż taka zmiana
w przypadku x3 (właściwie, ponieważ współczynnik przy zmiennej x3 jest ujemny, moglibyśmy wyłączyć ją z tych rozważań
 tylko jej wartości ujemne mogłyby wpływać na wzrost z, a te są  z definicji wykluczone). Warunki przepiszemy w
następującej postaci (pamiętając, że x3 ma pozostać niepodstawową  czyli wynoszącą zero):
42 5x2
x1
8
3x1 x4 60 16x2
5x1 x5 32 6x2
Jak widać pierwszego z warunków nie trzeba już przekształcać (odpowiednią formę uzyskał gdy wyznaczaliśmy wartość
zmiennej x1), pozostałe dwa przekształcamy (wykorzystując wyrażenie na x1):
42 5x2
x1
8
354 113x2
x4
8
46 23x2
x5
8
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 15
Ponownie sprawdzamy współczynniki (wolny wyraz podzielony przez współczynnik stojący przy x2)  najmniejszym
jest 46/23 (odpowiadający zmiennej x5)  czyli to zmienna x5 pierwsza osiągnie zero gdy x2 będzie wzrastać od zera. Zatem
zmienna x5 stanie się niepodstawową, zaś zmienna x2 stanie się podstawową   przenosi nas to z wierzchołka B do
wierzchołka A.
Znaną już procedurę powtarzamy w wierzchołku A. Znowu przekształcamy funkcję celu tak aby wyrażona była przez
zmienne niepodstawowe (w tym przypadku/wierzchołku są to x3 i x5). Wychodząc od ostatnio zmodyfikowanej postaci
funkcji celu otrzymujemy i przekształcając ostatnie z powyższych równań (wiążące x5 z x2):
z 2625 87.5x2 62.5x3
8x5 46
2625 87.5 62.5x3
23
700
2800 62.5x3 x5
23
Jak widać doszliśmy do momentu, gdy nie jesteśmy w stanie zwiększyć wartości funkcji celu z ponieważ obie
występujące we wzorze zmienne podstawowe mają ujemne współczynniki. Stwierdzamy więc, że znalez
liśmy rozwiązanie
optymalne  konkretnie z = 2800. Podstawienie x3=0 i x5=0 we wzorach na warunki zagadnienia daje:
8 x1 5x2 42
5x1 6x2 32
Skąd łatwo otrzymujemy x1=4 i x2=2 (czyli ten sam wynik, który otrzymaliśmy metodą graficzną).
Warto zwrócić tu uwagę, że w bardziej skomplikowanych przypadkach możliwe jest  zapętlenie rozwiązania (gdy
proces zacznie przebiegać sekwencję wierzchołków nie zwiększając wartości funkcji celu).
Początkowe przybliżenie rozwiązania
Problem znajdowania początkowego przybliżenia rozwiązania zilustrujemy na nieco zmodyfikowanym przykładzie
fabryki płaszczy. Dodamy jeszcze jeden warunek  fabryka musi produkować przynajmniej jedną partię płaszczy
przeciwdeszczowych dziennie.
Po uwzględnieniu tego dodatkowego warunku postać kanoniczna problemu jest następująca:
z 500x1 400x2
8x1 5x2 x3 42
3x1 16x2 x4 60
5x1 6x2 x5 32
x1 x6 1
x1,x2,x3,x4,x5,x6 0
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 16
A przedstawienie graficzne wygląda tak jak na rysunku obok. Jak widać zastosowane wcześniej rozwiązanie
początkowe, w którym obie zmienne występujące w oryginalnym sformułowaniu problamu wynoszą zero nie jest
dopuszczalne  dotychczasowy punkt C nie znajduje się już na krawędzi obszaru zainteresowania  założenie x1=0 łamie
warunek x6 0.
Aby rozpocząć jednak w jakiś sposób obliczenia zaczynamy od funkcji pseudocelu z . Funkcję tę wprowadzamy mając
na celu zwiększenie wartości zmiennej x6 i doprowadzenia jaj do wartości dodatniej. Przeprowadzimy w tym celu jeden
cykl metody simplex dla celu danego równaniem:
z = x6 = x1  1
i z pozostałymi warunkami oryginalnymi. Zgodnie z zasadami stwierdzamy, że zmiana x1 będzie mieć większy wpływ na
z niż zmiana x2. Wobec tego zmienna x1 staje się zmienną podstawową  zmiennymi podstawowymi są więc obecnie x3,
x4, x5 i x6, zatem utrzymując dla x2 wartość zero mamy:
x3 42 8x1
x4 60 3x1
x5 32 5x1
x6 1 x1
Widać, że przy wzroście x1 ponad zero x6 będzie pierwszą zmienną, która osiągnie zero i właśnie x6 stanie się
niepodstawową. Mamy teraz x2=0 i x6=0, skąd wynika, że x1=1, x3=34, x4=57 i x5=27.
To jest właśnie wyjściowe przybliżenie rozwiązania i od niego można rozpocząć dalsze rozwiązywanie zagadnienia
metodÄ… simplex.
Przeprowadzony jeden cykl z funkcją pseudocelu  przeniósł nas z punktu C do punktu Y (zaznaczonego na rysunku
powyżej).
Gdyby uzyskany punkt nie był dopuszczalnym rozwiązaniem musielibyśmy wybrać inną funkcję pseudocelu  mogłaby
to być np. suma wszystkich zmiennych niedodatnich.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Minima funkcji wielu zmiennych 17
Problemy z rozwiązaniami będącymi liczbami całkowitymi
W rozwiązanym przykładzie otrzymaliśmy wynik wyrażony liczbami całkowitymi  bardzo nas to satysfakcjonuje. W
wielu przypadkach rozwiązania będące liczbami całkowitymi są jedynymi dopuszczalnymi  nie sposób sobie wyobrazić
wyprodukowania ułamka zawartości wielkiego pieca stali czy  gruszki betonu  istnieją procesy produkcyjne, w których
nie można zostawić części produkcji  na jutro .
Znajdowanie rozwiązań całkowitych może znakomicie utrudnić rozwiązywania takich zagadnień. Można sobie z tym
jednak poradzić stosując następujące podejście:
Na początku rozwiązujemy zagadnienie tak samo jak w przedstawionych przykładach  czyli tak jaby zadowalało nas
znalezienie rozwiązań ułamkowych. Gdy znajdziemy już rozwiązanie optymalne sprawdzamy czy wyraża się ono liczbami
całkowitymi  jeśli tak  mamy problem z głowy . Jeśli nie  wybieramy tę zmienną, która w rozwiązaniu nie jest
całkowita (jeśli jest ich kilka to wybieramy którąkolwiek) i zaczynamy rozwiązywać dwa problemy  w każdym z nich
do oryginalnego zestawu warunków dodajemy jeden dodatkowy ograniczający zmienność tej właśnie zmiennej 
uniemożliwiający jej przyjęcie tej właśnie wartości ułamkowej.
Gdyby na przykład zmienna x1 w obliczonym rozwiązaniu przyjmowała wartość 3.3 przystąpilibyśmy do rozwiązywania
dwóch  nowych zagadnień  jedno z nich zostałoby  wzbogacone o warunek x1 3, a drugie o warunek x1 4. Procedurę
taką powtarzalibyśmy (zapisując otrzymywane rozwiązania całkowite  o ile oczywiście takie byśmy  po drodze
uzyskiwali) tak długo aż wyczerpalibyśmy wszystkie możliwości. Następnie spośród zanotowanych rozwiązań całkowitych
wybralibyśmy najlepsze (maksymalizujące funkcję celu).
Należy zwrócić uwagę, że dodanie warunków do nowych  podproblemów ogranicza  obszar zainteresowania dla
danego podproblemu  jeśli podczas rozwiązywania któregoś z podproblemów znalez
liśmy rozwiązanie optymalne dla
danego  podobszaru to rozwiązywanie kolejnych podproblemów dotyczących tego samego podobszaru mija się z celem
 umożliwia bowiem co najwyżej powtórne znalezienie tego samego rozwiązania.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej