Masarykova univerzita Přírodovědecká fakulta BAKALÁŘSKÁ PRÁCE Jan Kovář Algoritmické řešení úlohy lineárního programování v rovině Vedoucí práce: doc. RNDr. Martin Čadek, CSc. Studijní program: Aplikovaná matematika Studijní obor: Matematika - ekonomie 2011 Na tomto místě bych rád poděkoval doc. RNDr. Martinu Čadkovi, CSc. za ochotu, laskavý přístup a cenné rady poskytnuté během odborného vedení této práce. Prohlašuji, že jsem svou bakalářskou práci napsal samostatně a výhradně s použitím citovaných pramenů. V Brně dne Jan Kovář Název práce: Algoritmické řešení úlohy lineárního programování v rovině Autor: Jan Kovář Ústav matematiky a statistiky Přírodovědecké fakulty, MU Vedoucí bakalářské práce: doc. RNDr. Martin Čadek, CSc. Klíčová slova: lineární programování, náhodnostní algoritmus, optimalizace Abstrakt Práce se zabývá maximalizací lineární funkce na průniku konečně mnoha polorovin. Je odvozen náhodnostní algoritmus k řešení tohoto problému, který je založen na iteračním principu postupného přidávání polorovin. Následně je analyzována jeho očekávaná časová složitost. Je dokázáno, že roste lineárně vzhledem k počtu polorovin obsažených v úloze. Součástí práce je také implementace algoritmu. Title: Algoritmic approach to linear programming in the plane Author: Jan Kovář Department of Mathematics and Statistics, Faculty of Science, MU Supervisor: doc. RNDr. Martin Čadek, CSc. Keywords: linear programming, randomized algorithm, optimization Abstract The work deals with maximizing a linear function at the intersection of a finite number of half-planes. A randomized algorithm solving this problem is derived. It is based on the principle of iteration, when the half-planes are progressively added. Subsequent analysis of the expected time complexity is provided. The proof that the expected running time increases linearly due to the number of half-planes contained in the problem is given. The work also includes an implementation of the algorithm. Obsah Úvod 6 1 Základní pojmy 7 2 Jednodimenzionální úloha 11 3 Algoritmus pro omezenou úlohu lineárního programování 14 4 Neomezená úloha lineárního programování 21 5 Očekávaná časová složitost 35 6 Simplexová metoda 39 7 Úloha lineárního programování ve vyšších dimenzích 41 8 Implementace 43 Seznam použité literatury 46 5 Úvod Úloha lineárního programování je speciálním případem hledání vázaných extrémů reálné funkce více proměnných. Jedná se o situaci, kdy vyšetřovaná funkce je lineární a snažíme se ji maximalizovat na množině, která je zadána soustavou lineárních nerovnic. V rovině pak hledáme bod ležící v průniku konečně mnoha daných polorovin, ve kterém lineární funkce dvou proměnných nabývá maxima. Nejznámějším algoritmem k řešení úlohy lineárního programování je simplexová metoda. Zde se však budeme zabývat jiným algoritmem. Budeme vycházet především z [1, kapitola 4]. Jedná se o náhodnostní algoritmus. To znamená, že průběh algoritmu nezávisí pouze na jeho vstupu, ale také na nějaké náhodné události. Na začátku budou vysvětleny základní pojmy a souvislosti problému. Následně bude odvozen algoritmus pro jednodimenzionální úlohu lineárního programování, a to především proto, že jeho různé modifikace jsou použitelné při řešení úlohy lineárního programování v rovině. V další části se budeme zabývat omezenou úlohou v rovině, pro kterou bude uveden detailní postup odvození algoritmu. Posléze bude analyzován problém omezenosti úlohy. Odvodíme algoritmus k rozpoznání neomezené úlohy. Postup pro řešení omezené úlohy využijeme ke konstrukci algoritmu pro obecnou úlohu lineárního programování. Následně vysvětlíme koncept očekávané časové složitosti a dokážeme, že očekávaná časová složitost odvozeného algoritmu je lineární vzhledem k počtu polorovin obsažených v úloze. V dalších dvou kapitolách uvedeme srovnání se simplexovou metodou a možné rozšíření algoritmu pro úlohy lineárního programování, které neřešíme v rovině, ale v prostorech vyšších dimenzí. Přínosem práce oproti základnímu zdroji [1, kapitola 4] je detailní zpracování, a to jak postupného odvození, tak výsledných algoritmů. Uvedené algoritmy jsou proto snadno implementovatelné. Navíc byl odhalen případ neexistence lexikograficky nejmenšího řešení, který autoři [1] nevzali v úvahu. Součástí práce je také implementace algoritmu. Program byl vytvořen v prostředí Delphi 6.0 a zahrnuje mimo jiné přehledný grafický výstup. Důraz je kladen na demonstraci principu algoritmu, a proto je aplikace vhodná pro výukové účely. Funkce programu jsou popsány v poslední kapitole. 6 Kapitola 1 Základní pojmy Definice 1.1. Nechť I = {1, 2, . . . , n}, J = {1, 2, . . . , d}, I ⊆ I, J ⊆ J. Úlohou lineárního programování rozumíme úlohu najít bod (x1, . . . , xd) ∈ Rd , ve kterém nabývá lineární funkce c1x1 + c2x2 + · · · + cdxd svého maxima (minima) za podmínek ai1x1 + ai2x2 + · · · + aidxd ≤ bi, i ∈ I \ I, ai1x1 + ai2x2 + · · · + aidxd = bi, i ∈ I, xj ≥ 0, j ∈ J, kde aij, bi, ci ∈ R pro každé i ∈ I, j ∈ J a x1, x2, . . . , xd jsou proměnné. Lze jednoduše ukázat, že libovolnou úlohu lineárního programování lze převést na ekvivalentní úlohu lineárního programování ve tvaru1 dle definice 1.2. Definice 1.2. Úlohou lineárního programování rozumíme úlohu najít bod (x1, . . . , xd) ∈ Rd , ve kterém nabývá lineární funkce c1x1 + c2x2 + · · · + cdxd svého maxima za podmínek a11x1 + a12x2 + · · · + a1dxd ≤ b1 a21x1 + a22x2 + · · · + a2dxd ≤ b2 ... an1x1 + an2x2 + · · · + andxd ≤ bn, kde c1, . . . , cd, a11, . . . , and, b1, . . . , bd ∈ R a x1, x2, . . . , xd jsou proměnné. Označíme-li A = (aij) pro i = 1, . . . , n, j = 1, . . . , d, b = (b1, . . . , bn)T , c = (c1, . . . , cd), x = (x1, . . . , xd)T , pak lze pomocí vektorů úlohu zkráceně zapsat max {cx | Ax ≤ b} . 1 Požadujeme-li navíc, aby všechny proměnné byly nezáporné, nazývá se tento tvar kanonický. 7 KAPITOLA 1. ZÁKLADNÍ POJMY Funkce f(x) = f(x1, x2, . . . , xd) = c1x1 +c2x2 +· · ·+cdxd se nazývá účelová funkce. Řekneme, že x je přípustné řešení úlohy lineárního programování, jestliže platí Ax ≤ b. Optimálním řešením nazveme takové přípustné řešení x∗ , které maximalizuje účelovou funkci, tedy x∗ je optimální řešení, jestliže Ax∗ ≤ b a pro každé x splňující Ax ≤ b platí cx ≤ cx∗ . Zde se budeme zabývat maximalizační úlohou lineárního programování v rovině, tedy speciálním případem výše uvedené definice pro d = 2. Proměnné budeme značit x a y, složky vektoru c označme cx, cy. Pak účelová funkce je tvaru f(x, y) = cxx + cyy. V následujícím textu budeme automaticky předpokládat, že každá úloha lineárního programování je zadána ve tvaru dle definice 1.2. Navíc předpokládejme, že pro každé i ∈ {1, . . . , n} je alespoň jeden z koeficientů ai1, ai2 nenulový. V opačném případě by daná nerovnice buď byla splněna vždy, a proto byla nadbytečná, anebo by nebyla splněna nikdy, a tudíž by úloha neměla přípustné řešení. Dále předpokládejme, že účelová funkce je nenulová, jinak by úloha neměla smysl. Řešení (x, y)T můžeme reprezentovat jako body v dvourozměrném euklidovském prostoru. Z analytické geometrie víme, že rovnice ai1x + ai2y = bi zadává v R2 přímku a nerovnice ai1x+ai2y ≤ bi polorovinu, jejíž hraniční přímkou je tato přímka. Množina všech přípustných řešení úlohy lineárního programování je tedy průnikem n polorovin. Jelikož polorovina je konvexní oblast a průnik konvexních oblastí je konvexní oblast, je množina všech přípustných řešení konvexní oblastí. c x y f(x, y) = c2 x + c2 y > 0 f(x, y) = 0 cx cy Obrázek 1.1: Vrstevnice účelové funkce Účelová funkce f je lineární funkcí dvou proměnných. Zobrazíme-li v R2 její vrstevnice (tedy křivky dané rovnicemi cxx + cyy = k pro k ∈ R), dostaneme rovnoběžné přímky, jejichž normálovým vektorem je vektor c. Zřejmě cxx + cyy = 0 je rovnice vrstevnice procházející počátkem. V každém bodě ležícím na této přímce je hodnota účelové funkce nulová. Uvažme vrstevnici procházející bodem [cx, cy]. Její rovnice má tvar cxx+cyy = c2 x +c2 y. Hodnota účelové funkce v bodech ležících na této 8 KAPITOLA 1. ZÁKLADNÍ POJMY přímce je rovna c2 x+c2 y, a tedy je kladná, což znamená, že je větší než hodnota účelové funkce v bodech vrstevnice procházející počátkem. Vektor c proto udává směr, ve kterém účelová funkce roste, tedy směr, ve kterém se snažíme o maximalizaci. Pro každé i ∈ {1, . . . , n} označme poloroviny hi = {(x, y) ∈ R2 | ai1x+ai2y ≤ bi}. Nechť H = {h1, h2, . . . hn}. Pro úlohu lineárního programování pak budeme používat označení (H, c), zkráceně též pouze H, bude-li z kontextu jasné, o jaký vektor c se jedná. V závislosti na množině přípustných řešení a existenci optimálního řešení můžeme rozlišit následující 4 situace. • Úloha nemá přípustné řešení (tedy ani optimální řešení), což nastává, pokud je průnik polorovin prázdný. • Úloha má jediné optimální řešení. • Úloha má více optimálních řešení. Tato situace může nastat v případě, kdy je vektor c kolmý na přímku tvořící hranici konvexní množiny přípustných řešení. • Množina přípustných řešení je neprázdná, ale úloha nemá optimální řešení. Tato úloha lineárního programování se nazývá neomezená. V tomto případě existuje polopřímka ρ, která je celá obsažena v množině přípustných řešení a podél níž hodnota účelové funkce neomezeně roste. Příklady jednotlivých situací jsou znázorněny na obrázku 1.2. Na které straně od hraniční přímky polorovina leží, je vyznačeno vystínováním. V této práci budeme prezentovat různé algoritmy, stručně se proto podívejme na pojmy související se složitostí algoritmů. Vycházejme přitom z [2]. Každému vstupu můžeme přiřadit jeho velikost. V našem případě, kdy na vstupu je úloha lineárního programování, budeme za velikost vstupu považovat n, tedy počet nerovnic určujících množinu všech přípustných řešení. Definice 1.3. Délkou výpočtu algoritmu pro daný vstup rozumíme počet vykonaných elementárních operací. Nechť f : N → N je funkce. Řekneme, že algoritmus A má časovou složitost f(n), jestliže pro každé n ∈ N je f(n) maximální délka výpočtu pro všechny vstupy velikosti n. Při analýze časové složitosti algoritmu je však rozhodující asymptotické chování, zanedbáváme proto aditivní a multiplikativní konstanty. Zajímá nás tedy, jak rychle funkce roste s rostoucí velikostí vstupu. Definice 1.4. Řekneme, že funkce g : N → N roste nejvýše tak rychle jako funkce f : N → N, píšeme g = O(f), jestliže existuje c ∈ R a n0 ∈ N takové, že pro každé n ∈ N, n ≥ n0 platí g(n) ≤ cf(n). Při analýze časové složitosti pak říkáme, že algoritmus A má časovou složitost O(f(n)), což znamená, že časová složitost algoritmu A roste nejvýše tak rychle jako funkce f. Speciálně, jestliže algoritmus A má časovou složitost O(n), říkáme že má lineární časovou složitost; pokud má složitost O(n2 ), říkáme že má kvadratickou časovou složitost. 9 KAPITOLA 1. ZÁKLADNÍ POJMY c (a) Prázdná množina přípustných řešení c (b) Jediné optimální řešení c (c) Více optimálních řešení ρ p c (d) Neomezená úloha Obrázek 1.2: Varianty množin všech přípustných řešení a optimálních řešení 10 Kapitola 2 Jednodimenzionální úloha Než začneme odvozovat algoritmus pro řešení úlohy lineárního programování v rovině, podívejme se na jednodimenzionální úlohu, tedy úlohu na přímce. V následujícím textu se s touto úlohou několikrát setkáme, a proto je vhodné rozebrat ji zvlášť v této části. Nechť je dána jednodimenzionální úloha lineárního programování, tedy úloha najít bod x ∈ R, ve kterém účelová funkce f(x) = cx nabývá svého maxima za podmínek a1x ≤ b1, ... anx ≤ bn, kde c, a1, . . . , an, b1, . . . , bn ∈ R, c = 0. Množina všech přípustných řešení je průnikem polopřímek. Pro libovolné j ∈ {1, . . . , n} uvažme j-tou nerovnici ajx ≤ bj. Mohou nastat následující případy: • Jestliže aj > 0, pak můžeme nerovnici upravit do podoby x ≤ bj aj a jedná se o polopřímku omezenou zprava. • Jestliže aj < 0, pak můžeme nerovnici upravit do podoby x ≥ bj aj a jedná se o polopřímku omezenou zleva. • Jestliže aj = 0, pak v případě, že bj ≥ 0, je j-tá nerovnice splněna pro každé x a může být vynechána, v opačném případě, tj. je-li bj < 0, není j-tá nerovnice nikdy splněna a úloha tak nemá přípustné řešení. Pokud pro žádnou z nerovnic nenastane případ aj = 0 a bj < 0, uvažme následující dvě množiny indexů. Jl = {j ∈ N | 1 ≤ j ≤ n ∧ aj < 0} , Jr = {j ∈ N | 1 ≤ j ≤ n ∧ aj > 0} . 11 KAPITOLA 2. JEDNODIMENZIONÁLNÍ ÚLOHA Označíme-li xl = max bj aj j ∈ Jl , jestliže Jl = ∅, −∞, jestliže Jl = ∅, xr = min bj aj j ∈ Jr , jestliže Jr = ∅, ∞, jestliže Jr = ∅, pak množina všech přípustných řešení jednodimenzionální úlohy lineárního programování je ekvivalentní vztahům x ≥ xl a x ≤ xr. Ze všech polopřímek omezených zleva (respektive zprava) jsme vybrali tu, jejíž koncový bod je největší (resp. nejmenší), neboť je podmnožinou všech polopřímek omezených zleva (resp. zprava). Uvažme nejprve variantu, že obě uvedené množiny indexů jsou neprázdné, tedy xl = −∞ a xr = ∞. Pokud xl > xr, pak úloha nemá přípustné řešení. V opačném případě je množinou všech přípustných řešení úsečka s krajními body xl a xr. Vzhledem k linearitě účelové funkce je optimálním řešením buď bod xl, anebo bod xr. Stačí tedy porovnat hodnoty účelových funkcí v těchto bodech. Jestliže je některá z množin Jl nebo Jr prázdná, může být úloha neomezená a nemusí tedy existovat optimální řešení. Jednotlivé možnosti závisí na tom, zda je účelová funkce rostoucí, tj. c > 0, anebo klesající, tj. c < 0. Mohou nastat následující případy. • Jl = ∅ a Jr = ∅, tedy xl = −∞ a množina všech přípustných řešení je polopřímka omezená zprava. Je-li účelová funkce rostoucí, je optimálním řešením bod xr. Pokud je účelová funkce klesající, úloha je neomezená. • Jl = ∅ a Jr = ∅, tedy xr = ∞ a množina všech přípustných řešení je polopřímka omezená zleva. Jestliže je účelová funkce klesající, pak optimálním řešením je bod xl. V opačném případě je úloha neomezená. • Jl = ∅ a Jr = ∅, tedy xl = −∞, xr = ∞ a množina všech přípustných řešení je celá přímka. Úloha je neomezená v případě rostoucí i klesající účelové funkce. Uvedené poznatky o jednodimenzionální úloze lineárního programování můžeme shrnout v algoritmu 1. S tímto algoritmem se v různých obměnách setkáme při následné analýze úlohy lineárního programování v rovině. 12 KAPITOLA 2. JEDNODIMENZIONÁLNÍ ÚLOHA Algoritmus 1 Algoritmus k řešení jednodimenzionální úlohy OneDimLP(c, a1, . . . , an, b1, . . . , bn) Vstup: Reálná čísla c, a1, . . . , an, b1, . . . , bn, kde f(x) = cx je účelová funkce úlohy s množinou všech přípustných řešení danou nerovnicemi ajx ≤ bj, j ∈ {1, . . . , n}. Výstup: Optimální řešení jednodimenzionální úlohy, případně zpráva o tom, že úloha nemá přípustné řešení či je neomezená. 1: xl ← −∞; xr ← ∞; 2: for i = 1 to n do 3: if ai > 0 then 4: xr ← min(xr, bi ai ); 5: else if ai < 0 then 6: xl ← max(xl, bi ai ); 7: else 8: if bi < 0 then 9: return Úloha nemá přípustné řešení. 10: if xl = −∞ ∧ xr = ∞ then 11: if xl > xr then 12: return Úloha nemá přípustné řešení. 13: else 14: if c > 0 then 15: return xr 16: else 17: return xl 18: else 19: if xl = −∞ ∧ xr = ∞ ∧ c > 0 then 20: return xr 21: else if xl = −∞ ∧ xr = ∞ ∧ c < 0 then 22: return xl 23: else 24: return Úloha je neomezená. 13 Kapitola 3 Algoritmus pro omezenou úlohu lineárního programování Vraťme se nyní k úloze lineárního programování v rovině. V této části se budeme zabývat omezenou úlohou. Algoritmus, který zde popíšeme, je založen na postupném přidávání lineárních omezení (polorovin). V jednotlivých krocích jsou řešeny úlohy lineárního programování, jejichž množiny přípustných řešení jsou dány prvními i polorovinami. Účelová funkce všech těchto úloh je stejná. Ze známého optimálního řešení úlohy dané prvními i − 1 polorovinami algoritmus odvodí optimální řešení úlohy, která vznikne přidáním další poloroviny, tj. poloroviny hi. Avšak algoritmus vyžaduje, aby každá z postupně řešených úloh měla právě jedno optimální řešení, s výjimkou případu, kdy je množina přípustných řešení prázdná. Toho docílíme následujícími dvěma úmluvami. Protože chceme docílit toho, aby v každém kroku úloha nebyla neomezená, přidáme k úloze další dvě lineární omezení. Na jejich základě získáme počáteční řešení. Pokud je původní úloha omezená, existuje dostatečně velká konstanta M taková, že přidáním podmínek |x| ≤ M a |y| ≤ M nebude ovlivněno optimální řešení. Tento zápis nám však dává čtyři lineární omezení, při následující specifikaci v závislosti na vektoru c stačí přidat pouze dvě, a to m1 a m2. cx > 0 cx ≤ 0 cy > 0 m1 : x ≤ M m1 : x ≥ −M m2 : y ≤ M m2 : y ≤ M cy ≤ 0 m1 : x ≤ M m1 : x ≥ −M m2 : y ≥ −M m2 : y ≥ −M Viděli jsme, že může nastat situace, kdy má úloha lineárního programování více optimálních řešení (obrázek 1.2(d)). V takovém případě budeme hledat lexikograficky nejmenší řešení. Je ale možné, že úloha, která má více optimálních řešení, nemá lexikograficky nejmenší řešení. Touto situací se budeme zabývat v následující kapitole. V této části proto v případě více optimálních řešení předpokládejme existenci lexikograficky nejmenšího řešení. 14 KAPITOLA 3. ALGORITMUS PRO OMEZENOU ÚLOHU LINEÁRNÍHO PROGRAMOVÁNÍ Definice 3.1. Nechť a, b ∈ R2 , a = [ax, ay], b = [bx, by]. Řekneme, že bod a je lexikograficky menší než bod b, jestliže platí: ax < bx ∨ (ax = bx ∧ ay < by). Uvažujme úlohu lineárního programování (H, c), kde H = {h1, h2, . . . , hn} a m1, m2 jsou výše uvedená dodatečná lineární omezení. Zaveďme následující označení: • Hi = {m1, m2, h1, h2, . . . hi} pro každé i ∈ {0, . . . , n}, • Ci = Hi = m1 ∩ m2 ∩ h1 ∩ · · · ∩ hi pro každé i ∈ {0, . . . , n}, tedy Ci je množina všech přípustných řešení úlohy (Hi, c), • vi je optimální (lexikograficky nejmenší) řešení úlohy (Hi, c) pro každé i ∈ {0, . . . , n}, • li je hraniční přímka poloroviny hi pro každé i ∈ {1, . . . , n}. Pokud je úloha (H, c) omezená a m1, m2 vhodně zvolená dodatečná lineární omezení, tj. neovlivňující optimální řešení, pak jsou úlohy (Hn, c) a (H, c) ekvivalentní. Jestliže pro nějaké j ∈ {0, . . . , n} nemá úloha (Hj, c) přípustné řešení, tj. je-li Cj = ∅, pak nemá přípustné řešení ani žádná z úloh (Hi, c) pro i ∈ {j + 1, . . . , n}, neboť dle definice Ci zřejmě platí C0 ⊇ C1 ⊇ · · · ⊇ Cn. Speciálně pak nemá přípustné řešení ani úloha (H, c) v případě, že je omezená. Věta 3.1. Pro každé i ∈ {1, . . . , n} platí: (i ) jestliže vi−1 ∈ hi, pak vi = vi−1, (ii ) pokud vi−1 /∈ hi, pak buď vi ∈ li, nebo Ci = ∅. Důkaz. (i) Nechť vi−1 ∈ hi. Víme, že vi−1 ∈ Ci−1 a Ci = Ci−1 ∩ hi. Proto také vi−1 ∈ Ci. Jelikož Ci ⊆ Ci−1 a vi−1 je optimálním řešením na Ci−1, je vi−1 optimálním řešením také na Ci, a tedy vi = vi−1. (ii) Důkaz vedeme sporem: předpokládejme, že vi−1 /∈ hi, Ci = ∅ a vi /∈ li. Protože Ci je podmnožinou Ci−1, je vi prvkem Ci−1. Jelikož také vi−1 ∈ Ci−1 a současně Ci−1 je konvexní oblastí, leží v Ci−1 celá úsečka s krajními body vi a vi−1. Podle předpokladu leží vi−1 v polorovině opačné k hi. Naopak zřejmě vi leží v polorovině hi, neboť vi ∈ Ci a Ci = Ci−1 ∩ hi. Z toho vyplývá, že úsečka vivi−1 protíná li. Označme tento průsečík q. Protože q leží v polorovině hi, je také prvkem množiny Ci. Jelikož vi−1 je optimální řešení na Ci−1, musí být hodnota účelové funkce v bodě vi−1 větší, anebo rovna hodnotě účelové funkce v bodě vi. Protože účelová funkce je lineární na R2 , je také lineární, případně konstantní, na každé úsečce v R2 . Nejprve předpokládejme, že f(vi−1) > f(vi). Pak z linearity účelové funkce na úsečce vivi−1 dostáváme vztah f(vi−1) > f(q) > f(vi). Z druhé nerovnosti však dostáváme spor s optimalitou vi na Ci. 15 KAPITOLA 3. ALGORITMUS PRO OMEZENOU ÚLOHU LINEÁRNÍHO PROGRAMOVÁNÍ Ci−1 hi vi vi−1 q Ci li V případě, že f(vi−1) = f(vi), je účelová funkce na úsečce vivi−1 konstantní, a tedy f(vi−1) = f(q) = f(vi). Víme, že vi−1 je lexikograficky nejmenší optimální řešení na Ci−1, speciálně vi−1 je lexikograficky menší než vi. Proto je také q lexikograficky menší než vi, což je spor. Věta 3.1 nám dává částečný návod k nalezení optimálního řešení úlohy Hi. V případě (i) jsou optimální řešení úloh Hi a Hi−1 totožná. Nastane-li však případ (ii), věta 3.1 nám dává pouze nápovědu, kde příslušné řešení hledat. Při hledání se můžeme omezit pouze na přímku li. Předpokládejme nyní, že vi−1 /∈ hi. Hledané optimální řešení úlohy Hi musí být současně i přípustným řešením, tj. vi ∈ Ci, a proto musí ležet ve všech polorovinách1 h1, . . . , hi−1. Z uvedených omezení dostáváme novou úlohu lineárního programování v rovině. Označme ji K. Účelová funkce je stejná jako u původní úlohy, množina všech přípustných řešení je dána vztahy: ai1x + ai2y = bi, (3.1) a11x + a12y ≤ b1, a21x + a22y ≤ b2, ... ai−1,1x + ai−1,2y ≤ bi−1. Jedná se však o speciální případ, který můžeme snadno vyřešit, neboť úlohu lze převést na jednodimenzionální úlohu lineárního programování. Uvažujme nyní o průniku přímky li a poloroviny hj pro 1 ≤ j ≤ i − 1. Upravujme tedy nerovnici aj1x + aj2y ≤ bj za předpokladu, že platí ai1x + ai2y = bi. Rozlišme nyní dva případy. Předpokládejme, že ai2 = 0. Pak vyjádřením y z rovnice (3.1) dostáváme y = bi − ai1x ai2 . 1 V polorovině hi leží automaticky, neboť dle věty 3.1 leží na její hraniční přímce. 16 KAPITOLA 3. ALGORITMUS PRO OMEZENOU ÚLOHU LINEÁRNÍHO PROGRAMOVÁNÍ Dosaďme za y do nerovnice poloroviny hj a upravujme. aj1x + aj2y ≤ bj aj1x + aj2 bi − ai1x ai2 ≤ bj aj1ai2x + aj2bi − aj2ai1x ai2 ≤ bj x (aj1ai2 − aj2ai1) + aj2bi ai2 ≤ bj aj1ai2 − aj2ai1 ai2 x ≤ bj − aj2bi ai2 (3.2) Označme Dx j = aj1ai2 − aj2ai1 ai2 a Ex j = bj − aj2bi ai2 . Pak v závislosti na Dx j dostáváme následující omezení pro x: • Dx j > 0: x ≤ Ex j Dx j , • Dx j < 0: x ≥ Ex j Dx j , • Dx j = 0: Buď Ex j ≥ 0, a tedy nedostáváme žádné omezení pro x (platí li ⊂ hj), anebo Ex j < 0 (tj. li ∩ hj = ∅), pak úloha K nemá přípustné řešení a podle věty 3.1 je Ci prázdnou množinou. Předpokládejme, že ai2 = 0. Postupujeme analogicky jako v předchozím případě. Z rovnice (3.1) vyjádříme x, dosadíme do nerovnice poloroviny hj a po úpravách dostaneme aj2y ≤ bj − aj1bi ai1 . (3.3) Označme Dy j = aj2 a Ey j = bj − aj1bi ai1 . V závislosti na Dy j dostáváme následující omezení pro y: • Dy j > 0: y ≤ Ey j Dy j , • Dy j < 0: y ≥ Ey j Dy j , • Dy j = 0: Buď Ey j ≥ 0, a tedy nedostáváme žádné omezení pro y (platí li ⊂ hj), anebo Ey j < 0 (tj. li ∩ hj = ∅), pak úloha K nemá přípustné řešení a podle věty 3.1 je Ci prázdnou množinou. 17 KAPITOLA 3. ALGORITMUS PRO OMEZENOU ÚLOHU LINEÁRNÍHO PROGRAMOVÁNÍ Pro 1 ≤ j ≤ i − 1 zaveďme následující označení. Dj = Dx j , jestliže ai2 = 0, Dy j , jestliže ai2 = 0, Ej = Ex j , jestliže ai2 = 0, Ey j , jestliže ai2 = 0. Uvažme nyní výše uvedená omezení pro všechna j ∈ {1, . . . , i−1}. Pokud pro nějaké takové j platí Dj = 0 a Ej < 0, pak Ci = ∅. V opačném případě můžeme tato lineární omezení rozdělit do dvou množin na omezení shora a omezení zdola. Zaveďme tedy označení následujících množin indexů. Jl = {j ∈ N | 1 ≤ j ≤ i − 1 ∧ Dj < 0} , Jr = {j ∈ N | 1 ≤ j ≤ i − 1 ∧ Dj > 0} . Úloha K je tedy ekvivalentní jednodimenzionální úloze K s proměnnou z, jejíž množina všech přípustných řešení je dána vztahy z ≤ Ej Dj , j ∈ Jr, z ≥ Ej Dj , j ∈ Jl a jejíž účelová funkce je f(z) =    f z, bi−ai1z ai2 , jestliže ai2 = 0, f bi ai1 , z , jestliže ai2 = 0. Jestliže z je optimální řešení úlohy K, pak pro optimální řešení úlohy K, tedy vi, platí následující vztah: vi =    z, bi−ai1z ai2 , jestliže ai2 = 0, bi ai1 , z , jestliže ai2 = 0. Získali jsme tak jednodimenzionální úlohu lineárního programování, kterou již umíme řešit na základě kapitoly 2 tohoto textu. V případě, že Jl = ∅ a Jr = ∅, označíme zl = max Ej Dj j ∈ Jl a zr = min Ej Dj j ∈ Jr , množina všech přípustných řešení úlohy K je pak ekvivalentně dána vztahy z ≥ zl, z ≤ zr. Jestliže zl > zr, pak úloha K nemá přípustné řešení, a tedy podle věty 3.1 je Ci = ∅. Pokud zl ≤ zr, je množinou všech přípustných řešení úlohy K úsečka zlzr. Vzhledem 18 KAPITOLA 3. ALGORITMUS PRO OMEZENOU ÚLOHU LINEÁRNÍHO PROGRAMOVÁNÍ k linearitě účelové funkce je optimálním řešením jeden z krajních bodů této úsečky. Je jím bod zl, pokud f(zl) ≥ f(zr) (v případě rovnosti je zl lexikograficky menší). Jestliže f(zl) < f(zr), optimálním řešením úlohy K je bod zr. Jestliže však Jl = ∅ nebo Jr = ∅, na rozdíl od obecného případu jednodimenzionální úlohy analyzovaného v kapitole 2, nemůže nastat situace, kdy je úloha neomezená. To je zaručeno dodatečnými lineárními omezeními, která byla volena tak, že účelová funkce je shora omezená na množině C0 = m1 ∩m2, a tedy i na jejích podmnožinách C1, . . . , Cn. Na základě dosavadních poznatků formulujme algoritmus 2 pro řešení omezené úlohy lineárního programování. Věta 3.2. Časová složitost algoritmu 2 je O(n2 ). Důkaz. Časová složitost každé iterace cyklu for (ř. 2–35) závisí na vyhodnocení podmínky na řádku 3. První větev (ř. 4) má konstantní časovou složitost. Druhá část větvení (ř. 6–35) obsahuje 3 for cykly (ř. 7–8, 10–11, 13–20), z nichž každý má časovou složitost O(i). Ostatní části této větve mají konstantní časovou složitost. Celkem je tedy časová složitost i-té iterace cyklu for (ř. 2–35) v závislosti na splnění podmínky na řádku 3 buď konstantní, anebo O(i). V nejhorším případě může v každé iteraci dojít k vykonání druhé větve, a proto je celková složitost algoritmu 2 rovna O(1 + 2 + · · · + n) = O(n(n+1) 2 ) = O(n2 ). 19 KAPITOLA 3. ALGORITMUS PRO OMEZENOU ÚLOHU LINEÁRNÍHO PROGRAMOVÁNÍ Algoritmus 2 Algoritmus pro omezenou úlohu lineárního programování BoundedLP(H, c, m1, m2) Vstup: Úloha lineárního programování (H, c) a dodatečná lineární omezení m1, m2 taková, že úloha ({m1, m2}, c) je omezená. Výstup: Optimální řešení úlohy (H ∪ {m1, m2}, c), případně oznámení, že úloha nemá přípustné řešení. 1: v0 ← průnik hraničních přímek polorovin m1 a m2; 2: for i = 1 to n do 3: if vi−1 ∈ hi then 4: vi ← vi−1; 5: else 6: if ai2 = 0 then 7: for j = 1 to i − 1 do 8: Dj ← aj1ai2−aj2ai1 ai2 ; Ej ← bj − aj2bi ai2 ; 9: else 10: for j = 1 to i − 1 do 11: Dj ← aj2; Ej ← bj − aj1bi ai1 ; 12: zl ← −∞; zr ← ∞; 13: for j = 1 to i − 1 do 14: if Dj > 0 then 15: zr ← min(zr, Ej Dj ); 16: else if Dj < 0 then 17: zl ← max(zl, Ej Dj ); 18: else 19: if Ej < 0 then 20: return Nemá přípustné řešení 21: if zl = −∞ then 22: z ← zr; 23: else if zr = ∞ then 24: z ← zl; 25: else if zl > zr then 26: return Nemá přípustné řešení 27: else 28: if f(zl) ≥ f(zr) then 29: z ← zl; 30: else 31: z ← zr; 32: if ai2 = 0 then 33: vi ← z, bi−ai1z ai2 ; 34: else 35: vi ← bi ai1 , z ; 36: return vn 20 Kapitola 4 Neomezená úloha lineárního programování Doposud jsme se zabývali omezenou úlohou lineárního programování. Můžeme se však setkat i s úlohami, které jsou neomezené, jak jsme např. viděli na obrázku 1.2(d). V takovém případě nám algoritmus 2 nevrátí správný výsledek. Díky dodatečně přidaným lineárním omezením m1, m2 bude výstupem nějaký bod v rovině, přestože neomezená úloha nemá optimální řešení. Pokusíme se tedy odvodit algoritmus, pomocí kterého dokážeme rozpoznat, zda je daná úloha lineárního programování neomezená. Úloha lineárního programování je neomezená právě tehdy, když existuje polopřímka ρ, která je podmnožinou množiny všech přípustných řešení a podél níž účelová funkce neomezeně roste. Nechť se jedná o polopřímku ρ = {p + λs | λ ∈ R; λ ≥ 0} , kde p = [px, py] je její počáteční bod a s = (sx, sy) její směrový vektor. Pro každé i označme ni normálový vektor přímky li, který směřuje do poloroviny hi. Věta 4.1. Úloha lineárního programování (H, c) je neomezená právě tehdy, když existuje vektor s splňující s, c > 0, s, ni ≥ 0 pro každé i ∈ {1, · · · , n} a současně H = ∅, kde H = {hi ∈ H | s, ni = 0}. Důkaz. Nechť úloha (H, c) je neomezená. Pak existuje polopřímka ρ s výše uvedenými vlastnostmi a se směrovým vektorem s. Účelová funkce podél polopřímky ρ neomezeně roste, proto musí vektory s a c svírat úhel v intervalu 0, π 2 . Protože pro úhel ϕ svíraný vektory s a c platí cos ϕ = s, c s · c , z vlastností funkce cos dostáváme s, c > 0. Protože polopřímka ρ leží v množině všech přípustných řešení úlohy (H, c), musí vektor s svírat s každým z normálových vektorů ni úhel v intervalu 0, π 2 . Analogicky jako výše tedy dostáváme s, ni ≥ 0. 21 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Úloha (H, c) je neomezená, a proto je její množina přípustných řešení neprázdná, tj. H = ∅. Protože H ⊆ H, je také H = ∅. Vektor s je tedy hledaným vektorem s požadovanými vlastnostmi. Nechť nyní existuje vektor s, pro který platí s, c > 0, s, ni ≥ 0, i ∈ {1, · · · , n} a H = ∅. Z posledního vztahu vyplývá, že existuje bod p0 ∈ H . Uvažme polopřímku ρ0 = {p0 + λs | λ ∈ R; λ ≥ 0}. Pro každou polorovinu hi z H je normálový vektor ni její hraniční přímky li kolmý na vektor s (neboť s, ni = 0), a tedy polopřímka ρ0 je s přímkou li rovnoběžná. Proto leží celá polopřímka ρ0 v každé z polorovin hi z H . Pro každou polorovinu hi z H \ H je s, ni > 0, a tedy vektory s a ni svírají úhel z intervalu 0, π 2 . Z toho vyplývá, že pro každou polorovinu hi z H \ H existuje λi ∈ R, λi ≥ 0 takové, že (p0 + λs) ∈ hi pro každé λ ≥ λi. Označíme-li λmax = max{λi | hi ∈ H \ H }, pak (p0 + λs) ∈ H pro každé λ ≥ λmax. Zvolme p = p0 + λmaxs. Pak polopřímka ρ = {p + λs | λ ∈ R; λ ≥ 0} je podmnožinou množiny všech přípustných řešení úlohy (H, c) a zároveň účelová funkce podél ní neomezeně roste, neboť s, c > 0. Úloha (H, c) je tedy neomezená. Podaří-li se nám najít vektor splňující podmínky věty 4.1, je úloha neomezená. V takovém případě budeme požadovat, aby byla výstupem algoritmu polopřímka svědčící o neomezenosti úlohy, tj. polopřímka, která leží v množině všech přípustných řešení a podél níž účelová funkce neomezeně roste. Tuto polopřímku můžeme sestrojit postupem uvedeným v druhé části důkazu věty 4.1. Pokud dokážeme, že neexistuje vektor splňující podmínky věty 4.1, jedná se o omezenou úlohu lineárního programování. Podívejme se nyní, jak vypadají jednotlivé složky vektoru ni, který jsme zadefinovali jako normálový vektor hraniční přímky li poloroviny hi směřující do poloroviny hi. Víme, že přímka li je dána rovnicí ai1x + ai2y = bi. Protože se jedná o obecnou rovnici přímky, je vektor (ai1, ai2) normálovým vektorem přímky li. Normálovým vektorem je však také vektor (−ai1, −ai2). Jelikož nám nezáleží na velikosti vektoru, ale pouze na směru, stačí nám zjistit, který z těchto vektorů směřuje do poloroviny hi. Normálový vektor hraniční přímky nezávisí na posunutí poloroviny. Proto můžeme místo poloroviny hi = {(x, y) ∈ R2 | ai1x + ai2y ≤ bi} uvažovat polorovinu hi = {(x, y) ∈ R2 | ai1x + ai2y ≤ 0}. Normálový vektor zůstane stejný, navíc směřuje do poloroviny hi právě tehdy, když směřuje do poloroviny hi. Hraniční přímka li poloroviny hi prochází počátkem soustavy souřadnic. Umístíme-li počáteční bod vektoru (ai1, ai2), respektive (−ai1, −ai2), do počátku soustavy souřadnic, je koncovým bodem tohoto vektoru bod [ai1, ai2], respektive [−ai1, −ai2]. Přičemž vektor vycházející z počátku soustavy souřadnic směřuje do poloroviny hi právě tehdy, když v ní leží jeho koncový bod. Dosazením obou koncových bodů do 22 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ nerovnice poloroviny hi získáváme: a2 i1 + a2 i2 ≤ 0, −(a2 i1 + a2 i2) ≤ 0. První nerovnost není splněna nikdy1 , druhá platí vždy. Proto můžeme za normálový vektor ni směřující do poloroviny hi považovat vektor (−ai1, −ai2). Na základě věty 4.1 hledejme vektor s = (sx, sy) splňující s, c > 0, s, ni ≥ 0, i ∈ {1, . . . , n}. (4.1) S využitím výše odvozených normálových vektorů ni = (−ai1, −ai2) a rozepsáním skalárních součinů dostáváme ekvivalentní vztahy cxsx + cysy > 0, (4.2) −ai1sx − ai2sy ≥ 0, i ∈ {1, · · · , n}. (4.3) Zřejmě pokud vektor s splňuje tyto vztahy, pak je splňuje také vektor λs pro libovolné λ ∈ R, λ > 0. Jelikož nám nezáleží na velikosti vektorů, ale pouze na jejich směru, můžeme nerovnici (4.2) nahradit rovnicí cxsx + cysy = 1. (4.4) Tím nijak neomezíme možný směr vektoru s, omezíme pouze jeho velikost, a to tak, že požadujeme, aby jeho koncový bod ležel na přímce dané rovnicí (4.4) (počáteční bod leží v počátku soustavy souřadnic). Ke každému vektoru s splňujícímu nerovnice (4.2) a (4.3) existuje vhodné λ takové, že vektor λs splňuje rovnici (4.4) a nerovnice (4.3). Konkrétně při volbě λ = 1 cxsx + cysy dostáváme pro vektor u = λs cxux + cyuy = cx sx cxsx + cysy + cy sy cxsx + cysy = cxsx + cysy cxsx + cysy = 1. Vektor u tedy splňuje rovnici (4.4). Navíc z platnosti nerovnic (4.3) pro vektor s získáme vynásobením každé z nerovnic číslem λ platnost týž nerovnic pro vektor u. Naopak každý vektor s splňující rovnici (4.4) automaticky splňuje nerovnici (4.2). Nahrazení nerovnice (4.2) rovnicí (4.4) tedy neznamená žádné nové omezení pro směr vektoru s. Předpokládejme nyní, že cy = 0. Pak vyjádřením sy z rovnice (4.4) získáváme vztah sy = 1 cy − cx cy sx. (4.5) 1 Na začátku jsme požadovali, aby alespoň jedno z čísel ai1, ai2 bylo nenulové. 23 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Dosaďme za sy do každé z nerovnic (4.3) a upravujme. −ai1sx − ai2 cy + ai2cx cy sx ≥ 0 ai2cx cy − ai1 sx ≥ ai2 cy Označme Fx i = ai2cx cy − ai1, Gx i = ai2 cy . Pokud existuje j ∈ {1, . . . , n} takové, že Fx j = 0 a současně Gx j > 0, pak neexistuje vektor s, který by splňoval vztahy (4.1), a úloha (H, c) je podle věty 4.1 omezená. V opačném případě dostáváme následující omezení pro sx. sx ≥ Gx i Fx i , i ∈ Il, sx ≤ Gx i Fx i , i ∈ Ir, (4.6) kde Il a Ir jsou následující množiny indexů: Il = {i ∈ N | 1 ≤ i ≤ n ∧ Fx i > 0} , Ir = {i ∈ N | 1 ≤ i ≤ n ∧ Fx i < 0} . Označíme-li xl = max Gx i Fx i i ∈ Il , jestliže Il = ∅, −∞, jestliže Il = ∅, xr = min Gx i Fx i i ∈ Ir , jestliže Ir = ∅, ∞, jestliže Ir = ∅, pak všechny nerovnice (4.6) můžeme nahradit dvojicí jim ekvivalentních vztahů sx ≥ xl, sx ≤ xr. Pokud xr < xl, neexistuje vektor s, který by splňoval vztahy (4.1), a úloha (H, c) je podle věty 4.1 omezená. V opačném případě získáme vektor s vyhovující vztahům (4.1) volbou libovolného sx splňujícího xl ≤ sx ≤ xr a dopočtením sy ze vztahu (4.5). Např. volbou sx = xl (v případě, že xl = −∞) dostáváme vektor s = xl, 1 cy − cx cy xl . Za předpokladu, že cy = 0 dostáváme vyjádřením sx z rovnice (4.4) vztah sx = 1 cx . Jeho dosazením do každé z nerovnic (4.3) po úpravě získáme −ai2sy ≥ ai1 cx . 24 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Označíme-li Fy i = −ai2, Gy i = ai1 cx , můžeme dále postupovat zcela analogicky jako v případě, že cy = 0. Buď zjistíme, že vektor s s požadovanými vlastnostmi neexistuje a úloha (H, c) je podle věty 4.1 omezená, anebo získáme vektor s vyhovující vztahům (4.1). Jestliže jsme získali vektor s splňující vztahy (4.1), stále se ještě nemusí jednat o neomezenou úlohu lineárního programování. Dle věty 4.1 musí být splněna ještě jedna podmínka, a to H = ∅. Pokud tato podmínka splněna není, tj. je-li H = ∅, je také H = ∅, neboť H ⊆ H, a tedy úloha (H, c) nemá přípustné řešení. Označme I množinu indexů I = {i ∈ N | hi ∈ H }. Průnik všech polorovin z H je pak řešením nerovnic ai1x + ai2y ≤ bi, i ∈ I . Podle definice H je s, ni = −ai1sx − ai2sy = 0, i ∈ I . Normálový vektor hraniční přímky li libovolné poloroviny hi z H je tedy kolmý na vektor s. To znamená, že hraniční přímky všech polorovin z H jsou rovnoběžné s vektorem s. V případě, že sy = 0, tj. pokud není s rovnoběžný s osou x, můžeme místo polorovin hi, i ∈ I uvažovat jejich průniky s osou x, což jsou polopřímky určené nerovnicemi ai1x ≤ bi, i ∈ I . (4.7) Vzhledem k tomu, že hraniční přímky všech polorovin hi, i ∈ I jsou rovnoběžné, mají všechny tyto polopřímky neprázdný průnik právě tehdy, když mají neprázdný průnik všechny poloroviny hi, i ∈ I . Situace je opět zcela analogická té, se kterou jsme se již v tomto textu několikrát setkali. Hledáme však pouze množinu všech přípustných řešení jednodimenzionální úlohy lineárního programování, aniž by nás zajímala její účelová funkce či optimální řešení. Nerovnice (4.7) můžeme nahradit dvojicí nerovnic x ≥ xl, x ≤ xr, kde xl = max bi ai1 i ∈ I ∧ ai1 < 0 ∪ {−∞} , xr = min bi ai1 i ∈ I ∧ ai1 > 0 ∪ {∞} . Jestliže xl ≤ xr, je H = ∅ a podle věty 4.1 je úloha (H, c) neomezená. V opačném případě je H = ∅, proto také H = ∅, a tedy úloha (H, c) nemá přípustné řešení. V případě, že sy = 0, tj. s je rovnoběžný s osou x, postupujeme analogicky, s tím rozdílem, že uvažujeme průniky s osou y, což jsou polopřímky určené nerovnicemi ai2y ≤ bi, i ∈ I . 25 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Uvažujme nyní o omezené úloze lineárního programování. Vraťme se k situaci, kde jsme zjistili, že hledaný vektor s neexistuje a úloha je tedy omezená. Předpokládejme, že cy = 0, v opačném případě bychom postupovali zcela analogicky. Omezenost úlohy mohla být prokázána ve dvou odlišných situacích. První možností je, že omezenost úlohy jsme zjistili na základě platnosti vztahu xl > xr. Za této situace2 označme • hl takovou polorovinu hj, že výraz Gx i Fx i , i ∈ Il nabývá maxima pro i = j, • hr takovou polorovinu hj, že výraz Gx i Fx i , i ∈ Ir nabývá minima pro i = j. Uvážíme-li úlohu ({hl, hr}, c), získáme pro ni stejné závěry o neexistenci vhodného vektoru s jako pro úlohu (H, c). Podle věty 4.1 je tedy i úloha ({hl, hr}, c) omezená. Poloroviny hl a hr proto představují svědky omezenosti úlohy (H, c), které můžeme použít místo dodatečných lineárních omezení v algoritmu pro omezenou úlohu lineárního programování. Situace je znázorněna na obrázku 4.1(a). c hl hr (a) Poloroviny hl a hr jsou svědky c hj (b) Jediná polorovina hj je svědkem Obrázek 4.1: Svědci omezenosti úlohy (H, c) Druhý možný případ, kdy jsme mohli zjistit, že úloha je omezená, nastal, pokud existuje j, pro které je Fx j = 0 a Gx j > 0. Pak vyjádřením Fx j a Gx j po úpravě dostáváme: aj2cx − aj1cy = c, (aj2, −aj1) = 0, (4.8) aj2 cy > 0. (4.9) Víme, že (aj1, aj2) je normálový vektor přímky lj, a proto (aj2, −aj1) je její směrový vektor. Ze vztahu (4.8) tedy vyplývá, že vektor c je kolmý na přímku lj. Navíc na základě vztahu (4.9) platí, že cy a −aj2 mají opačná znaménka. Jelikož (−aj1, −aj2) je normálový vektor přímky lj směřující do poloroviny hj, vektor c nesměřuje do poloroviny hj. Nastává situace znázorněná na obrázku 4.1(b). 2 Tento případ mohl nastat pouze pokud xl = −∞ a xr = ∞. 26 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Z výše uvedeného postupu vyplývá, že ani pro úlohu lineárního programování ({hj}, c) neexistuje vektor s požadovaných vlastností, a proto je podle věty 4.1 i úloha ({hj}, c) omezená. Polorovina hj je tedy svědkem toho, že úloha (H, c) je omezená. Situace je však odlišná od té, kdy jsme měli k dispozici dva svědky. Může se stát, že úloha je sice omezená, ale nemá lexikograficky nejmenší řešení. Jelikož přímka lj je kolmá na vektor c, je hodnota účelové funkce ve všech bodech ležících na přímce lj stejná. Body na přímce lj však můžeme lexikograficky uspořádat. Pokud se nám podáří nalézt takovou polorovinu hk, že úloha ({hj, hk}, c) má lexikograficky nejmenší řešení, pak v případě, že množina všech přípustných řešení původní úlohy (H, c) je neprázdná, musí mít lexikograficky nejmenší řešení i původní úloha. V takovém případě můžeme místo dodatečných omezení m1, m2 v algoritmu 2 použít poloroviny hj, hk a za počáteční řešení v0 považovat průnik jejich hraničních přímek. Situace je znázorněna na obrázku 4.2(a). Ukážeme-li, že neexistuje polorovina hk ∈ H \ {hj} výše uvedených vlastností, znamená to, že úloha nemá lexikograficky nejmenší řešení, případně nemá přípustné řešení. Označme H množinu všech polorovin z H, jejichž hraniční přímky jsou rovnoběžné s přímkou lj (a tedy kolmé na vektor c). Jestliže je H = ∅, pak úloha H nemá přípustné řešení, neboť H ⊆ H. Pokud však H = ∅, úloha H má nekonečně mnoho přípustných řešení a nemá lexikograficky nejmenší optimální řešení. Příklad takové úlohy je znázorněn na obrázku 4.2(b). c hj hk v0 (a) Existuje lexikograficky nejmenší řešení c hj (b) Neexistuje lexikograficky nejmenší řešení Obrázek 4.2: Existence lexikograficky nejmenšího řešení v případě poloroviny s hraniční přímkou kolmou na vektor c Hledejme tedy polorovinu hk, jejíž hraniční přímka protíná přímku lj a současně platí, že úloha ({hj, hk}, c) má lexikograficky nejmenší řešení. Uvažme průniky přímky lj s polorovinami z H \ {hj}. V případě, že přímka lj není kolmá na osu x, můžeme nám známým způsobem řešit jednodimenzionální úlohu lineárního programování pro x-ové souřadnice těchto průniků. V opačném případě řešíme analogickou úlohu pro y-ové souřadnice. Postup je téměř stejný jako v případě odvození algoritmu pro omezenou úlohu lineárního programování, vycházíme proto ze vztahů (3.2) a (3.3). Nemusíme však řešit celou úlohu. Zajímá nás pouze, jestli 27 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ některý z průniků má podobu polopřímky omezené zleva, neboť hledáme lexikograficky nejmenší řešení. Jakmile najdeme polorovinu hk, jejíž průnik s přímkou lj je polopřímka omezená zleva, můžeme přejít k algoritmu 2. Pokud žádný z průniků nemá podobu polopřímky omezené zleva, ověříme, zda H = ∅. To je stejný problém jako zjištění, zda je H = ∅, s tím rozdílem, že místo vektoru s použijeme vektor (cy, −cx), který je kolmý na vektor c. Oba problémy můžeme řešit pomocí algoritmu 5. Uvedený postup ke zjištění existence lexikograficky nejmenšího řešení je využit v algoritmu 4. Pokud úloha má přípustné řešení, ale nemá lexikograficky nejmenší optimální řešení, budeme požadovat jako výstup lexikograficky největší optimální řešení. Pokud ani to neexistuje, výstupem bude nějaké optimální řešení. To získáme na základě algoritmu 5. Můžeme použít velmi podobný postup jako v předchozím odstavci. Hledáme tedy polorovinu hk, jejíž průnik s přímkou lj je polopřímka omezená zprava. K tomu můžeme použít algoritmus 4, který mírně modifikujeme, a to tak, že nerovnosti na řádcích 3 a 7 zaměníme za opačné. V případě, že nalezneme polorovinu hk, jejíž hraniční přímka protíná přímku lj a současně platí, že úloha ({hj, hk}, c) má lexikograficky největší řešení, můžeme využít upravený algoritmus 2 k řešení této úlohy. Záměnou neostré nerovnosti na řádku 28 algoritmu 2 za ostrou zajistíme, aby algoritmus hledal lexikograficky největší optimální řešení místo nejmenšího. Jelikož se jedná jen o drobné úpravy, nebudeme zde upravené algoritmy explicitně uvádět. Pro použití v následných algoritmech je zde pouze pojmenujme, a to FindLexBoundMax(H, c, j) a BoundedLPMax(H, c, m1, m2). Doposud získané poznatky o řešení úlohy lineárního programování v rovině můžeme shrnout ve schématu na obrázku 4.3. Zbývá tedy uvést algoritmy řešící jednotlivé podproblémy z tohoto schématu. Ty však můžeme snadno formulovat na základě výsledků získaných v této kapitole. Věta 4.2. Časová složitost algoritmu 3 je O(n). Důkaz. Algoritmus obsahuje tři for cykly, přičemž při vykonání algoritmu proběhnou vždy právě dva. Časová složitost každého z těchto for cyklů je O(n). Ostatní části algoritmu mají konstantní časovou složitost. Proto je celková časová složitost algoritmu 3 rovna O(n). Věta 4.3. Časová složitost algoritmu 4 je O(n). Důkaz. Zřejmě vždy proběhne maximálně jeden for cyklus časové složitosti O(n), a proto je i celková časová složitost algoritmu lineární. Věta 4.4. Časová složitost algoritmu 5 je O(n). Důkaz. Algoritmus obsahuje dva nevnořené cykly for. Časová složitost prvního cyklu je O(n). Mohutnost množiny I je nejvýše n, a proto i časová složitost cyklu for all je O(n). Ostatní části algoritmu mají konstantní časovou složitost. Celková časová složitost algoritmu 5 je tedy O(n). 28 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Existuje s splňující (4.1)? • ne, 2 svědci hl a hr → Použij svědky jako m1, m2 a spusť algoritmus 2. ◦ řešení ◦ úloha nemá přípustné řešení • ne, 1 svědek hj → Existuje hk tak, že ({hj, hk}, c) má lexikograficky nejmenší řešení? • ano → Použij hj, hk jako m1, m2 a spusť algoritmus 2. ◦ řešení ◦ úloha nemá přípustné řešení • ne → Existuje hk tak, že ({hj, hk}, c) má lexikograficky největší řešení? • ano → Použij hj, hk jako m1, m2 a spusť modifikovaný algoritmus 2. ◦ lexikograficky největší řešení (nejmenší neexistuje) ◦ úloha nemá přípustné řešení • ne → Je H = ∅? ◦ ano → úloha nemá přípustné řešení ◦ ne → nějaké optimální řešení (lex. nejm. ani nejv. neexistuje) • ano → Je H = ∅? ◦ ano → úloha nemá přípustné řešení ◦ ne → úloha je neomezená Obrázek 4.3: Schéma shrnující doposud získané poznatky 29 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Algoritmus 3 Algoritmus k nalezení vektoru s splňujícího podmínky (4.1) Findvector(H, c) Vstup: Úloha lineárního programování (H, c). Výstup: Vektor s splňující podmínky (4.1), pokud takový existuje. V opačném případě algoritmus vrací hodnotu nil, navíc do proměnných3 wit1, wit2 vloží indexy svědků omezenosti úlohy, v případě jednoho svědka vloží do wit2 hodnotu nil. 1: if cy = 0 then 2: for i = 1 to n do 3: Fi ← ai2cx cy − ai1; Gi ← ai2 cy ; 4: else 5: for i = 1 to n do 6: Fi ← −ai2; Gi ← ai1 cx ; 7: wl ← −∞; wr ← ∞; 8: wit1 ← nil; wit2 ← nil; 9: for i = 1 to n do 10: if (Fi > 0) ∧ Gi Fi > wl then 11: wl ← Gi Fi ; wit1 ← i; 12: else if (Fi < 0) ∧ Gi Fi < wr then 13: wr ← Gi Fi ; wit2 ← i; 14: else 15: if (Fi = 0) ∧ (Gi > 0) then 16: wit1 ← i; wit2 ← nil; 17: return nil 18: if wl > wr then 19: return nil 20: else 21: if wl = −∞ then 22: w ← wl; 23: else if wr = ∞ then 24: w ← wr; 25: else 26: w ← 0; 27: if cy = 0 then 28: return w; 1 cy − cx cy w 29: else 30: return 1 cx ; w 30 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Algoritmus 4 Algoritmus ke zjištění existence lexikograficky nejmenšího řešení FindLexBound(H, c, j) Vstup: Úloha lineárního programování (H, c), index j poloroviny hj s hraniční přímkou kolmou na vektor c. Výstup: Index k poloroviny hk takové, že úloha ({hj, hk}, c) má lexikograficky nejmenší řešení, nebo hodnota nil, pokud taková polorovina neexistuje. 1: if aj2 = 0 then 2: for all i ∈ {1, . . . , n} \ {j} do 3: if ai1aj2−ai2aj1 aj2 < 0 then 4: return i 5: else 6: for all i ∈ {1, . . . , n} \ {j} do 7: if ai2 < 0 then 8: return i 9: return nil Jestliže je úloha neomezená, požadujeme, aby výstupem algoritmu byla polopřímka ρ, která je celá obsažena v množině všech přípustných řešení a podél níž účelová funkce neomezeně roste. Předpokládejme tedy, že úloha je neomezená, na základě algoritmu 3 jsme získali vektor s = (sx, sy) splňující podmínky (4.1) a z algoritmu 5 máme bod p0 ∈ H , p0 = [p0x, p0y]. Vycházejme z důkazu věty 4.1. Pro každou polorovinu hi ∈ H \ H existuje λi takové, že p0 + λs leží v polorovině hi pro libovolné λ ≥ λi. Zřejmě toto splňuje takové λi, že p0 +λis je průsečík přímky li s polopřímkou ρ0 = {p0 +λs | λ ∈ R; λ ≥ 0}, pokud existuje. V opačném případě leží v polorovině hi celá polopřímka ρ0 neboť vektor s splňuje vlastnosti (4.1). Průsečík λi snadno vypočteme. Dosazením do rovnice přímky li získáme rovnici ai1(p0x + λisx) + ai2(p0y + λisy) = bi, z níž po úpravách dostaneme výraz λi = bi − (ai1p0x + ai2p0y) ai1sx + ai2sy . Poznamenejme, že pro libovolné i, pro něž hi ∈ H \ H , je výraz ai1sx + ai2sy nenulový. Poloroviny hi ∈ H, pro které je tento výraz nulový, mají hraniční přímky rovnoběžné s vektorem s, a tudíž patří do H . Hledanou polopřímkou je tedy polopřímka ρ = {p + λs | λ ∈ R; λ ≥ 0}, kde p = p0 + λmaxs, λmax = max{λi | i ∈ {1, . . . , n} \ I }. K nalezení polopřímky ρ formulujme algoritmus 6. Jeho časová složitost je zřejmě lineární. 3 Proměnné wit1, wit2 chápeme jako globální proměnné. 31 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Algoritmus 5 Algoritmus zjišťující, zda je průnik podmnožiny G množiny H prázdný IsEmpty(H, c, s) Vstup: Úloha lineárního programování (H, c) a vektor s zadávající podmnožinu G = {hi ∈ H | s, ni = 0}. Výstup: Jestliže G = ∅, je výstupem libovolný bod4 p0 ∈ G. Pokud G = ∅, výstupem je hodnota nil. 1: I ← ∅; 2: for i=1 to n do 3: if ai1sx + ai2sy = 0 then 4: I ← I ∪ {i}; 5: if sy = 0 then 6: j ← 1; c∗ ← cx; 7: else 8: j ← 2; c∗ ← cy; 9: wl ← −∞; wr ← ∞; 10: for all i ∈ I do 11: if aij > 0 then 12: wr ← min(wr, bi aij ); 13: else if aij < 0 then 14: wl ← max(wl, bi aij ); 15: if wl > wr then 16: return nil 17: else 18: if wl = −∞ ∧ wr = ∞ then 19: if c∗ ≤ 0 then 20: w ← wl; 21: else 22: w ← wr; 23: else if wl = −∞ then 24: w ← wl; 25: else if wr = ∞ then 26: w ← wr; 27: else 28: w ← 0; 29: if sy = 0 then 30: return [w, 0] 31: else 32: return [0, w] 32 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Algoritmus 6 Algoritmus k nalezení polopřímky, na níž je úloha neomezená FindRho(H, c, s, p0) Vstup: Neomezená úloha lineárního programování (H, c), vektor s splňující podmínky (4.1) a bod p0 ∈ H . Výstup: Bod p takový, že polopřímka ρ = {p + λs | λ ∈ R; λ ≥ 0} leží v množině všech přípustných řešení a účelová funkce podél ní neomezeně roste. 1: λmax ← 0; 2: for all i ∈ {1, . . . , n} \ I do 3: λmax ← max λmax, bi−(ai1p0x+ai2p0y) ai1sx+ai2sy ; 4: return p0 + λmaxs S využitím doposud odvozených algoritmů jsme schopni uvést algoritmus 7 k řešení obecné úlohy lineárního programování v rovině. Věta 4.5. Časová složitost algoritmu 7 je O(n2 ). Důkaz. Tvrzení přímo vyplývá z časové složitosti algoritmů využitých v algoritmu 7. 4 V situaci, kdy zjišťujeme, zda-li je H = ∅, je tento bod libovolným optimálním řešením. 33 KAPITOLA 4. NEOMEZENÁ ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ Algoritmus 7 Algoritmus k řešení úlohy lineárního programování v rovině LP(H, c) Vstup: Úloha lineárního programování (H, c). Výstup: Lexikograficky nejmenší optimální řešení úlohy (největší, pokud neexistuje, či nějaké, pokud ani nejv. neexistuje), zpráva o neomezenosti úlohy a polopřímka, na níž je úloha neomezená, nebo zpráva o tom, že úloha nemá přípustné řešení. 1: s ← FindVector(H, c); 2: if s = nil then 3: if wit2 = nil then 4: return BoundedLP(H \ {hwit1 , hwit2 }, c, hwit1 , hwit2 ) 5: else 6: k ← FindLexBound(H, c, wit1); 7: if k = nil then 8: return BoundedLP(H \ {hj, hk}, c, hwit1 , hk) 9: else 10: k ← FindLexBoundMax(H, c, wit1); 11: if k = nil then 12: return BoundedLPMax(H \ {hj, hk}, c, hwit1 , hk) 13: else 14: p ← IsEmpty(H, c, (cy, −cx)); 15: if p = nil then 16: return Úloha nemá přípustné řešení. 17: else 18: return p 19: else 20: p0 ← IsEmpty(H, c, s); 21: if p0 = nil then 22: return Úloha nemá přípustné řešení. 23: else 24: p ← FindRho(H, c, s, p0) 25: return Úloha je neomezená, a to na polopřímce ρ = {p + λs | λ ∈ R; λ ≥ 0}. 34 Kapitola 5 Očekávaná časová složitost Viděli jsme, že v nejhorším případě může být délka výpočtu algoritmu k řešení úlohy lineárního programování kvadratická vzhledem k velikosti vstupu. Délka výpočtu dané úlohy závisí na pořadí polorovin. Pro některá uspořádání může být kvadratická, pro jiná lineární. V této části upravíme algoritmus do podoby náhodnostního algoritmu a budeme se zabývat očekávanou časovou složitostí. Pro danou úlohu poloroviny náhodně uspořádáme tak, že všechny permutace polorovin budou mít stejnou pravděpodobnost zvolení. K nalezení náhodné permutace daného pole formulujme algoritmus 8. Předpokládáme, že máme k dispozici generátor náhodných čísel Random(k), jehož výstup je náhodné celé číslo z intervalu 1, k , a to každé se stejnou pravděpodobností. Časová složitost algoritmu 8 je zřejmě lineární. Algoritmus 8 Algoritmus pro náhodnou permutaci RndPerm(H, c) Vstup: Pole A = (A[1], . . . , A[n]) Výstup: Náhodná permutace pole A 1: for i = 1 to n − 1 do 2: rand ← n − Random(n − i + 1) + 1; 3: pom ← A[i]; 4: A[i] ← A[rand]; 5: A[rand] ← pom; 6: return A Pomocí procedury pro náhodnou permutaci upravme algoritmus 7 do podoby náhodnostního algoritmu. Jediná změna v algoritmu 9 oproti algoritmu 7 se nachází na řádcích 4, 8 a 12. Délka výpočtu algoritmu 9 závisí na náhodné volbě permutace polorovin. Příčinou možné kvadratické délky výpočtu je algoritmus pro omezenou úlohu lineárního programování na řádku 4, respektive 8 či 12. Analyzujme tedy délku výpočtu algoritmu 2 v případě náhodnostního algoritmu. Důvodem kvadratické časové složitosti algoritmu 2 je počet možných vykonání 35 KAPITOLA 5. OČEKÁVANÁ ČASOVÁ SLOŽITOST Algoritmus 9 Náhodnostní algoritmus k řešení úlohy lin. programování v rovině RndLP(H, c) Vstup: Úloha lineárního programování (H, c). Výstup: Lexikograficky nejmenší optimální řešení úlohy (největší, pokud neexistuje, či nějaké, pokud ani nejv. neexistuje), zpráva o neomezenosti úlohy a polopřímka, na níž je úloha neomezená, nebo zpráva o tom, že úloha nemá přípustné řešení. 1: s ← FindVector(H, c); 2: if s = nil then 3: if wit2 = nil then 4: return BoundedLP(RndPerm(H \ {hwit1 , hwit2 }), c, hwit1 , hwit2 ) 5: else 6: k ← FindLexBound(H, c, wit1); 7: if k = nil then 8: return BoundedLP(RndPerm(H \ {hj, hk}), c, hwit1 , hk) 9: else 10: k ← FindLexBoundMax(H, c, wit1); 11: if k = nil then 12: return BoundedLPMax(RndPerm(H \ {hj, hk}), c, hwit1 , hk) 13: else 14: p ← IsEmpty(H, c, (cy, −cx)); 15: if p = nil then 16: return Úloha nemá přípustné řešení. 17: else 18: return p 19: else 20: p0 ← IsEmpty(H, c, s); 21: if p0 = nil then 22: return Úloha nemá přípustné řešení. 23: else 24: p ← FindRho(H, c, s, p0) 25: return Úloha je neomezená, a to na polopřímce ρ = {p + λs | λ ∈ R; λ ≥ 0}. 36 KAPITOLA 5. OČEKÁVANÁ ČASOVÁ SLOŽITOST druhé větve podmínky na řádku 3. Časová složitost i-té iterace cyklu for (ř. 2– 35) je v závislosti na splnění podmínky na řádku 3 buď konstantní, anebo O(i). V případě náhodnostního algoritmu je splnění této podmínky náhodnou veličinou. Proto pro i ∈ {1, . . . , n} označme Xi náhodnou veličinu, která nabývá hodnoty 1, jestliže vi−1 ∈ hi (časová složitost i-té iterace je lineární), nebo hodnoty 0, jestliže vi−1 ∈ hi (časová složitost i-té iterace je konstantní). Celková časová složitost vykonání druhé větve během všech iterací dohromady je proto n i=1 XiO(i), zatímco celková časová složitost vykonání první větve během všech iterací je lineární, a tedy časová složitost algoritmu 2 je rovna náhodné veličině O(n) + n i=1 XiO(i). Definice 5.1. Nechť časová složitost algoritmu A pro daný vstup je náhodná veličina Y , pak očekávanou časovou složitostí algoritmu A rozumíme střední hodnotu náhodné veličiny Y, tj. E(Y ). Věta 5.1. Očekávaná časová složitost algoritmu 9 je O(n). Důkaz. Analyzujme očekávanou časovou složitost algoritmu 2. Jestliže ukážeme, že je lineární, bude lineární i očekávaná časová složitost algoritmus 9, neboť ostatní části algoritmu kromě řádku 4, respektive 8, mají lineární časovou složitost. Z předchozích výsledků vyplývá, že očekávaná časová složitost algoritmu 9 pro daný vstup velikosti n je1 E O(n) + n i=1 XiO(i) = O(n) + n i=1 E(Xi)O(i). Náhodná veličina Xi nabývá hodnot 0, nebo 1. Má tedy alternativní rozdělení pravděpodobnosti, a proto střední hodnota náhodné veličiny Xi je rovna pravděpodobnosti, že náhodná veličina Xi nabývá hodnoty 1, tj. E(Xi) = P(Xi = 1) = P(vi−1 ∈ hi). Podívejme se na algoritmus 2 po i iteracích hlavního cyklu. Máme tedy nalezeno optimální řešení vi úlohy (Hi, c). Bod vi je průsečíkem nejméně dvou z hraničních přímek polorovin z Hi. Odebereme-li jednu polorovinu, pak může dojít ke změně optimálního řešení pouze v případě, že bod vi leží na hraniční přímce této poloroviny. Pokud vi leží na hraničních přímkách právě dvou polorovin, pak odebráním jedné z nich dojde ke změně množiny všech přípustných řešení a může dojít ke změně 1 Využijme linearitu střední hodnoty náhodné veličiny, tedy že pro náhodnou veličinu Y a konstanty a, b ∈ R platí E(a + bY ) = a + bE(Y ). 37 KAPITOLA 5. OČEKÁVANÁ ČASOVÁ SLOŽITOST c vi (a) vi průsečíkem hraničních přímek právě dvou polorovin c vi (b) vi průsečíkem hraničních přímek více než dvou polorovin Obrázek 5.1: Bod vi jako průsečík hraničních přímek polorovin optimálního řešení. Situaci je znázorněna na obrázku 5.1(a). Jestliže vi je průsečíkem hraničních přímek více než dvou polorovin, pak pouze dvě z nich tvoří hranici množiny všech přípustných řešení, jak ilustruje obrázek 5.1(b). Pouze odebráním jedné z těchto dvou polorovin dojde ke změně množiny všech přípustných řešení a může dojít ke změně optimálního řešení. Odebrání jiné poloroviny, tedy takové, na jejíž hraniční přímce bod vi neleží, nevede v žádné z předchozích situací ke změně optimálního řešení. Pouze pro dvě poloroviny z Hi tedy platí, že odebráním jedné z nich může dojít ke změně optimálního řešení. Situace, kdy vi−1 ∈ hi nastane právě tehdy, když odebráním poloroviny hi dojde ke změně optimálního řešení (neboli po i − 1 iteracích je jiné optimální řešení než po i iteracích). Tedy vi−1 ∈ hi může nastat pouze tehdy, pokud hi je jednou ze dvou výše uvedených polorovin, jejichž odebráním může dojít ke změně optimálního řešení. Jelikož jsou poloroviny uspořádány náhodně, je pravděpodobnost, že hi je jednou z těchto polorovin, nejvýše 2 i . Proto E(Xi) = P(vi−1 ∈ hi) ≤ 2 i . Dosazením do vztahu pro očekávanou časovou složitost algoritmu 2 dostáváme O(n) + n i=1 E(Xi)O(i) ≤ O(n) + n i=1 2 i O(i) = O(n), a tedy očekávaná časová složitost je lineární. Zdůrazněme, že lineární časová složitost je očekávána pro každý vstup. Nejedná se o průměrnou časovou složitost přes všechny možné vstupy. Očekávání se vztahuje k náhodné volbě obsažené v algoritmu, nikoliv k možnému vstupu. Pro libovolnou úlohu je očekávaná časová složitost algoritmu 9 lineární. 38 Kapitola 6 Simplexová metoda K řešení úloh lineárního programování se většinou využívá simplexová metoda. Stručně ji proto v této části popíšeme a srovnáme s náhodnostním algoritmem odvozeným v této práci. Vycházejme z [3, kapitola 2], kde lze také nalézt podrobnější popis simplexové metody. Abychom mohli řešit úlohu lineárního programování simplexovou metodou, je třeba ji nejprve převést do tzv. standardního tvaru, tedy do tvaru max{cx | Ax = b, x ≥ 0}1 . Konkrétně v případě našeho výchozího tvaru (podle definice 1.2) úlohu max{cx | Ax ≤ b} převedeme na ekvivalentní úlohu max    (c, −c) x x (A | −A | E)   x x x∗   = b,   x x x∗   ≥ 0    . Původní úlohu jsme transformovali tak, že jsme přidali vektor doplňkových proměnných x∗ , x∗ ≥ 0 a převedli tak soustavu nerovnic Ax ≤ b na soustavu rovnic Ax + x∗ = b. Dále jsme každou z původních proměnných xi nahradili dvojicí nezáporných proměnných xi a xi tak, že xi = xi − xi . Vidíme, že transformací úlohy lineárního programování v rovině do standardního tvaru získáme úlohu s n+4 proměnnými. Chceme-li tedy řešit úlohu lineárního programování v rovině pomocí simplexové metody, musíme řešit úlohu vyšší dimenze. Množina všech přípustných řešení úlohy lineárního programování je průnikem konečně mnoha afinních poloprostorů2 a představuje proto polyedr. Simplexová metoda postupně prochází vrcholy tohoto polyedru a hledá optimální řešení. Základním krokem je nalezení výchozího řešení, tj. libovolného přípustného řešení, které představuje vrchol polyedru. Dále zkoumáme, jak by se měnila hodnota 1 Zápisem x ≥ 0 pro vektor x = (x1, . . . , xd)T rozumíme xi ≥ 0, i = 1, . . . , d. 2 Afinním poloprostorem nazýváme množinu vektorů {x ∈ Rd | ax ≤ b} pro nějaký vektor a ∈ Rd a číslo b ∈ R. 39 KAPITOLA 6. SIMPLEXOVÁ METODA účelové funkce, pokud bychom se pohybovali po některé z hran polyedru vycházející z aktuálního vrcholu. Pokud přírůstek hodnoty účelové funkce pro žádnou hranu není kladný, pak se nacházíme v optimálním řešení. V opačném případě se vydáme po některé hraně polyedru, podél níž hodnota účelové funkce roste. Pokud je tato hrana polopřímkou, pak je úloha neomezená. Jestliže se jedná o úsečku, pak dojdeme do jiného vrcholu polyedru, který představuje jiné přípustné řešení s vyšší hodnotou účelové funkce. Předchozí postup opakujeme pro nově nalezené řešení. Přestože při použití v praxi je simplexová metoda poměrně efektivní nástroj pro řešení úlohy lineárního programování, v nejhorším případě může délka výpočtu růst exponenciálně vzhledem k velikosti vstupu (viz např. [4, kapitola 5.2]). Časová složitost simplexové metody je tedy O(2n ). Jak simplexová metoda, tak předchozí náhodnostní algoritmus jsou iteračními algoritmy. Zatímco simplexová metoda v jednotlivých iteracích prochází přípustná řešení odpovídající vrcholům polyedru, řešení nalezená během iterací náhodnostního algoritmu nemusí být přípustnými řešeními původní úlohy. Každý algoritmus tak k problému přistupuje zcela odlišným způsobem, což vede také k odlišným závěrům o časových složitostech obou algoritmů. 40 Kapitola 7 Úloha lineárního programování ve vyšších dimenzích Doposud jsme se zabývali úlohou lineárního programování v rovině, popřípadě na přímce. Odvozený algoritmus lze však zobecnit a na stejném principu je možné řešit i úlohy ve vyšších dimenzích. Každá nerovnice v úloze dimenze d zadává afinní poloprostor v Rd . Stejně jako v dvourozměrném případě úlohu řešíme iteračně postupným přidáváním nerovnic. Pro i = 1, 2, . . . , n označme hi jednotlivé poloprostory, gi jejich hranice. Lze ukázat že i ve vyšších dimenzích platí analogie věty 3.1, viz [1, kapitola 4.6]. Je-li vi řešení úlohy dimenze d po i iteracích, pak se přidáním poloprostoru hi+1 optimální řešení nezmění, pokud vi v tomto poloprostoru leží. V opačném případě se nové řešení nachází na hranici gi+1 nově přidaného poloprostoru, anebo úloha nemá přípustné řešení. Tato hranice je nadrovinou dimenze d − 1. Řešení vi+1 musí ležet ve všech průnicích hranice gi+1 s poloprostory h1, . . . , hi. Na této množině tak hledáme bod, ve kterém účelová funkce nabývá maxima. Abychom tedy nalezli řešení úlohy po i + 1 iteracích, musíme vyřešit úlohu dimenze d − 1. Účelovou funkci této úlohy snadno získáme omezením původní účelové funkce na množinu gi+1. Uvažme úlohu lineárního programování v trojrozměrném euklidovském prostoru. Nejprve musíme zjistit, zda je úloha neomezená. Vycházíme přitom z analogie věty 4.1. Hledáme proto vhodný vektor s určující v množině všech přípustných řešení polopřímku, podél níž účelová funkce neomezeně roste. Požadujeme tedy, aby odchylka vektorů s a c byla v intervalu 0, π 2 a současně odchylka vektorů s a ni byla v intervalu 0, π 2 pro i = 1, 2, . . . , n, kde ni je normálový vektor hraniční roviny poloprostoru hi směřující do tohoto poloprostoru. Jelikož nám záleží pouze na směru vektoru s, můžeme požadovat, aby jeho koncový bod ležel v dané rovině kolmé na vektor c. Každému směru, ve kterém účelová funkce roste, odpovídá právě jeden bod této roviny. Vyjádřením výše uvedených požadavků na vektor s pak získáme úlohu lineárního programování v rovině. Tu umíme řešit v očekávaném lineárním čase. Jestliže nalezneme vektor s požadovaných vlastností, původní úloha je neomezená. V opačném případě můžeme během výpočtu získat svědky skutečnosti, že dvoudimenzionální úloha k nalezení vektoru s nemá přípustné řešení, a tedy svědky 41 KAPITOLA 7. ÚLOHA LINEÁRNÍHO PROGRAMOVÁNÍ VE VYŠŠÍCH DIMENZÍCH omezenosti původní úlohy. Až na singulární případy jsou tito svědci 3. Získáme je ve chvíli, kdy po přidání poloroviny hi jednodimenzionální úloha na přímce li k nalezení řešení vi nemá přípustné řešení. Konkrétně se jedná o polorovinu hi společně se svědky neexistence přípustného řešení jednodimenzionální úlohy. Singulární případy bychom řešili analogicky jako v kapitole 4 při zjišťování existence lexikograficky nejmenšího řešení. Máme-li k dispozici svědky omezenosti úlohy v R3 , můžeme přejít k řešení omezené úlohy. Vyjdeme z počátečního řešení, které je průnikem hraničních rovin svědků omezenosti. Postupně v náhodném pořadí přidáváme poloprostory. Předpokládejme že po i iteracích máme řešení vi. Pokud vi leží v poloprostoru hi+1, pak vi+1 = vi. V opačném případě řešíme úlohu lineárního programování v hraniční rovině gi+1 poloprostoru hi+1. Jednotlivé poloroviny této úlohy jsou určeny jako průniky roviny gi+1 s poloprostory h1, . . . , hi. Počet polorovin v této úloze je i, a proto její očekávaná časová složitost je O(i). Očekávaná délka výpočtu i-té iterace je tedy buď konstantní, pokud vi−1 ∈ hi, anebo O(i). Situace je analogická odvození očekávané časové složitosti v kapitole 5. Náhodnostní algoritmus k řešení úlohy lineárního programování v R3 má proto lineární očekávanou časovou složitost. Podobný postup lze využít i pro úlohy dimenze d. Zde získáme d svědků omezenosti úlohy. V jednotlivých iteracích řešíme úlohy dimenze d−1. Indukcí tak získáme stejný závěr o lineární očekávané časové složitosti náhodnostního algoritmu k řešení úlohy lineárního programování v Rd . Očekávaná časová složitost uvedeného algoritmu je lineární pro libovolnou fixní dimenzi d. To však neznamená, že délka výpočtu je stejná pro všechny dimenze. Lze ukázat, že očekávaná časová složitost algoritmu roste exponenciálně s rostoucí dimenzí, viz [1, kapitola 4.6]. Algoritmus je proto vhodný pouze pro úlohy nižších dimenzí. 42 Kapitola 8 Implementace Algoritmus k řešení úlohy lineárního programování v rovině odvozený v této práci byl implementován v prostředí Delphi 6.0. Vznikl tak program využívající teoretické poznatky získané v předchozích kapitolách. Předností aplikace je její názornost, a to díky grafickému výstupu. Je zaměřena na demonstraci principu algoritmu, a proto je vhodná pro výukové účely. Program umožňuje přehledné znázornění jednotlivých kroků výpočtu algoritmu. Aplikace je určena pro operační systém Windows a spouští se pomocí souboru LPinPlane.exe, který se nachází na přiloženém kompaktním disku. Po spuštění zadáme vektor c představující účelovou funkci a jednotlivé nerovnice úlohy. Pomocí tlačítek můžeme nerovnice uspořádat do námi požadovaného pořadí, anebo je uspořádat náhodně. Úlohu je možné uložit do souboru, ze kterého ji můžeme později otevřít, a nemusíme ji tak zadávat po jednotlivých nerovnicích. Úlohu si můžeme nechat vykreslit kliknutím na příslušné tlačítko. Poloroviny jsou zobrazovány stejným způsobem jako v této práci, tedy hraniční přímkou a vystínováním příslušné strany. Po kliknutí na některou z nerovnic je odpovídající hraniční přímka zvýrazněna. Režim zadávání nerovnic je znázorněn na obrázku 8.1. V nabídce Nástroje - Možnosti můžeme volit mezi omezenou a obecnou úlohou. V případě omezené úlohy se jedná o implementaci algoritmu 2, kde za dodatečná lineární omezení jsou považovány první dvě poloroviny. Uživatel tak musí sám vhodně vybrat první dvě nerovnice, a to tak, aby společně s vektorem c tvořily omezenou úlohu s lexikograficky nejmenším řešením. Naopak při volbě obecné úlohy si program najde svědky omezenosti sám a změní pořadí polorovin tak, že svědky umístí na první dvě pozice. Při volbě obecné úlohy se program dokáže vypořádat i s neomezenou úlohou. Jedná se o implementaci algoritmu 9. Po spuštění programu je zvolena obecná úloha. Kliknutím na tlačítko Demonstrace / Řešení přejdeme do demonstračního režimu. V případě, že v průběhu výpočtu dojde k volání algoritmu 2, můžeme zde procházet jednotlivé kroky algoritmu, a to buď ručně, anebo automaticky v nastaveném časovém intervalu. Průběžně jsou zobrazovány jednotlivé dílčí úlohy a jejich řešení. Hraniční přímka nově přidané poloroviny je zvýrazněna. Program v průběhu demonstrace je znázorněn na obrázku 8.2. V případě, že k volání algoritmu 2 nedo- 43 KAPITOLA 8. IMPLEMENTACE jde (např. u neomezené úlohy), je výsledek zobrazen ihned. Aplikace také umožňuje nastavení některých parametrů ovlivňujících její funkci a vzhled. Kromě již zmíněné volby mezi omezenou a obecnou úlohou lze nastavit časový interval mezi jednotlivými kroky při automatické demonstraci či barvy jednotlivých prvků grafického výstupu. Program umožňuje zkoumat, jak se mění průběh výpočtu algoritmu v závislosti na pořadí vstupních nerovnic. Je možné srovnávat různá uživatelem určená pořadí s náhodným pořadím, což může být přínosné pro objasnění časové složitosti algo- ritmu. Vytvořená aplikace tedy není pouhou strohou implementací algoritmu, ale zahrnuje také některé užitečné funkce usnadňující pochopení jeho principu, a to především prostřednictvím přehledného grafického výstupu. 44 KAPITOLA 8. IMPLEMENTACE Obrázek 8.1: Program LPinPlane.exe v režimu zadávání nerovnic Obrázek 8.2: Program LPinPlane.exe v průběhu demonstračního režimu 45 Seznam použité literatury [1] BERG, Mark, et al. Computational Geometry : Algorithms and Applications. Third Edition. Berlin : Springer-Verlag, 2008. 386 s. ISBN 978-3-540-77973-5. [2] CORMEN, Thomas H., et al. Introduction to algorithms. 2nd ed. Cambridge (Mass.) : MIT Press, 2001. 1180 s. ISBN 0-262-03293-7. [3] JIRSÍK, Tomáš. Simplexová metoda a její využití v ekonomii [online]. Brno : Masarykova univerzita, 2009. 44 s. Bakalářská práce. Masarykova univerzita, Přírodovědecká fakulta. Dostupné z WWW: . [4] KASANA, Harvir Singh; KUMAR, Krishna Dev. Introductory Operations Research : Theory and Applications. Berlin : Springer-Verlag, 2004. 581 s. ISBN 3-540-40138-5. 46