Simplex metoda

Do čitanja ove teme, već ste savladali strukturu linearnog programa i sastavljanje iste temeljem opisa slučaja. Osim toga, naučili ste grafičko rješavanje linearnog programa.

Podsjetimo se, opći oblik linearnog programa:

\[ \text{max ili min } Z=c_1x_1+c_2x_2+⋯+c_nx_n\]

uz ograničenja:

\[a_{11}x_1+a_{12}x_2+⋯+a_{1n}x_n≤b_1\] \[\dots\] \[a_{m1}x_1+a_{m2}x_2+⋯+a_{mn}x_n≤b_m\] \[xj≥0 \] \[\text{za sve j}\]

Također, znamo da se za modele s dvije varijable rješenje može tražiti grafički i da je optimalno rješenje na vrhovima poligona koji oblikuje prostor izvedivih rješenja (engl. feasible region).

Intuicija Simplex metode

Simplex je algoritam koji sustavno prelazi iz jednog vrha u drugi, sve dok ne dođe do najboljeg (optimalnog) – i to samo po rubovima prostora izvedivih rješenja.

To čini tako da:

U svakom koraku provjerava ima li susjedni vrh bolji cilj.
Ako ima, prelazi tamo.
Ako nema – stajemo, našli smo optimum.

Simplex zapravo provodi algebarsku verziju onoga što grafička metoda radi za varijantu modela s dvije varijable – šeta po vrhovima dok ne pronađe onaj u kojem jevrijednost funkcije cilja maksimalna (ili minimalna). Prednost Simplex metode je u tome što, za razliku od grafičkog pristupa, može rješavati i probleme s više od dvije varijable odluke.

Simplex metodu osmislio je George Dantzig, američki matematičar i istraživač, 1947. godine.

Dantzig je tada radio za U.S. Air Force, unutar RAND Corporation, gdje je rješavao probleme planiranja i logistike u sklopu poslova vojnog planiranja nakon Drugog svjetskog rata.

U to doba se identificirao problem: kako maksimalno iskoristiti ograničene resurse (npr. gorivo, kapaciteti, transport) da bi se postigao što bolji ishod.
Linearni modeli već su bili poznati, ali nije postojao učinkovit algoritam za njihovo rješavanje.
Dantzig je osmislio metodu kojom se iterativno dolazi do optimalnog rješenja – krećući se od vrha do vrha u izvedivom prostoru rješenja – i to je bio početak Simplex metode.

– Temeljna ideja je i kasnije ostala ista, ali su nastale razne varijante i nadogradnje, posebno za:

Degeneraciju i cikluse: rijetki, ali mogući problemi s beskonačnom petljim \(\rightarrow\) uvedeno Blandovo pravilo za sprječavanje ciklusa.
Artificijelne varijable - kad početno rješenje nije poznato \(\rightarrow\) razvijene su metode poput:
- Dvostruke Simplex metode
- Big-M metode
- Dvostupanjske metode (two-phase simplex)
Računalna implementacija:
- Klasični Simplex koristi tablični prikaz koji je pregledan ali neučinkovit na računalu.
- Revidirani Simplex (Revised Simplex) – koristi matricu i vektore, efikasniji u memoriji (koristi se u lpSolve paketu)
- Dualni Simplex – koristi se kad su dualne vrijednosti poznate ili kad se ograničenja mijenjaju (koristi se u linprog paketu).
Iako je Simplex vrlo učinkovit u praksi, u najgorem slučaju ima eksponencijalnu složenost. Zato se za vrlo velike probleme danas koriste:
- Interior Point metode (razvijene 1984., Karmarkar); (koristi se u linprog paketu)
- Hibridni pristupi - kombinacija Simplexa i heuristika (npr. u optimizaciji rute, strojnom učenju, planiranju proizvodnje)

Postupak možemo formalno zapisati na sljedeći način

Neka je zadan linearni program:

\[ \max \; z = c^T x \]

uz ograničenja:

\[ Ax \leq b, \quad x \geq 0, \]

gdje je
- \(A \in \mathbb{R}^{m \times n}\) matrica koeficijenata u ograničenjima,
- \(b \in \mathbb{R}^m\) vektor desnih strana ograničenja,
- \(c \in \mathbb{R}^n\) vektor koeficijenata funkcije cilja,
- \(x \in \mathbb{R}^n\) vektor varijabli odluka.

Definicija:
Simplex metoda je iterativni algoritam koji generira niz osnovnih dopustivih rješenja

\[ x^{(0)}, x^{(1)}, \dots, x^{(k)}, \]

pri čemu svako \(x^{(t)}\) zadovoljava ograničenja \(Ax = b, \; x \geq 0\), a vrijednost funkcije cilja se strogo povećava:

\[ c^T x^{(t+1)} > c^T x^{(t)}, \]

dok se ne postigne optimalno rješenje \(x^*\) za koje vrijedi:

\[ \bar{c}_j = c_j - c_B^T A_B^{-1} A_j \leq 0 \quad \forall j \in N, \]

gdje su \(\bar{c}_j\) reducirani troškovi.

Pojašnjenje. Da bi Simplex mogao “odlučiti” daje li sljedeći vrh bolju vrijednost funkcije cilja, mora izračunati reducirane troškove za sve varijable koje trenutno nisu u bazi (tzv. ne-bazične varijable). Reducirani trošak \(\bar{c}_j\) varijable \(x_j\) govori nam koliko bi se promijenila vrijednost funkcije cilja ako unesemo jednu jedinicu varijable \(x_j\) u rješenje. Ako je \(\bar{c}_j > 0\), uvođenje \(x_j\) vodi višim vrijednostima funkcije cilja (maksimizacija); a ako je \(\bar{c}_j \leq 0\), uvođenje \(x_j\) ne poboljšava cilj. Tako \(c_j - c_B^T A_B^{-1} A_j\) predstavlja neto doprinos varijable \(x_j\). Kad su svi reducirani troškovi \(\bar{c}_j \leq 0\), to znači da niti jedna nebazična varijabla ne može poboljšati cilj i uvjet optimalnosti je dosegnut.

Zanimljivost. Naziv reducirani troškovi može zbunjivati, posebno kod problema maksimizacije, jer potječe iz povijesnih razloga. Simplex metodu je George Dantzig izvorno razvio za rješavanje vojnih problema logistike gdje se minimizirala potrošnja resursa poput goriva ili troškova transporta, pa su koeficijenti u funkciji cilja doista predstavljali stvarne troškove po jedinici. Kada se metoda formalizirala i proširila na sve tipove linearnih programa, terminologija se zadržala iz tradicije, iako se u literaturi ponekad koriste i nazivi poput marginalna cijena, relativni doprinos ili jednostavno reducirani koeficijenti. U suštini, reducirani “trošak” zapravo pokazuje neto promjenu vrijednosti funkcije cilja kad uvedemo jednu jedinicu varijable \(x_j\) u rješenje - što kod maksimizacije više odgovara pojmu marginalnog dobitka nego troška, ali matematička zajednica je zadržala izvorno nazivlje.

Algoritamski pristup

Formulacija problema
Prevedi problem u standardni oblik: \(\max \; c^T x\), uz \(Ax = b, \; x \geq 0\).
Inicijalizacija
Odaberi početnu bazu \(B\) (ako postoji), te izračunaj bazno rješenje
\[ x_B = A_B^{-1} b, \quad x_N = 0. \]
Izračun redukcija
Za sve nebazične varijable \(j \in N\) izračunaj reducirane troškove:
\[ \bar{c}_j = c_j - c_B^T A_B^{-1} A_j. \]
Test optimalnosti
Ako \(\bar{c}_j \leq 0 \; \forall j \in N\), tada je rješenje optimalno \(\rightarrow\) zaustavi algoritam.
Odabir ulazne varijable
Odaberi varijablu \(j \in N\) s \(\bar{c}_j > 0\) koja ulazi u bazu.
Odabir izlazne varijable
Izračunaj smjer \(d = A_B^{-1} A_j\).
Odredi ograničavajući omjer:
\[ \theta = \min_{i : d_i > 0} \frac{(x_B)_i}{d_i}. \]
Odgovarajuća bazična varijabla napušta bazu.
Pivotiranje
Ažuriraj bazu: \(B \leftarrow (B \setminus \{i\}) \cup \{j\}\).
Izračunaj novo bazično rješenje.
Ponovi
Vrati se na korak 3.

Provedba Simplex metode

Svaka iteracija metode sastoji se od nekoliko koraka:

Odabir ulazne varijable (entering variable) određuje se koja varijabla će ući u bazu kako bi se ciljna funkcija poboljšala.
Odabir izlazne varijable (leaving variable) određuje se koja trenutno bazična varijabla izlazi iz baze kako bi se omogućilo uvođenje nove varijable, a pritom očuvalo dopušteno rješenje.
Pivotiranje (update tablice) - ažurira se Simplex tablica: ulazna varijabla postaje nova bazična varijabla, a svi retci se linearno transformiraju kako bi nova baza bila valjana.
Provjera optimalnosti – ako više nema negativnih koeficijenata u Z-retku (ili uvjeta za poboljšanje), algoritam se zaustavlja – rješenje je optimalno.

Ovisno o situaciji, moguće su i dodatne provjere:

je li rješenje degenerirano
je li problem neograničen (engl. unbounded)
je li potrebno razriješiti vezu (kod jednakih omjera)

U nastavku za prikazuju detalji za svaki od ovih koraka.

Priprema

Priprema je prvi korak u provedbi Simplex metode i uvijek se provodi na jednak način, neovisno o odabirima u sljedećim koracima. Odnosi se na svođenje na standardni oblik, odnosno kanonski zapis.

Što točno obuhvaća kanonski zapis?

Funkcija cilja mora biti zapisana kao jednadžba s nulom na desnoj strani
Sva ograničenja moraju biti:
- jednadžbe
- dopunjene slack, surplus i po potrebi artificielnim varijablama
- Sve varijable moraju imati nenegativne vrijednosti

(Napomena: Slack varijable (dodaju se za ≤ ograničenja) mjere neiskorištene kapacitete, surplus varijable (za ≥ ograničenja) mjere prekoračenje ciljne vrijednosti, a artificijelne varijable služe za pronalaženje početnog rješenja kada to nije moguće direktno.)

Kako to postići?

Na primjer, ako rješavamo model:

\[\text{max } 3x_1+5x_2\] \[2x_1 +x_2 \leq 100\] \[x_1 + 2x_2 \leq 80\] \[x_1, x_2 \geq 0\]

Za funkciju cilja u danom primjeru to znači:

\[ Z = 3x_1+5x_2\text{ }⇒\]

\[Z - 3x_1-5x_2=0\]

Pretvaramo sve nejednakosti u jednakosti dodavanjem slack (ili dopunskih varijabli; varijable koje pokrivaju višak) varijabli, npr.

\[2x_1+x_2≤100\text{ }⇒\text{ } 2x_1+x_2+s_1=100\] \[s_1≥0\] \(2x_1+x_2≤100\) podrazumijeva da izraz \(2x_1+x_2\) može poprimiti vrijednosti manje od 100, npr. 50 ili 3. Dakle, ne možemo direktno samo zamijeniti manje ili jednako od s jednakosti. Slack varijabla \(s_1\) pokriva tu moguću razliku, tj. kad \(s_1\) pridodamo postojećem izrazu na lijevoj strani, ona označava moguću razliku do vrijednosti na desnoj strani. Obratite pozornost na to da ne možemo unaprijed znati koliko će to iznositi, pa je \(s_1\) zapravo još jedna nepoznanica čiju vrijednost utvrđujemo modelom. Ograničena je (\(s_1≥0\)) na način da je najmanja vrijednost koju može poprimiti nula (\(s_1=0\) ako bi vrijedilo \(2x_1+x_2=100\)) ili više od toga (\(s_1>0\) za situacije u kojima je \(2x_1+x_2<100\)). Isti princip primijenjujemo na drugo ograničenje.

Tako će ograničenja glasiti:

\[2x_1 +1\cdot x_2 +1\cdot s_1 + 0 \cdot s_2 = 100\] \[x_1 + 2x_2 + 0\cdot s_1 + 1 \cdot s_2 = 80\]

Ako imamo, npr. 2 varijable i 3 ograničenja, onda će model imati ukupno 2 osnovne + 3 slack varijable = 5 varijabli.

Svaka slack varijabla mjeri neiskorištenost ograničenja.

Simplex tablica

Simplex se provodi pomoću tablice (Simplex tableau), u kojoj su zapisana:

Ograničenja (s koeficijentima za svaku varijablu, uključujući slackove) - u recima
Funkcija cilja (Z) - uobičajeno u zadnjem retku
Trenutna vrijednost funkcije cilja
Bazične varijable (BV) - broj bazičnih varijabli jednak je broju ograničenja.

Za problem ranije zapisan u kanonskom obliku, početna simplex tablica će izgledati ovako:

BV	x₁	x₂	s₁	s₂	RHS
s₁	2	1	1	0	100
s₂	1	2	0	1	80
Z	-3	-5	0	0	0

U bazi se nalaze slack varijable – one čine inicijalno rješenje (npr. sve osnovne varijable = 0, svi slackovi = b vrijednosti).

Što su bazične varijable?

Bazične varijable su one koje se nalaze u trenutnom bazičnom rješenju
U tablici ih prepoznajemo kao varijable čiji stupci sadrže jednu jedinicu - odnosno, 1 u jednom retku, 0 u ostalima
Samo bazične varijable imaju nenultu vrijednost u RHS, sve druge su trenutno 0

Prema početnim postavkama, samo \(s_1\) i \(s_2\) zadovoljavaju kriterij bazičnih varijabli. Zbog toga su i zapisane u stupcu BV.

Odabir ulazne varijable

Odabir ulazne varijable (eng. entering variable) u Simplex metodi je odluka koja može:

utjecati na brzinu konvergencije
pomoći u izbjegavanju ciklusa
odrediti redoslijed prelazaka kroz vrhove/ presjeke ograničenja (a time i broj iteracija).

Napomena: ulazne varijable biramo po stupcima.

Postoji nekoliko načina odabira ulazne varijable.

Klasični pristup

Odabire se varijabla s najmanjim koeficijentom u Z-retku (u cilju maksimizacije). To je varijabla koja trenutno najviše može povećati vrijednost ciljne funkcije po jedinici. Ovaj pristup primjenjuje se od samih početaka Simplex metode (1947) i poznat je pod nazivom Most Negative Rule ili pravilo najnegativnijeg koeficijenta. Dantzigova knjiga Linear Programming and Extensions iz 1963. sadrži prvi cjeloviti opis ove metode odabira ulazne varijable.

U našem primjeru: Z-red ima `\(-3, -5\) \(\rightarrow\) odabiremo \(x_2\) jer je \(-5\) “najgori”, tj. najviše može poboljšati Z).

Prednost: Brzo vodi prema optimalnosti u mnogim slučajevima.

Nedostatak: Može izazvati cikluse (ponavljanje istih vrhova), osobito u degeneriranim slučajevima.

U danom primjeru (za prvu iteraciju):

\(Z−3x_1−5x_2=0\) ⇒ odabiremo \(x_2\).

Blandovo pravilo

Odabire se prva varijabla (po indeksu; redom) s negativnim koeficijentom u Z-retku (Bland, 1977).

Prednost: Jamči završetak algoritma (dokazano: sprječava beskonačne cikluse).

Nedostatak: Može uzrokovati sporiju konvergenciju, jer ne bira mogućnost najvećeg rasta.

U danom primjeru:

\(x_1\) je prvi s negativnim koeficijentom ⇒ odabiremo \(x_1\)

Najveći porast Z-a

Ovaj pristup poznat je i kao Dantzigov kriterij „largest coefficient rule”.

Za svaku varijablu s negativnim Z-koeficijentom izračuna se

\[\Delta Z = \left| \frac{c_j}{a_{ij}} \right|\]

gdje je \(a_{ij}\) pivot element za tu varijablu i-ti red. Odabire se onaj \(a_{ij}\) koja najviše povećava Z u toj iteraciji. Ili, drugim riječima - računaju se potencijalna povećanja Z vrijednosti za svaki stupac, a odabire se stupac s najvećim pozitivnim potencijalnim rastom Z-a.

Prednost: Može brže doći do rješenja.

Nedostatak: Računanje je složenije i zahtijeva više operacija.

U danom primjeru (za prvu iteraciju):

> library(tidyr)
> library(dplyr)
> 
> z_row <- c(x1 = -3, x2 = -5)
> tableau <- data.frame(
+   s1 = c(x1 = 2, x2 = 1),
+   s2 = c(x1 = 1, x2 = 2)
+ )
> 
> ratios <- sapply(names(z_row), function(var) {
+   best <- -Inf
+   for (col in colnames(tableau)) {
+     a_ij <- tableau[var, col]
+     if (a_ij > 0) {
+       delta_z <- abs(z_row[[var]] / a_ij)
+       best <- max(best, delta_z)
+     }
+   }
+   return(best)
+ })
> 
> best_delta <- names(which.max(ratios))
> paste("Odabrana varijabla:", best_delta)

## [1] "Odabrana varijabla: x2"

Nasumični izbor

engl. Randomized Pivot Rule

Nasumično se bira neka od varijabli s negativnim Z-koeficijentom.

Prednost: Može pomoći kod problema sa simetrijama ili degeneracijom.

Nedostatak: Nepredvidiv broj iteracija; ponekad vrlo neučinkovito.

U danom primjeru (za prvu iteraciju):

> set.seed(0)
> rand_var <- sample(names(z_row), 1)
> paste("Odabrana varijabla:", rand_var)

## [1] "Odabrana varijabla: x2"

Steepest edge rule

Odabire se ulazna varijabla tako da se maksimizira nagib (gradijent) funkcije cilja u smjeru promjene varijable (Goldfarb & Reid, 1977; Forrest & Goldfarb, 1992). Ova metoda traži varijablu koja najviše poboljšava funkciju cilja po jedinici euklidske udaljenosti u prostoru izvedivih rješenja.

Drugim riječima, mjeri se:

\[ \Delta Z = \frac{|c_j|}{\|a_j\|} \]

gdje je:

\(c_j\) – koeficijent u funkciji cilja
\(a_j\) – stupac Simplex tablice za tu varijablu
\(\|a_j\|\) – euklidska norma vektora (tj. kvadratni korijen sume kvadrata)

Prednost: Potencijalno najmanji broj iteracija.

Nedostatak: Zahtijeva izračun normi rubova i stoga je računski zahtjevnije.

Za dani primjer (za prvu iteraciju):

> z_row <- c(x1 = -3, x2 = -5)
> 
> tableau <- data.frame(
+   s1 = c(x1 = 2, x2 = 1),
+   s2 = c(x1 = 1, x2 = 2)
+ )
> 
> # Računanje normi stupaca
> norme <- sapply(names(z_row), function(var) {
+   sqrt(sum(tableau[var, ]^2))
+ })
> 
> # Računanje ΔZ
> delta_z <- abs(z_row) / norme
> 
> # Tablica rezultata – bez zaokruživanja za usporedbu
> steepest_df <- data.frame(
+   Varijabla = names(z_row),
+   abs_cj = abs(z_row),
+   norma_aj = norme,
+   delta_Z = delta_z
+ )
> 
> # Odabir varijable s najvećim ΔZ
> max_delta <- max(steepest_df$delta_Z)
> steepest_df$Odabir <- ifelse(steepest_df$delta_Z == max_delta, "DA", "")
> 
> # Prikaz – s formatiranjem
> steepest_df$abs_cj <- round(steepest_df$abs_cj, 3)
> steepest_df$norma_aj <- round(steepest_df$norma_aj, 3)
> steepest_df$delta_Z <- round(steepest_df$delta_Z, 3)
> 
> names(steepest_df) <- c("Varijabla", "|c_j|", "||a_j||", "ΔZ = |c_j| / ||a_j||", "Odabir")
> 
> knitr::kable(steepest_df, caption = "Odabir ulazne varijable prema Steepest Edge pravilu")

Odabir ulazne varijable prema Steepest Edge pravilu
	Varijabla	\|c_j\|	\|\|a_j\|\|	ΔZ = \|c_j\| / \|\|a_j\|\|	Odabir
x1	x1	3	2.236	1.342
x2	x2	5	2.236	2.236	DA

Devex metoda

Ova metoda prvi je put objavljena 1973. kao “Pivot selection methods of the Devex LP code” autorice Paule M. J. Harris, koja formalno uvodi heuristiku Devex.

praktična aproksimacija steepest edge metode
Heuristika koja aproksimira steepest edge bez punih normi.
Koristi težine za varijable i ažurira ih tijekom iteracija.
Ona ne računa točnu normu svakog ruba, već koristi jednostavniji kriterij temeljen na iskustvu rješavanja problema.

Prednosti: – Bolje performanse od najnegativnijeg koeficijenta u praksi. – Koriste ga neki visokooptimizirani LP solveri.

Nedostatak: - Kao heuristika, nema jamstva da je uvijek bolja od drugih pristupa. - Može loše reagirati na degenerirane ili visoko korelirane probleme.

Za dani primjer (za prvu iteraciju):

> z_row <- c(x1 = -3, x2 = -5)
> 
> weights <- c(x1 = 1, x2 = 1)
> 
> # Računamo Devex omjer: |c_j| / sqrt(w_j)
> delta_devex <- abs(z_row) / sqrt(weights)
> 
> # Napravimo tablicu bez zaokruživanja za usporedbu
> devex_df <- data.frame(
+   Varijabla = names(z_row),
+   abs_cj = abs(z_row),
+   wj = weights,
+   delta_Z = delta_devex
+ )
> 
> # Odredimo maksimum i označimo odabir
> max_val <- max(devex_df$delta_Z)
> devex_df$Odabir <- ifelse(devex_df$delta_Z == max_val, "DA", "")
> 
> # Formatiramo brojeve za prikaz
> devex_df$`|c_j|` <- round(devex_df$abs_cj, 3)
> devex_df$`w_j` <- devex_df$wj
> devex_df$`ΔZ = |c_j| / sqrt(w_j)` <- round(devex_df$delta_Z, 3)
> 
> # Uklanjamo pomoćne stupce
> devex_df <- devex_df[, c("Varijabla", "|c_j|", "w_j", "ΔZ = |c_j| / sqrt(w_j)", "Odabir")]
> 
> # Prikaz tablice
> knitr::kable(devex_df, caption = "Odabir ulazne varijable prema Devex metodi")

Odabir ulazne varijable prema Devex metodi
	Varijabla	\|c_j\|	w_j	ΔZ = \|c_j\| / sqrt(w_j)	Odabir
x1	x1	3	1	3
x2	x2	5	1	5	DA

Kad odabrati pojedini pristup?

Pristup	Opis	Prednosti	Nedostaci
Klasičan pristup	Odabire se varijabla s najnegativnijim koeficijentom u Z-retku	Brzo vodi do rješenja Najprikladniji za male modele i ručni izračun	Može izazvati cikluse u degeneriranim slučajevima
Blandovo pravilo	Odabire se prva (po indeksu) negativna varijabla	Sprječava beskonačne cikluse	Može ignorirati najbržii rast; sporije konvergira
Najveći porast \(Z\)-a	Odabire se varijabla s najvećim omjerom \(\|c_j / a_{ij}\|\)	Fokusira se na varijablu koja najviše poboljšava \(Z\)	Zahtijeva dodatne izračune za sve redove i varijable
Nasumični izbor	Nasumični izbor među negativnim \(Z\)-koeficijentima	Jednostavan, pomaže u simetričnim slučajevima	Nepredvidivo ponašanje; često neučinkovit
Steepest edge rule	Odabir na temelju najvećeg gradijenta prema normi \(\\|a_j\\|\)	Obično daje najmanji broj iteracija	Zahtijeva normu svakog stupca; računski intenzivno
Devex metoda	Heuristika s težinama umjesto normi (aproksimacija steepest edge)	Brži od steepest edge, dobar kompromis	Nije precizan optimizacijski kriterij, nego heuristički pristup

Odabir izlazne varijable

Cilj je odrediti koja bazična varijabla će napustiti bazu – odnosno, koji red postaje pivot redak. Dakle, izlazne varijable biramo iz redaka.

Koja varijabla može izaći iz baze?

Kada ulazna varijabla uđe u bazu (retke), jedna od postojećih mora izaći
Iz baze može izaći samo neka od trenutnih bazičnih varijabli

Pri odabiru izlaznih varijabli također postoji više pristupa.

Pravilo minimuma omjera

Ovo je klasično pravilo koje se primjenjuje od samih početaka Simplex metode (Dantzig, 1947) — formalno objavljeno unutar Dantzigovog rada o simplex algoritmu iz 1947. godine.

Za sve redove (ograničenja), računamo:

\[\frac{b_i}{a_{ij}}\]

gdje je:

\(b_i\) ili RHS, vrijednost desne strane
\({a_{ij}}\) je koeficijent ulazne varijable u tom retku

Odabire se red s najmanjim pozitivnim omjerom - to jamči očuvanje ne-negativnosti rješenja.

Ovaj pristup se ne može koristiti ako su svi \(a_{ij}≤0\), jer bi to značilo da možemo beskonačno povećavati ulaznu varijablu bez kršenja ograničenja \(\rightarrow\) problem nema ograničeno rješenje (unbounded).

Prednost: Jednostavno, brzo, štiti od negativnih rješenja

Nedostatak: U degeneriranim slučajevima može dovesti do ciklusa (beskonačnog ponavljanja istih rješenja)

Na našem primjeru,

BV	x₁	x₂	s₁	s₂	RHS
s₁	2	1	1	0	100
s₂	1	2	0	1	80
Z	-3	-5	0	0	0

uz pretpostavku da smo odabrali \(x_2\) kao ulaznu varijablu:

> # Simplex tablica (za odabir izlazne varijable)
> tableau <- data.frame(
+   x1 = c(2, 1),
+   x2 = c(1, 2),
+   s1 = c(1, 0),
+   s2 = c(0, 1),
+   RHS = c(100, 80)
+ )
> rownames(tableau) <- c("s1", "s2")
> 
> # Pretpostavimo da je ulazna varijabla x2
> entering <- "x2"
> 
> # Povlačimo stupac x2 i računamo omjere RHS / x2
> aij <- tableau[[entering]]
> rhs <- tableau$RHS
> 
> ratio <- ifelse(aij > 0, rhs / aij, Inf)
> 
> # Tablica s omjerima
> out_table <- data.frame(
+   BV = rownames(tableau),
+   `a_ij` = aij,
+   RHS = rhs,
+   `RHS / a_ij` = round(ratio, 3)
+ )
> 
> # Odabir minimalnog pozitivnog omjera
> min_row <- which.min(ratio)
> out_table$Odabir <- ""
> out_table$Odabir[min_row] <- "DA"
> 
> knitr::kable(out_table, caption = paste("Odabir izlazne varijable za ulaznu varijablu:", entering))

Odabir izlazne varijable za ulaznu varijablu: x2
BV	a_ij	RHS	RHS…a_ij	Odabir
s1	1	100	100
s2	2	80	40	DA

Blandovo pravilo

Ako postoji više redova s istim minimalnim omjerom, odabire se onaj s varijablom koja ima najmanji indeks (redoslijedno/ abecedno), s ciljem sprječavanja ciklusa i ostvarivanja konačnosti algoritma (Bland, 1977).

Prednost: sprječava cikluse

Nedostatak: može biti sporije (jer ne bira najagresivniji korak).

U ovom primjeru pretpostavljamo da je kao ulazna varijabla odabrana x2, i da početna Simplex tablica izgleda ovako:

BV	x₁	x₂	s₁	s₂	RHS
s₁	2	1	1	0	100
s₂	1	2	0	1	80
Z	-3	-5	0	0	0

> tableau <- data.frame(
+   x1 = c(2, 1),
+   x2 = c(1, 2),
+   s1 = c(1, 0),
+   s2 = c(0, 1),
+   RHS = c(100, 80)
+ )
> rownames(tableau) <- c("s1", "s2")
> entering <- "x2"
> 
> # Dohvat koeficijenata i desne strane
> aij <- tableau[[entering]]
> rhs <- tableau$RHS
> ratio <- ifelse(aij > 0, rhs / aij, Inf)
> 
> # Kandidati s minimalnim (pozitivnim) omjerom
> min_ratio <- min(ratio)
> bland_candidates <- which(ratio == min_ratio)
> 
> # Dohvati imena redova kandidata i uzmi abecedno prvo ime
> candidate_names <- rownames(tableau)[bland_candidates]
> bland_row_name <- sort(candidate_names)[1]
> 
> # Dohvati redni broj tog reda
> bland_row <- which(rownames(tableau) == bland_row_name)
> bland_row

## [1] 2

Dakle \(s_2\) će biti izlazna varijabla.

Tie-breaking rules

Metoda se koristi kada pravilo minimuma omjera daje isti omjer u više redaka. Da se izbjegne ciklus, primjenjuju se dodatna pravila, npr. odabir one s najmanjim indeksom. Prethodno razmatrano Blandovo pravilo je najpoznatije iz ove skupine.

Mogućnosti su:

odabir reda s najmanjom vrijednosti \(a_{ij}\)
ili s najvećim povećanjem Z
ili nasumično (ako sve ostalo zakaže).

Dakle, pristup Tie-breaking rules koristi se kad više redova ima isti (minimalni) omjer RHS / a₍ᵢⱼ₎, i tada treba dodatno pravilo za odabir pivot reda (izlazne varijable).

U danom primjeru to nije slučaj, ali će se prikazati kako bi postupak izgledao:

> tableau <- data.frame(
+   x1 = c(2, 1),
+   x2 = c(1, 2),
+   s1 = c(1, 0),
+   s2 = c(0, 1),
+   RHS = c(100, 80)
+ )
> rownames(tableau) <- c("s1", "s2")
> 
> entering <- "x2"
> 
> # Koeficijenti u stupcu i RHS
> aij <- tableau[[entering]]
> rhs <- tableau$RHS
> ratio <- ifelse(aij > 0, rhs / aij, Inf)
> 
> # Kandidati s jednakim minimalnim omjerom
> min_ratio <- min(ratio)
> tie_candidates <- which(ratio == min_ratio)
> 
> # Iz njih odaberi onaj s najmanjim a_ij (manji korak)
> tie_break_row <- tie_candidates[which.min(aij[tie_candidates])]
> tie_break_row

## [1] 2

Dakle \(s_2\) će biti izlazna varijabla.

Steepest-edge

U nekim varijantama se i za izlaznu varijablu koristi geometrijski pristup – bira se red koji najviše smanjuje normu vektora osnovnih varijabli, tj. “najstrmiji pad” (Forrest & Goldfarb, 1992).

Treba napomenuti da je to rjeđe korišteni pristup i da se češće Steepest-edge pravilo koristi za izbor ulaznih varijabli, u kombinaciji s pravilom minimuma omjera za izlazne varijable. No, postoji varijanta za odabir izlaznih varijabli, koja je primjerena za korištenje u određenim situacijama.

Prednost: brže konvergira u velikim problemima

Nedostatak: složenije za izračun.

Za dani primjer, to izgleda ovako:

> tableau <- data.frame(
+   x1 = c(2, 1),
+   x2 = c(1, 2),
+   s1 = c(1, 0),
+   s2 = c(0, 1),
+   RHS = c(100, 80)
+ )
> rownames(tableau) <- c("s1", "s2")
> 
> entering <- "x2"
> 
> # Koeficijenti i RHS
> aij <- tableau[[entering]]
> rhs <- tableau$RHS
> 
> # Računamo euklidsku normu svakog reda
> row_norms <- apply(tableau[, c("x1", "x2", "s1", "s2")], 1, function(row) sqrt(sum(row^2)))
> 
> # Računamo ΔZ_i = RHS / norma retka (samo ako a_ij > 0)
> delta_steep <- ifelse(aij > 0, rhs / row_norms, Inf)
> 
> # Odabir reda s najmanjim ΔZ_i
> steep_row <- which.min(delta_steep)
> steep_row

## [1] 2

Ponovo, \(s_2\) treba biti izlazna varijabla.

Devex heuristika

Nakon što je Devex izabrao ulaznu varijablu, iste težine se koriste i pri izboru izlazne, preferirajući varijable čiji omjeri i težine sugeriraju značajan pomak (Harris, 1973).

> tableau <- data.frame(
+   x1 = c(2, 1),
+   x2 = c(1, 2),
+   s1 = c(1, 0),
+   s2 = c(0, 1),
+   RHS = c(100, 80)
+ )
> rownames(tableau) <- c("s1", "s2")
> 
> entering <- "x2"
> 
> # Početne težine za redove (bazične varijable) – sve 1 u prvoj iteraciji
> devex_weights <- c(s1 = 1, s2 = 1)
> 
> # Koeficijenti ulazne varijable
> aij <- tableau[[entering]]
> rhs <- tableau$RHS
> 
> # Devex omjer: RHS / sqrt(weight), samo ako a_ij > 0
> delta_devex <- ifelse(aij > 0, rhs / sqrt(devex_weights), Inf)
> 
> # Odabir reda s najmanjim Devex omjerom
> devex_row <- which.min(delta_devex)
> 
> # Tablica za prikaz
> devex_df <- data.frame(
+   BV = rownames(tableau),
+   a_ij = aij,
+   RHS = rhs,
+   w = devex_weights,
+   ΔZ = round(delta_devex, 3),
+   Devex = ifelse(1:2 == devex_row, "Da", "")
+ )
> 
> knitr::kable(devex_df, caption = "Odabir izlazne varijable prema Devex heuristici")

Odabir izlazne varijable prema Devex heuristici
	BV	a_ij	RHS	w	ΔZ	Devex
s1	s1	1	100	1	100
s2	s2	2	80	1	80	Da

Dualni Simplex

Dualni Simplex koristi se kad rješenje nije primarno dopustivo (neki RHS negativni), ali je dualno dopustivo (Lemke, 1954). Ovdje se radi o promjeni perspektive: dualni Simplex promatra rješenje iz dualne perspektive – izlazna varijabla se bira prema najnegativnijem RHS, a ulazna prema minimalnom omjeru dualnih koeficijenata. Ovdje neće biti Z-redak taj koji pokazuje koja varijabla ulazi, nego se traži redak s negativnim RHS (\(b_i\)) – izlazna varijabla u dualnom smislu. Ovo vrijedi u slučaju kad:

se promijeni RHS nakon promjene modela,
ili imamo nevaljano početno rješenje (negativni RHS).

Dualni Simplex je metoda u kojoj je rješenje dualno dopustivo (Z-redak je “uredan”, tj. svi koeficijenti funkcije cilja su ≥ 0), ali ne primarno (neki RHS < 0). Ulazi se u bazu tako da se poboljša primarna dopustivost.

Dakle, u našem primjeru nema smisla koristiti ovu metodu, pa ćemo za potrebe primjera promijeniti drugo ograničenje u \(-80\).

> # Simulirani tableau s negativnim RHS
> tableau <- data.frame(
+   x1 = c(2, 1),
+   x2 = c(1, 2),
+   s1 = c(1, 0),
+   s2 = c(0, 1),
+   RHS = c(100, -80)
+ )
> rownames(tableau) <- c("s1", "s2")
> 
> z_row <- c(x1 = 0, x2 = 0, s1 = 0, s2 = 0)
> 
> # 1. Odabir izlazne varijable $\rightarrow$ najnegativniji RHS
> rhs <- tableau$RHS
> out_idx <- which.min(rhs)
> pivot_row <- tableau[out_idx, 1:4]
> 
> # 2. Odabir ulazne varijable $\rightarrow$ najmanji (cj / aij) za aij > 0
> ratios <- sapply(names(pivot_row), function(var) {
+   aij <- pivot_row[[var]]
+   cj <- z_row[[var]]
+   if (aij > 0) return(cj / aij) else return(Inf)
+ })
> 
> in_var <- names(which.min(ratios))
> 
> # Rezultat:
> cat("Izlazna varijabla:", rownames(tableau)[out_idx], "\n")

## Izlazna varijabla: s2

> cat("Ulazna varijabla:", in_var, "\n")

## Ulazna varijabla: x1

Kad odabrati pojedini pristup?

Pristup	Kada koristiti
Minimum omjera	Standardno
Blandovo pravilo	Za izbjegavanje ciklusa
Tie-breaking pravila	Kod veza
Steepest-edge	Kod velikih modela
Devex	Ako se već koristi Devex za ulaznu
Dualni Simplex	Kod negativnih RHS vrijednosti

Update tablice

U ovom koraku uvijek se koristi varijanta Gauss-Jordanove eliminacije, ali postoje nekoliko pristupa u načinu implementacije i računanja, ovisno o:

vrsti Simplex metode (standardni, dualni, dvostupanjski),
veličini problema (ručno vs. računalno),
numeričkoj stabilnosti i učinkovitosti u računalnim rješenjima.

Standardni pristup (ručna verzija) – Gauss-Jordanova eliminacija

Pivot element se postavlja na 1 (dijeljenjem cijelog retka)
Svi ostali elementi u pivot stupcu postaju 0 (linearnim kombinacijama s pivot retkom)

Ovo se koristi:

u nastavi
u ručnom rješavanju tablica
u didaktičkom smislu (omogućuje razumijevanje osnovne ideje procedura)

Na danom primjeru, uz \(x_2\) kao ulaznu varijablu i \(s_2\) kao izlaznu varijablu, provedba funkcionira na sljedeći način:

BV	x₁	x₂	s₁	s₂	RHS
s₁	2	1	1	0	100
s₂	1	2	0	1	80
Z	\(-3\)	\(-5\)	0	0	0

Postupak kreće s identifikacijom pivot elementa, koji je ovdje \(a_{22}\), obojan ljubičasto. \(x_2\) ulazi u bazične varijable, a \(s_2\) izlazi, tj. \(x_2\) dolazi na mjesto \(s_2\). Potom se cijeli redak normalizira na način da \(a_{22}\) poprimi vrijednost 1. To dobivamo tako što ćemo \(a_{22}\) podijeliti sa samim sobom, a potom i sve ostale elemente u tom retku dijelimo s \(a_{22}\).

BV	x₁	x₂	s₁	s₂	RHS
s₁	2	1	1	0	100
x₂	\(\frac{1}{2}\)	\(\frac{2}{2}=1\)	\(\frac{0}{2}=0\)	\(\frac{1}{2}\)	\(\frac{80}{2}=40\)
Z	\(-3\)	\(-5\)	0	0	0

Sljedeći je korak sve ostale koeficijente u stupcu ulazne varijable \(x_2\) postaviti na nulu, a da vrijednost u pivot ćeliji ostane 1. Pa ćemo to za prvi redak provesti tako što ćemo vrijednosti u prvom retku umanjiti za vrijednosti u drugom retku bez drugih intervencija, jer je \(a_{12}=a_{22}\) (pa je dovoljno izravno oduzeti vrijednosti kako bismo dobili nule u stupcu \(x_2\), uz dosljednu primjenu na ostale vrijednosti u retku). Za treći redak, potrebno je eliminirati \(a_{32}=-5\), pa ćemo drugi redak prvo pomnožiti s 5, a potom ga pridodati trećem retku, slijedeći istu logiku.

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(2-\frac{1}{2}\)	\(1 - 1\)	\(1 - 0\)	\(0 - \frac{1}{2}\)	\(100-40\)
x₂	\(\frac{1}{2}\)	\(1\)	\(0\)	\(\frac{1}{2}\)	\(40\)
Z	\(-3 + 5 \cdot \frac{1}{2}\)	\(-5 + 5 \cdot 1\)	\(0 + 5 \cdot 0\)	\(0 + 5 \cdot \frac{1}{2}\)	\(0 + 5 \cdot 40\)

Ove intervencije rezultiraju sljedećim izgledom tablice:

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(\frac{3}{2}\)	0	1	\(-\frac{1}{2}\)	60
x₂	\(\frac{1}{2}\)	1	0	\(\frac{1}{2}\)	40
Z	\(-\frac{1}{2}\)	0	0	\(\frac{5}{2}\)	200

U ovom koraku Simplex metode dobili smo novo bazično rješenje, što znači da smo došli do novog vrha poligona izvedivih rješenja. U 2D prostoru, svaki korak Simplex metode odgovara skoku s jednog vrha na susjedni.

Rješenja možemo iščitati samo za one varijable uz koje stoji 1, a u ostatku stupca stoje 0-e (bazične varijable), pa u ovoj iteraciji to daje rješenje \((0, 40)\). Ako te vrijednosti uvrstimo u funkciju cilja, to daje \(3 \cdot 0 + 5 \cdot 40 = 200\) (vrijednost koju uočavamo u trećem retku, kao trenutnu vrijednost funkcije cilja).

Ovo je bila prva iteracija. Nastavljamo s postupkom, tako što ćemo, prema klasičnom pristupu odabrati \(x_1\) kao ulaznu varijablu, a prema pravilu minimalnih omjera (\(60 : \frac{3}{2}= \frac{\frac{60}{1}}{\frac{3}{2}}=40\), \(40 : \frac{1}{2}=\frac{\frac{40}{1}}{\frac{1}{2}}=80\)) odabiremo \(s_1\) kao izlaznu varijablu.

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(\frac{3}{2}\)	0	1	\(-\frac{1}{2}\)	60
x₂	\(\frac{1}{2}\)	1	0	\(\frac{1}{2}\)	40
Z	\(-\frac{1}{2}\)	0	0	\(\frac{5}{2}\)	200

Pivot element je \(a_{11}=\frac{3}{2}\) i započinjemo tako što ćemo sve elemente u tom retku podijeliti s \(\frac{3}{2}\) (odnosno pomnožiti s \(\frac{2}{3}\)).

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(\frac{3}{2} \cdot \frac{2}{3}=1\)	\(0 \cdot \frac{2}{3}=0\)	\(1 \cdot \frac{2}{3} = \frac{2}{3}\)	\(-\frac{1}{2} \cdot \frac{2}{3} = -\frac{1}{3}\)	\(60 \cdot \frac{2}{3}=40\)
x₂	\(\frac{1}{2}\)	1	0	\(\frac{1}{2}\)	40
Z	\(-\frac{1}{2}\)	0	0	\(\frac{5}{2}\)	200

Tablica sad izgleda ovako:

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(1\)	\(0\)	\(\frac{2}{3}\)	\(-\frac{1}{3}\)	\(40\)
x₂	\(\frac{1}{2}\)	1	0	\(\frac{1}{2}\)	40
Z	\(-\frac{1}{2}\)	0	0	\(\frac{5}{2}\)	200

Preostaje još koeficijente \(\frac{1}{2}\) i \(-\frac{1}{2}\) postaviti na 0. To ćemo učiniti tako da ćemo prvi redak pomnožiti s koeficijentom ispred \(x_1\) u retku (tj. \(\frac{1}{2}\)) i oduzeti ga od drugog retka te pribrojiti trećem retku.

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(1\)	\(0\)	\(\frac{2}{3}\)	\(-\frac{1}{3}\)	\(40\)
x₂	\(\frac{1}{2}- 1 \cdot\frac{1}{2}=0\)	\(1- 0 \cdot\frac{1}{2}=1\)	\(0- \frac{2}{3} \cdot\frac{1}{2}=-\frac{1}{3}\)	\(\frac{1}{2}- \left(- \frac{1}{3} \right) \cdot\frac{1}{2} = \frac{2}{3}\)	\(40 - 40 \cdot\frac{1}{2}=20\)
Z	\(-\frac{1}{2}+ 1 \cdot\frac{1}{2}=0\)	\(0 + 0 \cdot\frac{1}{2}=0\)	\(0+ \frac{2}{3} \cdot\frac{1}{2}=\frac{1}{3}\)	\(\frac{5}{2} + \left(-\frac{1}{3} \right) \cdot \frac{1}{2}=\frac{7}{3}\)	\(200+20=220\)

Ove transformacije rezultiraju sljedećom tablicom:

BV	x₁	x₂	s₁	s₂	RHS
s₁	\(1\)	\(0\)	\(\frac{2}{3}\)	\(-\frac{1}{3}\)	\(40\)
x₂	\(0\)	\(1\)	\(-\frac{1}{3}\)	\(\frac{2}{3}\)	\(20\)
Z	\(0\)	\(0\)	\(\frac{1}{3}\)	\(\frac{7}{3}\)	\(220\)

U točci \((40, 20)\), vrijednost funkcije cilja je \(40 \cdot 3 + 20 \cdot 5 = 220\), što ukazuje na povećanje u vrijednosti funckije cilja uslijed druge iteracije. Grafički prikaz pomaže vizualno pratiti Simplex metodu, ali je primjenjiv isključivo na probleme s dvije varijable odluke.

U Simplex metodi (maksimizacija), promatramo \(Z-redak\) (ili redak funkcije cilja) i:

ako su svi koeficijenti varijabli u Z-redu ≥ 0, onda:
- nije moguće dodatno povećati funkciju cilja,
- ne postoji negativan “put” prema većoj vrijednosti,
- pa je trenutno rješenje optimalno.

U tablici kreiranoj drugom iteracijom, ovi su uvjeti zadovoljeni, pa zaključujemo da je to optimalno rješenje. Dakle, postupak je završen, jer ne postoji više niti jedna negativna vrijednost u Z-redu (uz varijable koje mogu “ući”), što ujedno znači da više nema mogućnosti pomaka prema višim vrijednostima funkcije cilja.

Računalni pristupi (optimizirani)

U praktičnim računalnim rješenjima koriste se:

Revidirani Simplex (Revised Simplex)
Product Form of the Inverse (PFI)
LU dekompozicija bazne matrice

Revidirani Simplex (Revised Simplex)

Ne ažurira cijelu tablicu
Radi sa baznom matricom B i njezinom inverzom \(B^{−1}\)
Ažurira samo potrebne vektore
Mnogo brži i stabilniji za velike probleme

Za razliku od klasičnog Simplex algoritma, Revidirani Simplex izbjegava ažuriranje cijele tablice. Ovaj pristup koristi samo bazne varijable i matricu \(B^{−1}\), čime štedi memoriju i računsko vrijeme — osobito kod velikih problema.

Za dani problem, u nastavku je prikazan postupak koristeći ovaj pristup. Započinjemo s pripremom i unosom vektora i matrica. Krećemo sa standardnim zapisom linearnog programa pri rješavanju u R-u, koji se oslanja na opći zapis linearnog progama u obliku:

\[ \text{max } z=c^Tx\]

\[Ax=b\]

\[x≥0\]

gdje su:

\(z\) - funkcija cilja,
\(c^T\) - transponirani vektor koeficijenata funkcije cilja,
\(A\) - matrica dimenzija \(n×m\) i predstavlja koeficijente lijeve strane ograničenja,
\(b\) - vektor desne strane ograničenja,
\(x\) - vektor varijabli odluka.

Također, dodajemo slack varijable \(s_1\) i \(s_2\) kako bismo koristili kanonski zapis.

> library(MASS)
> library(knitr)
> library(kableExtra)
> 
> cvec <- c(3, 5, 0, 0) # koeficijenti uz varijable u funkciji cilja, vektor c prije transponiranja
> 
> A <- matrix(c(2,1,1,0,
+               1,2,0,1), nrow=2, byrow=TRUE) # matrica koeficijenata na lijevoj strani ograničenja
> colnames(A) <- c("x1", "x2", "s1", "s2")
> 
> b <- c(100, 80) # desna strana ograničenja

Iz početnih postavki, izvlačimo dijelove potrebne za daljnje izračune:

> B <- A[, 3:4]  # matrica koeficijenata uz s1 i s2 u ograničenjima (lijeva strana ograničenja)
> B

##      s1 s2
## [1,]  1  0
## [2,]  0  1

> N <- A[, 1:2]  # matrica koeficijenata uz x1 i x2 u ograničenjima (lijeva strana ograničenja)
> N

##      x1 x2
## [1,]  2  1
## [2,]  1  2

> cb <- cvec[3:4]  # koeficijenti uz s1 i s2 u funkciji cilja
> cb

## [1] 0 0

> cN <- cvec[1:2]  # koeficijenti uz x1 i x2 u funkciji cilja
> cN

## [1] 3 5

Inverz matrice \(B\), \(B^{−1}\) je:

> B_inv <- solve(B)
> B_inv

##    [,1] [,2]
## s1    1    0
## s2    0    1

Digresija. solve() će dati rješenje samo ako je matrica kvadratna i regularna (tj. \(det(A) ≠ 0\))

Ako matrica nije invertibilna, dobivamo poruku: “Lapack routine dgesv: system is exactly singular”

Ako je u pitanju takva situacija ili je potreban numerički stabilniji ili poseban pristup (npr. za velike sustave, rjeđe situacije, pseudoinverz), tada u ozbir dolazi paket MASS i funkcija ginv(), npr:

> library(MASS)
> ginv(A)

Osim toga, u istom paketu dostupne su funkcije:

qr.solve(A, b) — koristi QR dekompoziciju za rješavanje sustava
chol2inv(chol(A)) — koristi Cholesky dekompoziciju (za simetrične pozitivno-definirane matrice)

U 1. koraku transformiraju se vektori \(b_i\) (desna strana ograničenja) množenjem s \(B^{−1}\). Potom se dobivena vrijednost množi s transponiranim vektorom koeficijenata uz bazične varijable i tako se dobiva nova vrijednost funkcije cilja (\(Z\)). U sljedećem koraku izračunava se pomoćni vektor \(z_j\), tako što se transponirani vektor uz bazične varijable množi s \(B^{−1}\) i matricom koeficijenata uz varijable odluke u ograničenjima (lijeva strana ograničenja). \(z_j\) je zapravo marginalni doprinos varijable \(x_j\) izračunat preko trenutne baze, tj. koliko bi se funkcija cilja promijenila ako bismo tu varijablu uključili u bazu.

Razlika vektora koeficijenata uz varijable odluke u funkciji cilja (\(c_j\)) i vektora \(z_j\) daje izraz \(c_j−z_j\). Te se vrijednosti nazivaju još i reduced cost ili reducirani trošak i označavaju:

koliko bi se vrijednost funkcije cilja promijenila kada bi varijabla izvan baze \(x_j\) postala aktivna
odnosi se na promjenu koeficijenta u funkciji cilja (\(c_j\))
intuitivna interpretacija jest da nam taj izraz govori koliko košta ubacivanje te varijable u bazu — tj. je li korisna za funkciju cilja.

Naziv reduced cost dolazi iz činjenice da je ukupni trošak (ili doprinos) \(c_j\) reduciran za ono što baza već implicitno doprinosi funkciji cilja. Drugim riječima:

Ako je varijabla izvan baze, njezina vrijednost je trenutno 0.
Ako je \(c_j−z_j>0\) kod maksimizacije, tada uključivanje te varijable umanjuje funkciju cilja \(\rightarrow\) nije korisna.
Ako je\(c_j−z_j<0\), tada varijabla može povećati funkciju cilja ako uđe u bazu \(\rightarrow\) dobar kandidat za ulazak.
Ako je \(c_j−z_j=0\), uključivanje varijable ne mijenja funkciju cilja \(\rightarrow\) imamo višestruka optimalna rješenja (degeneracija).

> xB <- B_inv %*% b
> xB

##    [,1]
## s1  100
## s2   80

> Z <- t(cb) %*% xB
> Z

##      [,1]
## [1,]    0

> zj <- t(cb) %*% B_inv %*% N
> zj

##      x1 x2
## [1,]  0  0

> cj_zj <- cN - zj
> cj_zj

##      x1 x2
## [1,]  3  5

> cat("Vrijednosti varijabli nakon 1. iteracije:\n")

## Vrijednosti varijabli nakon 1. iteracije:

> solution1 <- rep(0, 4)
> names(solution1) <- colnames(A)
> solution1[c("s1", "s2")] <- round(xB, 2)
> print(solution1)

##  x1  x2  s1  s2 
##   0   0 100  80

> cat(sprintf("Vrijednost funkcije cilja: Z = %.2f\n", Z))

## Vrijednost funkcije cilja: Z = 0.00

Druga iteracija započinje odabirom ulazne i izlazne varijable. Vrijednosti razlike \(c_j−z_j\) predstavljaju koliko bi iznosilo povećanje funkcije cilja kad bismo uveli pojedinu varijablu u bazu. Ako su sve te vrijednosti \(≤ 0\), tada više nema prostora za poboljšanje i postignuto je optimalno rješenje. Najveća vrijednost \(cj - zj\) pojavljuje se za \(x_2 = 5\), pa \(x_2\) ulazi u bazu.

Računamo smjer promjene \(d = B^{-1}a_j\), gdje je \(a_j\) stupac matrice A za \(x_2\).

> a_j <- A[, 2]  # stupac za x2
> d <- B_inv %*% a_j
> d

##    [,1]
## s1    1
## s2    2

> ratio <- ifelse(d > 0, xB / d, Inf)
> leaving_index <- which.min(ratio)
> leaving_var <- c("s1", "s2")[leaving_index]
> 
> list("d" = round(d, 2), "theta" = round(ratio, 2), "Izlazna varijabla" = leaving_var)

## $d
##    [,1]
## s1    1
## s2    2
## 
## $theta
##    [,1]
## s1  100
## s2   40
## 
## $`Izlazna varijabla`
## [1] "s2"

Sljedeći korak u drugoj iteraciji je ažuriranje baze.

> # x2 ulazi, s2 izlazi
> B <- A[, c(3,2)]  # s1 i x2
> cb <- cvec[c(3,2)]
> B_inv <- solve(B)
> xB <- B_inv %*% b
> Z <- t(cb) %*% xB
> N <- A[, c(1,4)]
> cN <- cvec[c(1,4)]
> zj <- t(cb) %*% B_inv %*% N
> cj_zj <- cN - zj
> 
> result2 <- data.frame(
+   "Bazicne varijable" = c("s1", "x2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> result2

##    Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..x1
## s1                s1             60                        0.5
## x2                x2             40                        0.5
##    Reduced.costs..cj...zj..s2
## s1                       -2.5
## x2                       -2.5

> cat("Vrijednosti varijabli nakon 2. iteracije:\n")

## Vrijednosti varijabli nakon 2. iteracije:

> solution2 <- rep(0, 4)
> names(solution2) <- colnames(A)
> solution2[c("s1", "x2")] <- round(xB, 2)
> print(solution2)

## x1 x2 s1 s2 
##  0 40 60  0

> cat(sprintf("Vrijednost funkcije cilja: Z = %.2f\n", Z))

## Vrijednost funkcije cilja: Z = 200.00

Ovdje je dosegnut drugi vrh poligona, no nije postignuto optimalno rješenje. To znači da treba provesti još jednu iteraciju. Ponovo se odabiru ulazna i izlazna varijabla.

> a_j <- A[, 1]  # stupac za x1
> d <- B_inv %*% a_j
> 
> ratio <- ifelse(d > 0, xB / d, Inf)
> leaving_index <- which.min(ratio)
> leaving_var <- c("s1", "x2")[leaving_index]
> 
> list("d" = round(d, 2), "theta" = round(ratio, 2), "Izlazna varijabla" = leaving_var)

## $d
##    [,1]
## s1  1.5
## x2  0.5
## 
## $theta
##    [,1]
## s1   40
## x2   80
## 
## $`Izlazna varijabla`
## [1] "s1"

> B <- A[, c(1,2)]  # x1 i x2
> cb <- cvec[c(1,2)]
> B_inv <- solve(B)
> xB <- B_inv %*% b
> Z <- t(cb) %*% xB
> N <- A[, c(3,4)]
> cN <- cvec[c(3,4)]
> zj <- t(cb) %*% B_inv %*% N
> cj_zj <- cN - zj
> 
> result3 <- data.frame(
+   "Bazicne varijable" = c("x1", "x2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> result3

##    Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..s1
## x1                x1             40                      -0.33
## x2                x2             20                      -0.33
##    Reduced.costs..cj...zj..s2
## x1                      -2.33
## x2                      -2.33

> cat("Vrijednosti varijabli nakon 3. iteracije:\n")

## Vrijednosti varijabli nakon 3. iteracije:

> solution <- rep(0, 4)
> names(solution) <- colnames(A)
> solution[c("x1", "x2")] <- round(xB, 2)
> print(solution)

## x1 x2 s1 s2 
## 40 20  0  0

> cat(sprintf("Vrijednost funkcije cilja: Z = %.2f\n", Z))

## Vrijednost funkcije cilja: Z = 220.00

Budući da su svi \(cj - zj \leq 0\), optimalno rješenje je postignuto:

\(x_1 = 40\)
\(x_2 = 20\)
\(Z = 220\)

Time završavamo postupak Revidiranog Simplex algoritma.

Računanje inverza matrice je računski zahtjevno i numerički nestabilno kod velikih matrica. Zbog toga se koriste tehnike poput Product Form of the Inverse (PFI) ili LU dekompozicije.

Product Form of the Inverse (PFI)

Umjesto ponovnog računanja \(B^{−1}\) nakon svake iteracije, koristi produkt jednostavnih transformacija
Brži za neke specifične implementacije
Numerički stabilniji kod velikih problema gdje se više puta radi zamjena u bazi
Pogodan za implementaciju u programskim jezicima koji mogu iskoristiti strukturu transformacijskih matrica (kao što su sparse reprezentacije)

Simplex metoda u svojem izvornom obliku koristi tablični prikaz u kojem se pri svakoj iteraciji ažurira cijela tablica (matrica) koristeći Gauss-Jordanove eliminacije. To je prikladno za manje probleme i obrazovne svrhe, ali za velike probleme može postati neefikasno jer je potrebno višestruko računanje nepotrebnih vrijednosti.

Revidirani Simplex algoritam poboljšava učinkovitost izbjegavanjem potrebe za ažuriranjem cijele tablice. Umjesto toga, radi s manjim skupom podataka – posebno s baznom matricom B i njezinim inverzom. Time se znatno smanjuje memorijska i računska složenost.

Product Form of the Inverse dodatno optimizira Revidirani Simplex tako što se inverz matrice B ne računa ispočetka u svakoj iteraciji, već se postupno ažurira korištenjem elementarnih transformacijskih matrica E koje reflektiraju promjene baze. Svaka transformacija odgovara jednoj zamjeni u bazi (ulazna i izlazna varijabla), a nova inverzna matrica se dobiva kao umnožak prethodne inverzne i matrice E: \(B_{novo}−1=E⋅B_{staro}^{−1}\)

Početni, pripremni koraci, provode se na isti način kao i pri Revidiranom Simplexu.

> library(MASS)
> library(knitr)
> library(kableExtra)
> 
> cvec <- c(3, 5, 0, 0) # koeficijenti uz varijable u funkciji cilja, vektor c prije transponiranja
> 
> A <- matrix(c(2,1,1,0,
+               1,2,0,1), nrow=2, byrow=TRUE) # matrica koeficijenata na lijevoj strani ograničenja
> colnames(A) <- c("x1", "x2", "s1", "s2")
> 
> b <- c(100, 80) # desna strana ograničenja

Iz početnih postavki, izvlačimo dijelove potrebne za daljnje izračune:

> B <- A[, 3:4]  # matrica koeficijenata uz s1 i s2 u ograničenjima (lijeva strana ograničenja)
> B

##      s1 s2
## [1,]  1  0
## [2,]  0  1

> N <- A[, 1:2]  # matrica koeficijenata uz x1 i x2 u ograničenjima (lijeva strana ograničenja)
> N

##      x1 x2
## [1,]  2  1
## [2,]  1  2

> cb <- cvec[3:4]  # koeficijenti uz s1 i s2 u funkciji cilja
> cb

## [1] 0 0

> cN <- cvec[1:2]  # koeficijenti uz x1 i x2 u funkciji cilja
> cN

## [1] 3 5

Inverz matrice \(B\), \(B^{−1}\) je:

> B_inv <- diag(2)  # Inverz identiteta jer B = I
> B_inv

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

Inicijalna bazna matrica

> xB <- B_inv %*% b
> Z <- t(cb) %*% xB
> zj <- t(cb) %*% B_inv %*% N
> cj_zj <- cN - zj
> 
> result1 <- data.frame(
+   "Bazicne varijable" = c("s1", "s2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> result1

##   Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..x1
## 1                s1            100                          3
## 2                s2             80                          3
##   Reduced.costs..cj...zj..x2
## 1                          5
## 2                          5

U prvoj iteraciji \(x_2\) ulazi, \(s_2\) izlazi iz bazičnih varijabli.

> a_j <- A[, 2]  # stupac za x2
> d <- B_inv %*% a_j
> ratio <- ifelse(d > 0, xB / d, Inf)
> leaving_index <- which.min(ratio)
> leaving_var <- c("s1", "s2")[leaving_index]
> 
> E <- diag(2)
> E[, leaving_index] <- -d / d[leaving_index]
> E[leaving_index, leaving_index] <- 1 / d[leaving_index]
> B_inv <- E %*% B_inv
> 
> # Ažuriramo bazu: s1, x2
> cb <- c(0, 5)
> xB <- B_inv %*% b
> Z <- t(cb) %*% xB
> N <- A[, c(1,4)]
> cN <- cvec[c(1,4)]
> zj <- t(cb) %*% B_inv %*% N
> cj_zj <- cN - zj
> 
> result2 <- data.frame(
+   "Bazicne varijable" = c("s1", "x2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> result2

##   Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..x1
## 1                s1             60                        0.5
## 2                x2             40                        0.5
##   Reduced.costs..cj...zj..s2
## 1                       -2.5
## 2                       -2.5

U drugoj iteraciji, \(x_1\) ulazi, \(s_1\) izlazi iz bazičnih varijabli.

> a_j <- A[, 1]  # stupac za x1
> d <- B_inv %*% a_j
> ratio <- ifelse(d > 0, xB / d, Inf)
> leaving_index <- which.min(ratio)
> leaving_var <- c("s1", "x2")[leaving_index]
> 
> E <- diag(2)
> E[, leaving_index] <- -d / d[leaving_index]
> E[leaving_index, leaving_index] <- 1 / d[leaving_index]
> B_inv <- E %*% B_inv
> 
> # Ažuriramo bazu: x1, x2
> cb <- c(3, 5)
> xB <- B_inv %*% b
> Z <- t(cb) %*% xB
> N <- A[, c(3,4)]
> cN <- cvec[c(3,4)]
> zj <- t(cb) %*% B_inv %*% N
> cj_zj <- cN - zj
> 
> result3 <- data.frame(
+   "Bazicne varijable" = c("x1", "x2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> result3

##   Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..s1
## 1                x1             40                      -0.33
## 2                x2             20                      -0.33
##   Reduced.costs..cj...zj..s2
## 1                      -2.33
## 2                      -2.33

Kriterij zaustavljanja je isti kao i za Revidirani Simplex, pa ovdje možemo zaključiti (ponovo) da (40, 20) čine optimalno rješenje.

LU dekompozicija bazne matrice

LU dekompozicija bazne matrice koristi se kao numerički učinkovit način za rješavanje sustava linearnih jednadžbi bez eksplicitnog računanja inverza matrice \(B^{-1}\).

U kontekstu Revidiranog Simplex algoritma, u svakoj iteraciji potrebno je računati:

\(x_B = B^{-1} b\)
\(z_j = c_B^T B^{-1} a_j\) za sve nebazne varijable

Umjesto da se za svaki novi \(B\) eksplicitno računa inverz, možemo primijeniti LU dekompoziciju: \[PB=LU\]

gdje je:

\(L\) donja trokutasta matrica,
\(U\) gornja trokutasta matrica,
\(P\) permutacijska matrica (pivotiranje radi numeričke stabilnosti).

LU dekompozicija razlaže matricu B na donju (L) i gornju (U) trokutastu matricu, čime se efikasnije rješava sustav linearnih jednadžbi bez eksplicitnog računanja inverzne matrice.

Sustav \(B x = b\) tada rješavamo u dva koraka:

\(L y = P b\)
\(U x = y\)

Ova metoda omogućuje veću numeričku stabilnost i bolje performanse u većim problemima. Zbog toga se i koristi se u industrijskim solverima (npr. CPLEX, Gurobi).

> library(MASS)
> library(kableExtra)
> library(Matrix)
> 
> A <- matrix(c(2,1,1,0,
+               1,2,0,1), nrow=2, byrow=TRUE)
> colnames(A) <- c("x1", "x2", "s1", "s2")
> b <- c(100, 80)
> cvec <- c(3, 5, 0, 0) # funkcija cilja

LU dekompozicija početne bazne matrice

> B <- A[, 3:4]  # Bazna matrica: s1 i s2
> Lu <- lu(B)
> LU <- expand(Lu)
> L <- LU$L
> U <- LU$U
> P <- LU$P
> 
> # Rješavanje L y = P b
> y <- forwardsolve(L, P %*% b)
> # Rješavanje U xB = y
> xB <- backsolve(U, y)
> 
> cb <- cvec[3:4]
> Z <- t(cb) %*% xB
> 
> N <- A[, 1:2]
> cN <- cvec[1:2]
> 
> # izračun B^{-1} * N koristeći LU
> Binv_N <- solve(B, N)
> zj <- t(cb) %*% Binv_N
> cj_zj <- cN - zj
> 
> result1 <- data.frame(
+   "Bazicne varijable" = c("s1", "s2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> result1

##   Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..x1
## 1                s1            100                          3
## 2                s2             80                          3
##   Reduced.costs..cj...zj..x2
## 1                          5
## 2                          5

Za nastavak metode, LU dekompozicija bi se koristila za svaku novu baznu matricu \(B\), bez potpunog invertiranja. Umjesto toga, koristi se faktoriranje jednom, a sustavi se rješavaju korak po korak pomoću forwardsolve i backsolve. Ovo može biti numerički stabilnija i računski efikasnija alternativa za revidirani simplex kada se koristi u programskom okruženju.

Slijedi ažuriranje tablice te druga i treća iteracija.

> # Nova bazna matrica: s1 i x2
> B <- A[, c(3, 2)]  # s1 i x2
> B_sparse <- Matrix(B)
> LU <- lu(B_sparse)
> LU_parts <- expand(LU)
> L <- LU_parts$L
> U <- LU_parts$U
> P <- LU_parts$P
> 
> cb <- cvec[c(3,2)]
> 
> # Novi xB
> y <- forwardsolve(L, as.vector(P %*% b))
> xB <- backsolve(U, y)
> Z <- t(cb) %*% xB
> 
> # Nova nebazna matrica: x1 i s2
> N <- A[, c(1, 4)]
> cN <- cvec[c(1, 4)]
> 
> # izračun B^{-1} * N
> Binv_a1 <- backsolve(U, forwardsolve(L, as.vector(P %*% A[,1])))
> Binv_a2 <- backsolve(U, forwardsolve(L, as.vector(P %*% A[,4])))
> Binv_N <- cbind(Binv_a1, Binv_a2)
> 
> zj <- t(cb) %*% Binv_N
> cj_zj <- cN - zj
> 
> result2 <- data.frame(
+   "Bazicne varijable" = c("s1", "x2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> 
> result2

##   Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..Binv_a1
## 1                s1             60                             0.5
## 2                x2             40                             0.5
##   Reduced.costs..cj...zj..Binv_a2
## 1                            -2.5
## 2                            -2.5

> # Nova bazna matrica: x1 i x2
> B <- A[, c(1, 2)]  # x1 i x2
> B_sparse <- Matrix(B)
> LU <- lu(B_sparse)
> LU_parts <- expand(LU)
> L <- LU_parts$L
> U <- LU_parts$U
> P <- LU_parts$P
> 
> cb <- cvec[c(1,2)]
> 
> # Novi xB
> y <- forwardsolve(L, as.vector(P %*% b))
> xB <- backsolve(U, y)
> Z <- t(cb) %*% xB
> 
> # Nova nebazna matrica: s1 i s2
> N <- A[, c(3, 4)]
> cN <- cvec[c(3, 4)]
> 
> Binv_a1 <- backsolve(U, forwardsolve(L, as.vector(P %*% A[,3])))
> Binv_a2 <- backsolve(U, forwardsolve(L, as.vector(P %*% A[,4])))
> Binv_N <- cbind(Binv_a1, Binv_a2)
> 
> zj <- t(cb) %*% Binv_N
> cj_zj <- cN - zj
> 
> result3 <- data.frame(
+   "Bazicne varijable" = c("x1", "x2"),
+   "Vrijednosti xB" = round(xB, 2),
+   "Reduced costs (cj - zj)" = round(cj_zj, 2)
+ )
> 
> result3

##   Bazicne.varijable Vrijednosti.xB Reduced.costs..cj...zj..Binv_a1
## 1                x1             40                           -0.33
## 2                x2             20                           -0.33
##   Reduced.costs..cj...zj..Binv_a2
## 1                           -2.33
## 2                           -2.33

Primjenom LU dekompozicije za računanje rješenja sustava \(Bx = b\) i \(B^{-1}N\), izbjegavamo direktno računanje inverzne matrice. To pridonosi stabilnosti i učinkovitosti u numeričkoj obradi, posebno kod većih problema. LU pristup posebno je koristan kada se koristi u kombinaciji s revidiranim simplexom, ali se može koristiti i samostalno kao pouzdan metod za rješavanje linearnih sustava unutar simplex koraka.

Paketi dostupni u R-u

npr. paketi dostupni u R-u podržavaju sljedeće pristupe:

Paket	Metoda Simplexa	Vrsta ažuriranja	Napomena
`lpSolve`	Klasični Simplex	Tableau (Gauss-Jordan)	Didaktički, jednostavan, pouzdan
`linprog` (R)	Revidirani Simplex	B-inverz pristup	Manje poznat, ali numerički precizan
`ROI.plugin.lpsolve`	Klasični Simplex	Tableau	Isto kao `lpSolve`
`ROI.plugin.highs`	HiGHS	Revidirani dualni Simplex / IP	Moderna, brza, skalabilna

Primjeri ručnog rješavanja i upotrebe paketa

U nastavku će se prikazati postupci rješavanja za nekoliko primjera iz lekcije Uvod u linearno programiranje - strukturiranje problema.

1. primjer

Funkcija cilja

\[max 40x_1+50x_2\]

Ograničenja

\[x_2≤8\]

\[x_1≤24\]

\[x_1+x_2≤28\]

\[x_1≥0\]

\[x_2≥0\]

kanonski zapis:

\[ \text{max } Z = 40x_1 + 50x_2 \] uz ograničenja

\[ \begin{aligned} x_2 + s_1 &= 8 \\ x_1 + s_2 &= 24 \\ x_1 + x_2 + s_3 &= 28 \\ x_1, x_2, s_1, s_2, s_3 &\ge 0 \end{aligned} \]

prikladan pristup za ovakav problem - koristimo klasični Simplex (tableau) ručno ili putem paketa lpSolve.
1. iteracija

BV	x₁	x₂	s₁	s₂	s₃	RHS
s₁	0	1	1	0	0	8
s₂	1	0	0	1	0	24
s₃	1	1	0	0	1	28
Z	-40	-50	0	0	0	0

Ulazna varijabla: \(x_2\) (najnegativniji koeficijent -50)
Omjeri: 8/1=8; 24/0=∞; 28/1=28 \(\rightarrow\) minimalni = 28 \(\rightarrow\) \(s_1\) izlazi
pivot ćelija je (s₁, x₂)=1
normaliziramo red i eliminiramo ostale koeficijente u stupcu te dobivamo:

BV	x₁	x₂	s₁	s₂	s₃	RHS
x₂	0	1	1	0	0	8
s₂	1	0	0	1	0	24
s₃	1	0	-1	0	1	20
Z	-40	0	50	0	0	400

1. iteracija
Ulazna varijabla: \(x_1\) (jedini negativni Z-koeficijent = -40)
Omjeri: s₂: 24/1=24, s₃: 20/1=20, x₂: \(\inf\) \(\rightarrow\) najmanji je 20 \(\rightarrow\) s₃ izlazi
Pivot na ćeliji (\(s_3\),\(x_1\)).
normaliziramo red i eliminiramo ostale koeficijente u stupcu te dobivamo:

BV	x₁	x₂	s₁	s₂	s₃	RHS
x₂	0	1	1	0	0	8
s₂	0	0	1	1	-1	4
x₁	1	0	-1	0	1	20
Z	0	0	10	0	40	1200

Nema više negativnih koeficijenata u Z-retku \(\rightarrow\) optimalno rješenje:

\[(x_1,x_2)=(20,8)\] \[Z=1200\]

Upravo se takav postupak provodi kad koristimo naredbu solveLP iz paketa linprog. Slijedi minimalni kod koji je dovoljan da bismo došli do rješenja.

> library(lpSolve)
> 
> c <- c(40, 50)
> A <- matrix(c(0, 1,
+               1,0,
+               1, 1),
+             nrow = 3, byrow = TRUE)
> b <- c(8, 24, 28)
> dir <- rep("<=", 3)
> 
> rjesenje <- lp(direction = "max", objective.in = c, const.mat = A, const.dir = dir, const.rhs = b)
> rjesenje

## Success: the objective function is 1200

> rjesenje$solution

## [1] 20  8

2. primjer

Funkcija cilja:
\[ \text{Max } 90x_1 + 100x_2 \quad \]

Ograničenja:
\[ \begin{aligned} 0.7x_1 + x_2 &\leq 630 \quad \\ 0.5x_1 + 0.8x_2 &\leq 600 \quad \\ x_1 + 0.6x_2 &\leq 708 \quad \\ 0.1x_1 + 0.25x_2 &\leq 135 \quad \\ x_1, x_2 &\geq 0 \quad \end{aligned} \]

Model u kanonskom obliku:

Maximiziramo \[Z=90x_1+100x_2\] uz \[ \begin{aligned} 0.7x_1 + x_2 + s_1 &= 630 \\ 0.5x_1 + 0.8x_2 + s_2 &= 600 \\ x_1 + 0.6x_2 + s_3 &= 708 \\ 0.1x_1 + 0.25x_2 + s_4 &= 135 \\ x_1, x_2, s_i &\ge0 \end{aligned} \]

prikladan pristup za ovakav problem - Klasični Simplex (rjeđe revidirani putem lpSolve).
Početna tablica:

BV	x₁	x₂	s₁	s₂	s₃	s₄	RHS
s₁	0.70	1.00	1	0	0	0	630
s₂	0.50	0.80	0	1	0	0	600
s₃	1.00	0.60	0	0	1	0	708
s₄	0.10	0.25	0	0	0	1	135
Z	-90	-100	0	0	0	0	0

Ulazna varijabla: x₂ (najnegativniji koeficijent -100)
Omjeri: s₁:630/1=630; s₂:600/0.8=750; s₃:708/0.6≈1180; s₄:135/0.25=540 \(\rightarrow\) najmanji = 540 \(\rightarrow\) s₄ izlazi
pivot na (s₄, x₂)
1. iteracija

BV	x₁	x₂	s₁	s₂	s₃	s₄	RHS
s₁	0.30	0	1	0	0	-4.00	90
s₂	0.18	0	0	1	0	-3.20	168
s₃	0.76	0	0	0	1	-2.40	384
x₂	0.40	1	0	0	0	4.00	540
Z	-50	0	0	0	0	400.00	54000

Privremeno rješenje: Z(0, 540) = 54000
Ulazna varijabla: x₁ (jedini negativni koef. -50)
Omjeri: x₂:540/0.40=1350; s₁:90/0.30=300; s₂:168/0.18≈933; s₃:384/0.76≈505 \(\rightarrow\) najmanji = 300 \(\rightarrow\) s₁ izlazi
1. iteracija

BV	x₁	x₂	s₁	s₂	s₃	s₄	RHS
x₁	1	0	3.33	0	0	-13.33	300
s₂	0	0	-0.60	1	0	-0.80	114
s₃	0	0	-2.53	0	1	7.73	156
x₂	0	1	-1.33	0	0	9.33	420
Z	0	0	166.67	0	0	-266.67	69000

nema negativnih koef. u Z-retku \(\rightarrow\) dosegnuto optimalno rješenje
Rješenje: \(x_1=300\), \(x_2=420\), \(Z=69000\)
1. iteracija
U Z retku je još negativan koeficijent: \(s_4=−266.67\) → ulazi \(s_4\)
Omjer RHS/stupaca: 156/7.73≈20.17, 420/9.33=45 → izlazi \(s_3\) s3
pivot na \((s_3,s_4)=7.73\)

BV	x₁	x₂	s₁	s₂	s₃	s₄	RHS
x₁	1.00	0.00	-1.03	0	1.72	0	568.97
s₂	0.00	0.00	-0.86	1	0.10	0	130.14
s₄	0.00	0.00	-0.33	0	0.13	1	20.17
x₂	0.00	1.00	1.72	0	-1.21	0	231.72
Z	0.00	0.00	79.31	0	34.48	0	74379.31

Više nema negativnih vrijednosti u Z retku i rješenje je dostignuto. Preostaje još iščitati ih i zapisati.

\[x_1 = 568.9 \\ x_2 = 231.72 \\ s_1 = 0 \\ s_2 = 130.14\\ s_3 = 0\\ s_4 = 20.17\\ Z= 74379.31\]

> library(linprog)
> 
> c <- c(90, 100)
> A <- matrix(c(0.7,1,
+               0.5, 0.8,
+               1, 0.6,
+               0.1, 0.25),
+             nrow = 4, byrow = TRUE)
> b <- c(630, 600, 708, 135)
> dir <- rep("<=", 4)
> 
> rjesenje <- solveLP(cvec = c, bvec = b, Amat = A, maximum = TRUE, const.dir = dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Maximum): 74379.3 
## 
## Iterations in phase 1: 0
## Iterations in phase 2: 2
## Solution
##       opt
## 1 568.966
## 2 231.724
## 
## Basic Variables
##          opt
## 1   568.9655
## 2   231.7241
## S 2 130.1379
## S 4  20.1724
## 
## Constraints
##    actual dir bvec     free    dual dual.reg
## 1 630.000  <=  630   0.0000 79.3103 134.4000
## 2 469.862  <=  600 130.1379  0.0000 130.1379
## 3 708.000  <=  708   0.0000 34.4828 156.0000
## 4 114.828  <=  135  20.1724  0.0000  20.1724
## 
## All Variables (including slack variables)
##          opt cvec    min.c    max.c     marg marg.reg
## 1   568.9655   90   70.000 166.6667       NA       NA
## 2   231.7241  100   54.000 128.5714       NA       NA
## S 1   0.0000    0     -Inf  79.3103 -79.3103    134.4
## S 2 130.1379    0 -333.333  92.0000   0.0000       NA
## S 3   0.0000    0     -Inf  34.4828 -34.4828    156.0
## S 4  20.1724    0 -266.667 242.1053   0.0000       NA

Ovdje se može primijetiti da output daje opširniji izvještaj nego što dobivamo ručnim pristupom. Time ćemo se baviti u sljedećoj lekciji i, za sad, nećemo ulaziti u detalje.

3. primjer

Funkcija cilja

\(\text{max } Z = 3200x_1 + 500x_2 + 700x_3 + 1000x_4\)

Ograničenja:

\(x_1 + x_2 + x_3 + x_4 = 100\)

\(x_3 \geq 20\)

\(x_4 \leq 40\)

\(x_3 - x_2 \geq 10\)

\(x_2 \leq 35\)

\(x_1 \leq 45\)

\(x_1, x_2, x_3, x_4 \geq 0\)

kanonski zapis:

\(\text{max } Z = 3200x_1 + 500x_2 + 700x_3 + 1000x_4\)

uz ograničenja:

\(\begin{aligned} x_1 + x_2 + x_3 + x_4 &= 100 \\ x_3 - s_2 &= 20, \quad s_2 \geq 0 \\ x_4 + s_3 &= 40, \quad s_3 \geq 0 \\ x_3 - x_2 - s_4 &= 10, \quad s_4 \geq 0 \\ x_2 + s_5 &= 35, \quad s_5 \geq 0 \\ x_1 + s_6 &= 45, \quad s_6 \geq 0 \\ x_1, x_2, x_3, x_4 &\geq 0 \end{aligned}\)

prikladan pristup za ovakav problem - Dvostupanjski Simplex (zbog jednakosti i \(\geq\) ograničenja), Big‑M ili two‑phase. Ovdje će se prikazati dvostupanjski simplex.
početna tablica:

BV	x₁	x₂	x₃	x₄	s₂	s₃	s₄	s₅	s₆	a₁	a₂	a₃	RHS
a₁	1	1	1	1	0	0	0	0	0	1	0	0	100
a₂	0	0	1	0	-1	0	0	0	0	0	1	0	20
s₃	0	0	0	1	0	1	0	0	0	0	0	0	40
a₃	0	-1	1	0	0	0	-1	0	0	0	0	1	10
s₅	0	1	0	0	0	0	0	1	0	0	0	0	35
s₆	1	0	0	0	0	0	0	0	1	0	0	0	45
W	-1	0	-3	-1	1	0	1	0	0	0	0	0	-130

započinjemo s Fazom I – uklanjanje (minimizacija broja) umjetnih varijabli; cilj je minimizirati \(W = a_1 + a_2 + a_3\). Ako je na kraju \(W^* = 0\), prelazimo na Fazu II.

Faza I – Inicijalne postavke (\(C_j\) za \(w\) prikazane kao 0 za sve varijable koje nisu artifijelne; Artificijelne varijable se rješavaju putem eliminacije retka)

BV	x₁	x₂	x₃	x₄	s₂	s₃	s₄	s₅	s₆	a₁	a₂	a₄	RHS
a₁	1	1	1	1	0	0	0	0	0	1	0	0	100
a₂	0	0	1	0	-1	0	0	0	0	0	1	0	20
s₃	0	0	0	1	0	1	0	0	0	0	0	0	40
a₄	0	-1	1	0	0	0	-1	0	0	0	0	1	10
s₅	0	1	0	0	0	0	0	1	0	0	0	0	35
s₆	1	0	0	0	0	0	0	0	1	0	0	0	45
w	-1	0	-3	-1	1	0	1	0	0	0	0	0	-130

Ulazni stupac: \(x₃\) (najnegativniji u w-retku, -3).
Omjeri (RHS / pozitivni koef. u stupcu \(x_3\)): a₁:100/1=100, a₂:20/1=20, a₄: 10/1=10 ⇒ izlazi a₄ (minimum 10).
Pivot: element \((a_4,x_3)=1\)

BV	x₁	x₂	x₃	x₄	s₂	s₃	s₄	s₅	s₆	a₁	a₂	a₄	RHS
a₁	1	1	0	1	0	0	1	0	0	1	0	-1	90
a₂	0	0	0	0	-1	0	-1	0	0	0	1	-1	10
s₃	0	0	0	1	0	1	0	0	0	0	0	0	40
x₃	0	-1	1	0	0	0	-1	0	0	0	0	1	10
s₅	0	1	0	0	0	0	0	1	0	0	0	0	35
s₆	1	0	0	0	0	0	0	0	1	0	0	0	45
w	-1	-1	0	-1	1	0	0	0	0	0	0	3	-100

Ulazni stupac: x₄ (npr. −1 najnegativniji; i \(x_1\) je -1, ali uzimamo \(x_4\), jer ima pogodniju vrijednost na poziciji za pivotiranje).

Omjeri (pozitivni u \(x_4\) stupcu): a₁: 90/1=90, s₃: 40/1=40 ⇒ izlazi s₃.

Pivot: \((s_3,x_4)=1\)

BV	x₁	x₂	x₃	x₄	s₂	s₃	s₄	s₅	s₆	a₁	a₂	a₄	RHS
a₁	1	1	0	0	0	-1	1	0	0	1	0	-1	50
a₂	0	0	0	0	-1	0	-1	0	0	0	1	-1	10
x₄	0	0	0	1	0	1	0	0	0	0	0	0	40
x₃	0	-1	1	0	0	0	-1	0	0	0	0	1	10
s₅	0	1	0	0	0	0	0	1	0	0	0	0	35
s₆	1	0	0	0	0	0	0	0	1	0	0	0	45
w	-1	-1	0	0	1	1	0	0	0	0	0	3	-60

Ulazni stupac: x₁ (−1)
Omjeri (pozitivni u \(x_1\) stupcu): a₁: 50/1=50, s₆: 45/1=45 ⇒ izlazi s₆.
Pivot: \((s_6,x_1)=1\)

BV	x₁	x₂	x₃	x₄	s₂	s₃	s₄	s₅	s₆	a₁	a₂	a₄	RHS
a₁	0	1	0	0	0	-1	1	0	-1	1	0	-1	5
a₂	0	0	0	0	-1	0	-1	0	0	0	1	-1	10
x₄	0	0	0	1	0	1	0	0	0	0	0	0	40
x₃	0	-1	1	0	0	0	-1	0	0	0	0	1	10
s₅	0	1	0	0	0	0	0	1	0	0	0	0	35
x₁	1	0	0	0	0	0	0	0	1	0	0	0	45
w	0	-1	0	0	1	1	0	0	1	0	0	3	-15

Još je negativan koeficijent u w-retku: \(x₂ = -1\) ⇒ ulazni stupac x₂.

Omjeri (pozitivni u \(x_2\) stupcu): a₁: 5/1=5, s₅: 35/1=35 ⇒ izlazi a₁.
Pivot: \((a_1,x_2)=1\)

Ovime završava Faza I, sve umjetne varijable su maknute.

BV	x₁	x₂	x₃	x₄	s₂	s₃	s₄	s₅	s₆	a₁	a₂	a₄	RHS
x₂	0	1	0	0	0	-1	1	0	-1	1	0	-1	5
a₂	0	0	0	0	-1	0	-1	0	0	0	1	-1	10
x₄	0	0	0	1	0	1	0	0	0	0	0	0	40
x₃	0	-1	1	0	0	0	-1	0	0	0	0	1	10
s₅	0	1	0	0	0	0	0	1	0	0	0	0	35
x₁	1	0	0	0	0	0	0	0	1	0	0	0	45
w	0	0	0	0	0	0	0	0	0	0	0	0	0

U zadnjem koraku Faze I, red \(a_2\) se lako zamijeni s varijablom x₃ (jer se stupac \(x_3\)može promatrati kao jedinični vektor), čime se \(a_2\) izbacuje iz baze. Rezultat: bazično dopušteno rješenje bez umjetnih varijabli.
Uzmemo prethodnu bazu {\(x_1,x_2,x_3,x_4,s_5\)} i formiramo Z-redak. Započinjemo Fazu II.

BV	x₁	x₂	x₃	x₄	s₃	s₄	s₅	s₆	RHS
x₂	0	1	0	0	-1	1	0	-1	5
x₄	0	0	0	1	1	0	0	0	40
x₃	0	-1	1	0	0	-1	0	0	10
s₅	0	1	0	0	0	0	1	0	35
x₁	1	0	0	0	0	0	0	1	45
Z	-3200	-500	-700	-1000	0	0	0	0	0

Redukcijom Z-retka prema trenutačnoj bazi, vidimo da je najisplativije povećavati \(x_1\) (najnegativniji −3200), zatim \(x_4\) (−1000), zatim \(x_3\) (−700), pa \(x_2\) (−500).
Uz ograničenja, isplativo je „napuniti” gornje granice za skuplje varijable i poštovati minimum za \(x_3\) i relaciju \(x_3−x_2≥10\).
također, mora biti zadovoljen minimum za \(x_3\), \(x_3≥20\), poštujući sumu \(x_1+x_2+x_3+x_4=100\) i relaciju \(x_3−x_2≥10\), najviše preostalih jedinica (nakon \(x_1\) i \(x_3\)) isplati se dodijeliti \(x_4\) (jer \(1000>700>500\)), te držati \(x2\) što manjim.
Jedno takvo bazično rješenje koje je zgodno za nastavak je: \((x_1,x_2,x_3,x_4)=(30,10,20,40)\),

Trenutno bazično rješenje (na početku Faze II)

Var	Vrijednost
x₁	30
x₂	10
x₃	20
x₄	40
Z	155000

Aktivna (vezana) ograničenja

x₁ + x₂ + x₃ + x₄ = 100
x₃ = 20 (donja granica)
x₄ = 40 (gornja granica)
x₃ − x₂ = 10

Intuicija za ulaznu varijablu: najveći koeficijent u funkciji cilja je \(x_1\) (3200), zatim \(x_4\) (1000), \(x_3\) (700), \(x_2\) (500). U ovom trenutku je \(x4\) već na maksimumu (40), \(x_3\) je na minimumu (20) i ne smijemo ga smanjiti. Povećavati se “najviše isplati” \(x_1\).

Ako povećamo \(x_1\) za \(t\), zbog jednadžbe zbroja i aktivnih ograničenja \(x_3=20\), \(x_4=40\), moramo smanjiti \(x_2\) za \(t\) (da bi zbroj ostao 100). Ograničenja koja limitiraju \(t\) (Ratio-test ograničenja na porast t za x₁):

Ograničenje	Kako utječe kad x₁ ↑ t	Gornja granica t
x₂ ≥ 0	x₂ → 10 − t	t ≤ 10
x₁ ≤ 45	x₁ → 30 + t	t ≤ 15
x₃ = 20 (donja granica)	x₃ ostaje 20	—
x₄ = 40 (gornja granica)	x₄ ostaje 40	—
x₃ − x₂ ≥ 10	(20) − (10 − t) ≥ 10	t ≥ 0 (zadovoljen)

Najmanja gornja granica je \(t=10\) ⇒ izlazi \(x_2\)

Var	Prije	Promjena	Poslije
x₁	30	+10	40
x₂	10	−10	0
x₃	20	0	20
x₄	40	0	40
Z	155000	+27000	182000

Sada je \(x_2=0\) i ne može dalje padati. Ako želimo još povećati \(x_1\), jednadžba zbroja nalaže da netko drugi mora pasti. \(x_3\) je na donjoj granici (20) i ne smije pasti, pa “prebacujemo” jedinice s \(x_4\) na \(x_1\). Ako \(x_1\) poraste za \(t\), tada \(x_4\) pada za \(t\). Ratio-test (ograničenja na porast t za x₁):

Ograničenje	Kako utječe kad x₁ ↑ t	Gornja granica t
x₁ ≤ 45	x₁ → 40 + t	t ≤ 5
x₄ ≥ 0	x₄ → 40 − t	t ≤ 40
x₃ = 20 (donja granica)	x₃ ostaje 20	—
x₃ − x₂ ≥ 10 (x₂=0)	20 − 0 ≥ 10	✓ (zadovoljen)

Najmanja gornja granica je \(t=5\) (zbog \(x_1≤45\)) ⇒ izlazi \(x_4\). Nakon pivota (\(t=5\)):

Var	Prije	Promjena	Poslije
x₁	40	+5	45
x₂	0	0	0
x₃	20	0	20
x₄	40	−5	35
Z	182000	+11000	193000

Test optimalnosti

Sada su ključne granice aktivne:

\(x_1=45\) (gornja granica) — ne može rasti.
\(x_3=20\) (donja granica) — ne može pasti.
Ako bismo htjeli povećati \(x_4\), morali bismo nekoga smanjiti (zbroj=100). Jedini koji ima manji koeficijent od \(x_4\) je \(x_3\) (700), ali \(x_3\) je na minimumu i ne može se smanjiti. Smanjivati \(x_1\) (3200) ili povećavati \(x_2\) (500) bi smanjilo Z. Dakle, nema dopuštenog smjera poboljšanja ⇒ optimalno rješenje je:

Var	Vrijednost
x₁	45
x₂	0
x₃	20
x₄	35
Z	193000

Može se uočiti da, što su problemi kompleksniji, iako je procedura šablonska, postaje ih teže (dugotrajnije) računati ih ručno. Naravno da nemamo vremena za to u današnjem ubrzanom poslovnom svijetu, pa preferiramo računanje uz računalnu potporu. Mi ćemo za to koristiti R, kao što ste mogli naslutiti u prethodno prikazanim primjerima. Primjeri u nastavku će se rješavati koristeći pakete linprog i lpSolve - dva jednostavna i najčešće korištena paketa. Toplo preporučam pročitati dokumentaciju paketa.

4. primjer

4. primjer a)

Funkcija cilja:

\[ \text{Min } 8x_1 + 3x_2 \]

Ograničenja:

\[ \begin{aligned} 50x_1 + 100x_2 &\leq 1200000 \\ 0.1 \cdot 50 \cdot x_1 + 0.04 \cdot 100 \cdot x_2 &\geq 60000 \\ x_2 &\geq \frac{300000}{100} \\ x_1, x_2 &\geq 0 \end{aligned} \]

Kanonski zapis:

\[ z = 8x_1 + 3x_2\]

Ograničenja	Jednadžbe u kanonskom obliku	Napomena
(1)	\(50x₁ + 100x₂ + s₁ = 1 200 000\)	s₁ ≥ 0 (slack, jer ≤)
(2)	\(5x₁ + 4x₂ − s₂ + a₂ = 60 000\)	s₂ ≥ 0 (surplus), a₂ ≥ 0 (artificial za ≥)
(3)	\(x₂ − s₃ + a₃ = 3 000\)	s₃ ≥ 0 (surplus), a₃ ≥ 0 (artificial za ≥)
	\(x₁, x₂, s₁, s₂, s₃, a₂, a₃ ≥ 0\)

> c <- c(8, 3)
> 
> A <- matrix(c(50, 100,
+               5,   4,
+               0,   1),
+             nrow = 3, byrow = TRUE)
> 
> dir <- c("<=", ">=", ">=")
> 
> b <- c(1200000, 60000, 3000)
> 
> rjesenje <- solveLP(c, b, A, maximum = FALSE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Minimum): 62000 
## 
## Iterations in phase 1: 2
## Iterations in phase 2: 1
## Solution
##     opt
## 1  4000
## 2 10000
## 
## Basic Variables
##       opt
## 1    4000
## 2   10000
## S 3  7000
## 
## Constraints
##    actual dir    bvec free      dual dual.reg
## 1 1200000  <= 1200000    0 0.0566667   420000
## 2   60000  >=   60000    0 2.1666667    42000
## 3   10000  >=    3000 7000 0.0000000     7000
## 
## All Variables (including slack variables)
##       opt cvec       min.c max.c      marg marg.reg
## 1    4000    8 -12.2500000    NA        NA       NA
## 2   10000    3          NA   6.4        NA       NA
## S 1     0    0  -0.0566667   Inf 0.0566667   420000
## S 2     0    0  -2.1666667   Inf 2.1666667    42000
## S 3  7000    0          NA   3.4 0.0000000       NA

4. primjer b)

Funkcija cilja:

\[ \text{Max } 0.1 \cdot 50 \cdot x_1 + 0.04 \cdot 100 \cdot x_2 \]

Ograničenja:

\[ \begin{aligned} 50x_1 + 100x_2 &\leq 1200000 \\ 0.1 \cdot 50 \cdot x_1 + 0.04 \cdot 100 \cdot x_2 &\geq 60000 \\ x_2 &\geq \frac{300000}{100} \\ x_1, x_2 &\geq 0 \end{aligned} \]

> c <- c(0.1*50, 0.04*100)
> 
> A <- matrix(c(50, 100,
+               5,   4,
+               0,   1),
+             nrow = 3, byrow = TRUE)
> 
> dir <- c("<=", ">=", ">=")
> 
> b <- c(1200000, 60000, 3000)
> 
> rjesenje <- solveLP(c, b, A, maximum = TRUE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Maximum): 102000 
## 
## Iterations in phase 1: 2
## Iterations in phase 2: 1
## Solution
##     opt
## 1 18000
## 2  3000
## 
## Basic Variables
##       opt
## 1   18000
## 2    3000
## S 2 42000
## 
## Constraints
##    actual dir    bvec  free dual dual.reg
## 1 1200000  <= 1200000     0  0.1   420000
## 2  102000  >=   60000 42000  0.0    42000
## 3    3000  >=    3000     0  6.0     7000
## 
## All Variables (including slack variables)
##       opt cvec min.c max.c marg marg.reg
## 1   18000    5     2   Inf   NA       NA
## 2    3000    4  -Inf  10.0   NA       NA
## S 1     0    0  -Inf   0.1 -0.1   420000
## S 2 42000    0    -1   Inf  0.0       NA
## S 3     0    0  -Inf   6.0 -6.0     7000

5. primjer

Funkcija cilja:

\[\text{max } 65x + 76y + 75z + 68w\]

Ograničenja:

\[9x + 8y + 10z + 11w \leq 7200\]

\[3x + 3y + 4z \leq 2000\]

\[x + y + z + w \leq 800\]

\[z \geq 100\]

\[x, y, z, w \geq 0\]

> c <- c(65, 76, 75, 68)
> 
> A <- matrix(c(9, 8, 10, 11,
+               3, 3, 4, 0,
+               1, 1, 1, 1,
+               0, 0, 1, 0),
+             nrow = 4, byrow = TRUE)
> 
> dir <- c("<=", "<=", "<=", ">=")
> 
> b <- c(7200, 2000, 800, 100)
> 
> rjesenje <- solveLP(c, b, A, maximum = TRUE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Maximum): 59366.7 
## 
## Iterations in phase 1: 1
## Iterations in phase 2: 2
## Solution
##       opt
## 1   0.000
## 2 533.333
## 3 100.000
## 4 166.667
## 
## Basic Variables
##         opt
## 2   533.333
## 3   100.000
## 4   166.667
## S 1 100.000
## 
## Constraints
##   actual dir bvec free     dual dual.reg
## 1   7100  <= 7200  100  0.00000 100.0000
## 2   2000  <= 2000    0  2.66667 100.0000
## 3    800  <=  800    0 68.00000 166.6667
## 4    100  >=  100    0  3.66667  33.3333
## 
## All Variables (including slack variables)
##         opt cvec    min.c    max.c      marg marg.reg
## 1     0.000   65     -Inf 76.00000 -11.00000 100.0000
## 2   533.333   76 73.25000      Inf        NA       NA
## 3   100.000   75     -Inf 78.66667        NA       NA
## 4   166.667   68  0.00000 76.00000        NA       NA
## S 1 100.000    0 -1.22222  6.18182   0.00000       NA
## S 2   0.000    0     -Inf  2.66667  -2.66667 100.0000
## S 3   0.000    0     -Inf 68.00000 -68.00000 166.6667
## S 4   0.000    0     -Inf  3.66667  -3.66667  33.3333

6. primjer

Funkcija cilja:

\(min\ 0.02x_1 + 0.01x_2 + 0.005x_3 + 0.03x_4 + 0.025x_5\)

Ograničenja:

\(0.14x_1 + 0.25x_2 + 0.037x_3 + 0.02x_4 + 0.15x_5 \geq 25\)

\(4x_1 + 3.3x_2 + 0.55x_3 + 2x_4 + 1.5x_5 \leq 600\)

\(0.05x_3 + 0.10x_4 \geq 15\)

\(x_3 \geq 100\)

\(x_4 \geq 50\)

\(x_1 + x_2 + x_3 + x_4 + x_5 \leq 500\)

\(0.02x_1 + 0.01x_2 + 0.005x_3 + 0.03x_4 + 0.025x_5 \leq 7\)

\(x_2 \geq 50\)

\(x_5 \geq 60\)

\(x_1, x_2, x_3, x_4, x_5 \geq 0\)

Napomena: bit će Vam korisno sastaviti kanonski zapis strukture problema prije pisanja koda (obratite pozornost na nule u matrici lijeve strane ograničenja). Kad izvježbate, ovaj korak ćete vjerojatno preskakati.

> c <- c(0.02, 0.01, 0.005, 0.03, 0.025)
> 
> A <- matrix(c(0.14, 0.25, 0.037, 0.02, 0.15,
+               4, 3.3, 0.55, 2, 1.5,
+               0, 0, 0.05, 0.10, 0,
+               0, 0, 1, 0, 0,
+               0, 0, 0, 1, 0,
+               1, 1, 1, 1, 1,
+               0, 1, 0, 0, 0,
+               0, 0, 0, 0, 1,
+               0.02, 0.01,0.005, 0.03, 0.025),
+             ncol = 5, byrow = TRUE)
> 
> dir <- c(">=", "<=", ">=", ">=", ">=", "<=", "<=", ">=", ">=")
> 
> b <- c(25, 600, 15, 100, 50, 500, 7, 50, 60)
> 
> rjesenje <- solveLP(c, b, A, maximum = FALSE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Minimum): 4.81544 
## 
## Iterations in phase 1: 6
## Iterations in phase 2: 5
## Solution
##         opt
## 1   2.02913
## 2   7.00000
## 3 390.97087
## 4  50.00000
## 5  50.00000
## 
## Basic Variables
##           opt
## 1     2.02913
## 2     7.00000
## 3   390.97087
## 4    50.00000
## 5    50.00000
## S 2 178.74951
## S 3   9.54854
## S 4 290.97087
## S 9  55.18456
## 
## Constraints
##     actual dir bvec      free       dual  dual.reg
## 1  25.0000  >=   25   0.00000 0.14563107   5.33658
## 2 421.2505  <=  600 178.74951 0.00000000 178.74951
## 3  24.5485  >=   15   9.54854 0.00000000   9.54854
## 4 390.9709  >=  100 290.97087 0.00000000 290.97087
## 5  50.0000  >=   50   0.00000 0.02747573  88.51538
## 6 500.0000  <=  500   0.00000 0.00038835 140.50000
## 7   7.0000  <=    7   0.00000 0.02601942   7.00000
## 8  50.0000  >=   50   0.00000 0.00354369   1.84956
## 9 115.1846  >=   60  55.18456 0.00000000  55.18456
## 
## All Variables (including slack variables)
##           opt  cvec       min.c       max.c       marg  marg.reg
## 1     2.02913 0.020 -0.02108108 0.023230088         NA        NA
## 2     7.00000 0.010          NA 0.036019417         NA        NA
## 3   390.97087 0.005 -0.04150000 0.005285714         NA        NA
## 4    50.00000 0.030 -0.05747573         Inf         NA        NA
## 5    50.00000 0.025 -0.02854369         Inf         NA        NA
## S 1   0.00000 0.000 -0.14563107         Inf 0.14563107   5.33658
## S 2 178.74951 0.000 -0.00125000 0.000563380 0.00000000        NA
## S 3   9.54854 0.000 -0.65813953 0.005714286 0.00000000        NA
## S 4 290.97087 0.000 -0.03650000 0.000285714 0.00000000        NA
## S 5   0.00000 0.000 -0.02747573         Inf 0.02747573  88.51538
## S 6   0.00000 0.000 -0.00038835         Inf 0.00038835 140.50000
## S 7   0.00000 0.000 -0.02601942         Inf 0.02601942   7.00000
## S 8   0.00000 0.000 -0.00354369         Inf 0.00354369   1.84956
## S 9  55.18456 0.000          NA 1.000000000 0.00000000        NA

7. primjer

Funkcija cilja:

\[min\ x_1+x_2+x_3+x_4+x_5\]

Ograničenja:

\[x_1+x_2+x_3+x_4+x_5≤2500\]

\[x_1≥600\]

\[x_2\geq 200\]

\[x_3≤400\]

\[x_4\geq100\]

\[x_4≤500\]

\[x_5≥150\]

\[x_1+x_2≤1000\]

\[x_3+x_4<800\]

\[x_1+x_2+x_5≤1500\]

\[x_1+x_2+x_3+x_4+x_5≤2300\]

\[x_1, x_2, x_3, x_4, x_5 \geq 0\]

> c <- c(1, 1, 1, 1, 1)
> 
> A <- matrix(
+   c(
+     1, 1, 1, 1, 1,   # (1)
+     1, 0, 0, 0, 0,   # (2)
+     0, 1, 0, 0, 0,   # (3)
+     0, 0, 1, 0, 0,   # (4)
+     0, 0, 0, 1, 0,   # (5)
+     0, 0, 0, 1, 0,   # (6)
+     0, 0, 0, 0, 1,   # (7)
+     1, 1, 0, 0, 0,   # (8)
+     0, 0, 1, 1, 0,   # (9)
+     1, 1, 0, 0, 1,   # (10)
+     1, 1, 1, 1, 1    # (11)
+   ),
+   nrow = 11, byrow = TRUE
+ )
> 
> dir <- c(
+   "<=",  # (1)
+   ">=",  # (2)
+   ">=",  # (3)
+   "<=",  # (4)
+   ">=",  # (5)
+   "<=",  # (6)
+   ">=",  # (7)
+   "<=",  # (8)
+   "<=",  # (9)  (iz < pretvoreno u <= s malim razmakom); alternativno koristimo lpSolve
+   "<=",  # (10)
+   "<="   # (11)
+ )
> 
> b <- c(
+   2500,          # (1)
+   600,           # (2)
+   200,           # (3)
+   400,           # (4)
+   100,           # (5)
+   500,           # (6)
+   150,           # (7)
+   1000,          # (8)
+   800 - 1e-9,    # (9)   strogo < 800
+   1500,          # (10)
+   2300           # (11)
+ )
> 
> rjesenje <- solveLP(c, b, A, maximum = FALSE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Minimum): 1050 
## 
## Iterations in phase 1: 4
## Iterations in phase 2: 0
## Solution
##   opt
## 1 600
## 2 200
## 3   0
## 4 100
## 5 150
## 
## Basic Variables
##       opt
## 1     600
## 2     200
## 4     100
## 5     150
## S 1  1450
## S 4   400
## S 6   400
## S 8   200
## S 9   700
## S 10  550
## S 11 1250
## 
## Constraints
##    actual dir bvec free dual dual.reg
## 1    1050  <= 2500 1450    0     1450
## 2     600  >=  600    0    1      200
## 3     200  >=  200    0    1      200
## 4       0  <=  400  400    0      400
## 5     100  >=  100    0    1      400
## 6     100  <=  500  400    0      400
## 7     150  >=  150    0    1      550
## 8     800  <= 1000  200    0      200
## 9     100  <=  800  700    0      700
## 10    950  <= 1500  550    0      550
## 11   1050  <= 2300 1250    0     1250
## 
## All Variables (including slack variables)
##       opt cvec min.c max.c marg marg.reg
## 1     600    1    -2   Inf   NA       NA
## 2     200    1    -2   Inf   NA       NA
## 3       0    1    99    77    1      400
## 4     100    1    -2   Inf   NA       NA
## 5     150    1    -2   Inf   NA       NA
## S 1  1450    0    NA     1    0       NA
## S 2     0    0    -1   Inf    1      200
## S 3     0    0    -1   Inf    1      200
## S 4   400    0    NA     1    0       NA
## S 5     0    0    -1   Inf    1      400
## S 6   400    0    NA     1    0       NA
## S 7     0    0    -1   Inf    1      550
## S 8   200    0    NA     1    0       NA
## S 9   700    0    NA     1    0       NA
## S 10  550    0    NA     1    0       NA
## S 11 1250    0    NA     1    0       NA

8. primjer

Funkcija cilja:

\[max (0.05CF + 0.08VC + 0.03MK + 0.07PA + 0.04OB)\]

Ograničenja:

\[CF + VC + MK + PA + OB = 1000000\]

\[0.03VC + 0.015PA ≤ 300000\]

\[CF \leq 200000\]

\[VC \leq 500000\]

\[MK \leq 150000\]

\[PA \leq 250000\]

\[OB \leq 400000\]

\[CF, VC, MK, PA, OB \geq 0\]

> c <- c(0.05, 0.08, 0.03, 0.07, 0.04)
> 
> A <- matrix(
+   c(
+     1,     1,     1,     1,     1,      # (1)
+     0,   0.03,    0,   0.015,   0,      # (2)
+     1,     0,     0,     0,     0,      # (3)
+     0,     1,     0,     0,     0,      # (4)
+     0,     0,     1,     0,     0,      # (5)
+     0,     0,     0,     1,     0,      # (6)
+     0,     0,     0,     0,     1       # (7)
+   ),
+   nrow = 7, byrow = TRUE
+ )
> 
> dir <- c(
+   ">=",   # (1) kako bismo dobili rješenje koristeći solveLP(), u ovom slučaju smijemo "=" zamijeniti s ">=", jer će minimizacija i ostala ograničenja prisiliti vrijednost na donju granicu, tj. na "="; ovakve situacije uvijek moramo specifično naznačiti i objasniti
+   "<=",  # (2)
+   "<=",  # (3)
+   "<=",  # (4)
+   "<=",  # (5)
+   "<=",  # (6)
+   "<="   # (7)
+   )
> 
> b <- c(
+   1000000,  # (1)
+   300000,   # (2)
+   200000,   # (3)
+   500000,   # (4)
+   150000,   # (5)
+   250000,   # (6)
+   400000    # (7)
+   )
> 
> rjesenje <- solveLP(c, b, A, maximum = FALSE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Minimum): 48000 
## 
## Iterations in phase 1: 3
## Iterations in phase 2: 3
## Solution
##      opt
## 1 200000
## 2      0
## 3 150000
## 4 250000
## 5 400000
## 
## Basic Variables
##        opt
## 1   200000
## 3   150000
## 4   250000
## 5   400000
## S 2 296250
## S 4 500000
## S 6      0
## 
## Constraints
##    actual dir    bvec   free dual dual.reg
## 1 1000000  >= 1000000      0 0.07 19750000
## 2    3750  <=  300000 296250 0.00   296250
## 3  200000  <=  200000      0 0.02   200000
## 4       0  <=  500000 500000 0.00   500000
## 5  150000  <=  150000      0 0.04   150000
## 6  250000  <=  250000      0 0.00       NA
## 7  400000  <=  400000      0 0.03   400000
## 
## All Variables (including slack variables)
##        opt cvec min.c     max.c marg marg.reg
## 1   200000 0.05    NA  0.070000   NA       NA
## 2        0 0.08 99.00 77.000000 0.01   250000
## 3   150000 0.03    NA  0.070000   NA       NA
## 4   250000 0.07 -0.09  0.080000   NA       NA
## 5   400000 0.04    NA  0.070000   NA       NA
## S 1      0 0.00 -0.07       Inf 0.07 19750000
## S 2 296250 0.00    NA  0.666667 0.00       NA
## S 3      0 0.00 -0.02       Inf 0.02   200000
## S 4 500000 0.00    NA  0.010000 0.00       NA
## S 5      0 0.00 -0.04       Inf 0.04   150000
## S 6      0 0.00 -0.01  0.020000 0.00       NA
## S 7      0 0.00 -0.03       Inf 0.03   400000

9.primjer

Funkcija cilja:

\[ \text{Max } 65x_1 + 90x_2 + 40x_3 + 60x_4 + 20x_5 \]

Ograničenja:

\[ \begin{aligned} x_1 &\leq 15 \quad \\ x_2 &\leq 10 \quad \\ x_3 &\leq 25 \quad \\ x_4 &\leq 4 \quad \\ x_5 &\leq 30 \quad \end{aligned} \]

\[ 1500x_1 + 3000x_2 + 400x_3 + 1000x_4 + 100x_5 \leq 30000 \]

\[ x_1 + x_2 \geq 10 \]

\[ 1500x_1 + 3000x_2 \leq 18000 \]

\[ 1000x_1 + 2000x_2 + 1500x_3 + 2500x_4 + 300x_5 \geq 50000 \]

\[ x_1, x_2, x_3, x_4, x_5 \geq 0 \]

> c <- c(65, 90, 40, 60, 20)
> 
> A <- matrix(
+   c(
+     1,    0,    0,    0,    0,       # (1)
+     0,    1,    0,    0,    0,       # (2)
+     0,    0,    1,    0,    0,       # (3)
+     0,    0,    0,    1,    0,       # (4)
+     0,    0,    0,    0,    1,       # (5)
+     1500, 3000, 400,  1000, 100,     # (6)
+     1,    1,    0,    0,    0,       # (7)
+     1500, 3000, 0,    0,    0,       # (8)
+     1000, 2000, 1500, 2500, 300      # (9)
+     ),   nrow = 9, byrow = TRUE
+   )
> 
> dir <- c(
+   "<=",  # (1)
+   "<=",  # (2)
+   "<=",  # (3)
+   "<=",  # (4)
+   "<=",  # (5)
+   "<=",  # (6)
+   ">=",  # (7)
+   "<=",  # (8)
+   ">="   # (9)
+   )
> 
> b <- c(
+   15,     # (1)
+   10,     # (2)
+   25,     # (3)
+   4,      # (4)
+   30,     # (5)
+   30000,  # (6)
+   10,     # (7)
+   18000,  # (8)
+   50000   # (9)
+   )
> 
> 
> rjesenje <- solveLP(c, b, A, maximum = TRUE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Maximum): 2370 
## 
## Iterations in phase 1: 4
## Iterations in phase 2: 4
## Solution
##   opt
## 1  10
## 2   0
## 3  25
## 4   2
## 5  30
## 
## Basic Variables
##       opt
## 1      10
## 3      25
## 4       2
## 5      30
## S 1     5
## S 2    10
## S 4     2
## S 8  3000
## S 9 11500
## 
## Constraints
##   actual dir  bvec  free  dual    dual.reg
## 1     10  <=    15     5  0.00     5.00000
## 2      0  <=    10    10  0.00    10.00000
## 3     25  <=    25     0 16.00     5.00000
## 4      2  <=     4     2  0.00     2.00000
## 5     30  <=    30     0 14.00    20.00000
## 6  30000  <= 30000     0  0.06  2000.00000
## 7     10  >=    10     0 25.00     1.33333
## 8  15000  <= 18000  3000  0.00  3000.00000
## 9  61500  >= 50000 11500  0.00 11500.00000
## 
## All Variables (including slack variables)
##       opt cvec        min.c    max.c   marg   marg.reg
## 1      10   65   0.00000000  90.0000     NA         NA
## 2       0   90         -Inf 155.0000 -65.00    1.33333
## 3      25   40  24.00000000      Inf     NA         NA
## 4       2   60  43.33333333 100.0000     NA         NA
## 5      30   20   6.00000000      Inf     NA         NA
## S 1     5    0 -25.00000000  65.0000   0.00         NA
## S 2    10    0 -65.00000000      Inf   0.00         NA
## S 3     0    0         -Inf  16.0000 -16.00    5.00000
## S 4     2    0 -40.00000000  16.6667   0.00         NA
## S 5     0    0         -Inf  14.0000 -14.00   20.00000
## S 6     0    0         -Inf   0.0600  -0.06 2000.00000
## S 7     0    0         -Inf  25.0000 -25.00    1.33333
## S 8  3000    0  -0.01666667      Inf   0.00         NA
## S 9 11500    0  -0.00909091      Inf   0.00         NA

10. primjer

Funkcija cilja:

\[\text{max}\ x_1 + x_2\]

Ograničenja:

\[f_1 + f_2 \leq 2\]

\[g_1 + g_2 \leq 0.6\]

\[k_1 + k_2 \leq 0.15\]

\[b_1 + b_2 \leq 0.15\] \[0.5x_1 - f_1 \leq 0\]

\[f_1 - 0.7x_1 \leq 0\]

\[0.3x_2 - f_2 \leq 0\]

\[f_2 - 0.5x_2 \leq 0\]

\[0.2x_1 - g_1 \leq 0\]

\[g_1 - 0.4x_1 \leq 0\]

\[0.4x_2 - g_2 \leq 0\]

\[g_2 - 0.6x_2 \leq 0\]

\[ 0.05x_1 - k_1 = 0\]

\[ 0.05x_2 - k_2 = 0\]

\[0.05x_1 - b_1 \leq 0\]

\[0.1x_2 - b_2 \leq 0\]

\[b_2 - 0.15x_2 \leq 0\]

\[x_1=f_1+g_1+k_1+b_1\]

\[x_2=f_2+g_2+k_2+b_2\]

\[x_1 \geq 0.5\]

\[x_2 \geq 0.5\]

\[x_1, x_2, f_1, f_2, g_1, g_2, k_1, k_2, b_1, b_2 \geq 0\]

> c <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0)
> 
> A <- matrix(
+   c(
+     0,   0,   1,   1,   0,   0,   0,   0,   0,   0,    # (1)
+     0,   0,   0,   0,   1,   1,   0,   0,   0,   0,    # (2)
+     0,   0,   0,   0,   0,   0,   1,   1,   0,   0,    # (3)
+     0,   0,   0,   0,   0,   0,   0,   0,   1,   1,    # (4)
+     0.5, 0,  -1,   0,   0,   0,   0,   0,   0,   0,    # (5)
+    -0.7, 0,   1,   0,   0,   0,   0,   0,   0,   0,    # (6)
+     0,  0.3,  0,  -1,   0,   0,   0,   0,   0,   0,    # (7)
+     0, -0.5,  0,   1,   0,   0,   0,   0,   0,   0,    # (8)
+     0.2, 0,   0,   0,  -1,   0,   0,   0,   0,   0,    # (9)
+    -0.4, 0,   0,   0,   1,   0,   0,   0,   0,   0,    # (10)
+     0,  0.4,  0,   0,   0,  -1,   0,   0,   0,   0,    # (11)
+     0, -0.6,  0,   0,   0,   1,   0,   0,   0,   0,    # (12)
+     0.05,0,   0,   0,   0,   0,  -1,   0,   0,   0,    # (13)
+     0, 0.05,  0,   0,   0,   0,   0,  -1,   0,   0,    # (14)
+     0.05,0,   0,   0,   0,   0,   0,   0,  -1,   0,    # (15)
+     0,  0.1,  0,   0,   0,   0,   0,   0,   0,  -1,    # (16)
+     0, -0.15, 0,   0,   0,   0,   0,   0,   0,   1,    # (17)
+     1,   0,  -1,   0,  -1,   0,  -1,   0,  -1,   0,    # (18)
+     0,   1,   0,  -1,   0,  -1,   0,  -1,   0,  -1,    # (19)
+     1,   0,   0,   0,   0,   0,   0,   0,   0,   0,    # (20)
+     0,   1,   0,   0,   0,   0,   0,   0,   0,   0     # (21)
+    ),  nrow = 21, byrow = TRUE
+   )
> 
> dir <- c(
+   "<=",  # (1)
+   "<=",  # (2)
+   "<=",  # (3)
+   "<=",  # (4)
+   "<=",  # (5)
+   "<=",  # (6)
+   "<=",  # (7)
+   "<=",  # (8)
+   "<=",  # (9)
+   "<=",  # (10)
+   "<=",  # (11)
+   "<=",  # (12)
+   "=",   # (13)
+   "=",   # (14)
+   "<=",  # (15)
+   "<=",  # (16)
+   "<=",  # (17)
+   "=",   # (18)
+   "=",   # (19)
+   ">=",  # (20)
+   ">="   # (21)
+   )
> 
> b <- c(
+   2,       # (1)
+    0.6,     # (2)
+    0.15,    # (3)
+    0.15,    # (4)
+    0,       # (5)
+    0,       # (6)
+    0,       # (7)
+    0,       # (8)
+    0,       # (9)
+    0,       # (10)
+    0,       # (11)
+    0,       # (12)
+    0,       # (13)
+    0,       # (14)
+    0,       # (15)
+    0,       # (16)
+    0,       # (17)
+    0,       # (18)
+    0,       # (19)
+    0.5,     # (20)
+    0.5      # (21)
+   )
> 
> 
> rjesenje <- solveLP(c, b, A, maximum = TRUE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Maximum): 2.5 
## 
## Iterations in phase 1: 8
## Iterations in phase 2: 1
## Solution
##     opt
## 1  2.00
## 2  0.50
## 3  1.00
## 4  0.15
## 5  0.40
## 6  0.20
## 7  0.00
## 8  0.00
## 9  0.10
## 10 0.05
## 
## Basic Variables
##        opt
## 1    2.000
## 2    0.500
## 3    1.000
## 4    0.150
## 5    0.400
## 6    0.200
## 9    0.100
## 10   0.050
## S 1  0.850
## S 3  0.150
## S 4  0.000
## S 6  0.400
## S 8  0.100
## S 10 0.400
## S 12 0.100
## S 13 0.000
## S 14 0.000
## S 17 0.025
## S 18 0.000
## S 19 0.000
## S 20 1.500
## 
## Constraints
##    actual dir bvec  free dual dual.reg
## 1   1.150  <= 2.00 0.850    0    0.850
## 2   0.600  <= 0.60 0.000    5    0.300
## 3   0.000  <= 0.15 0.150    0    0.150
## 4   0.150  <= 0.15 0.000    0       NA
## 5   0.000  <= 0.00 0.000    0    0.400
## 6  -0.400  <= 0.00 0.400    0    0.400
## 7   0.000  <= 0.00 0.000    0    0.100
## 8  -0.100  <= 0.00 0.100    0    0.100
## 9   0.000  <= 0.00 0.000    5    0.200
## 10 -0.400  <= 0.00 0.400    0    0.400
## 11  0.000  <= 0.00 0.000    5    0.100
## 12 -0.100  <= 0.00 0.100    0    0.100
## 13  0.000   = 0.00 0.000    0       NA
## 14  0.000   = 0.00 0.000    0       NA
## 15  0.000  <= 0.00 0.000    0      Inf
## 16  0.000  <= 0.00 0.000    0    0.025
## 17 -0.025  <= 0.00 0.025    0    0.025
## 18  0.000   = 0.00 0.000    0       NA
## 19  0.000   = 0.00 0.000    0       NA
## 20  2.000  >= 0.50 1.500    0    1.500
## 21  0.500  >= 0.50 0.000    1    0.750
## 
## All Variables (including slack variables)
##        opt cvec min.c    max.c marg marg.reg
## 1    2.000    1   0.5      Inf   NA       NA
## 2    0.500    1  -Inf  2.00000   NA       NA
## 3    1.000    0  -1.0      Inf   NA       NA
## 4    0.150    0  -Inf  3.33333   NA       NA
## 5    0.400    0  -2.5      Inf   NA       NA
## 6    0.200    0  -Inf  2.50000   NA       NA
## 7    0.000    0  -Inf  0.00000    0    0.150
## 8    0.000    0  -Inf  0.00000    0    0.150
## 9    0.100    0 -10.0      Inf   NA       NA
## 10   0.050    0  -Inf 10.00000   NA       NA
## S 1  0.850    0    NA  1.42857    0       NA
## S 2  0.000    0  -Inf  5.00000   -5    0.300
## S 3  0.150    0  -Inf      Inf    0       NA
## S 4  0.000    0    NA 20.00000    0       NA
## S 5  0.000    0  -Inf  0.00000    0    0.400
## S 6  0.400    0  -2.5      Inf    0       NA
## S 7  0.000    0  -Inf  0.00000    0    0.100
## S 8  0.100    0  -Inf  5.00000    0       NA
## S 9  0.000    0  -Inf  5.00000   -5    0.200
## S 10 0.400    0  -2.5      Inf    0       NA
## S 11 0.000    0  -Inf  5.00000   -5    0.100
## S 12 0.100    0  -5.0  5.00000    0       NA
## S 13 0.000    0  -Inf      Inf    0       NA
## S 14 0.000    0  -Inf      Inf    0       NA
## S 15 0.000    0  -Inf  0.00000    0      Inf
## S 16 0.000    0  -Inf  0.00000    0    0.025
## S 17 0.025    0  -Inf 20.00000    0       NA
## S 18 0.000    0  -Inf      Inf    0       NA
## S 19 0.000    0  -Inf      Inf    0       NA
## S 20 1.500    0  -0.5      Inf    0       NA
## S 21 0.000    0  -Inf  1.00000   -1    0.750

11. primjer

Funkcija cilja:

\(min\ ∑C_i \cdot x_i\)

Ograničenja:

\[x_1\geq100\] \[x_2\geq2\] \[x_3\geq50\] \[x_4\geq50\] \[x_5\geq20\] \[x_6\geq2\] \[x_7\geq50\] \[∑C_i×x_i\leq5,000,000\] \[x_2+0.5x_6\geq4\] \[x_1\leq500\] \[x_4\leq400\]

> c <- c(2000, 150000, 100000, 5000, 120000, 20000, 110000)
> 
> A <- matrix(
+   c(
+     1, 0, 0, 0, 0, 0, 0,                 # (1)
+     0, 1, 0, 0, 0, 0, 0,                 # (2)
+     0, 0, 1, 0, 0, 0, 0,                 # (3)
+     0, 0, 0, 1, 0, 0, 0,                 # (4)
+     0, 0, 0, 0, 1, 0, 0,                 # (5)
+     0, 0, 0, 0, 0, 1, 0,                 # (6)
+     0, 0, 0, 0, 0, 0, 1,                 # (7)
+     2000, 150000, 100000, 5000, 120000, 20000, 110000,   # (8)
+     0, 1, 0, 0, 0, 0.5, 0,               # (9)
+     1, 0, 0, 0, 0, 0, 0,                 # (10)
+     0, 0, 0, 1, 0, 0, 0                  # (11)
+     ),  nrow = 11, byrow = TRUE
+   )
> 
> dir <- c(
+   ">=",  # (1)
+   ">=",  # (2)
+   ">=",  # (3)
+   ">=",  # (4)
+   ">=",  # (5)
+   ">=",  # (6)
+   ">=",  # (7)
+   "<=",  # (8)
+   ">=",  # (9)
+   "<=",  # (10)
+   "<="   # (11)
+   )
> 
> b <- c(
+   100,       # (1)
+   2,         # (2)
+   50,        # (3)
+   50,        # (4)
+   20,        # (5)
+   2,         # (6)
+   50,        # (7)
+   5000000,   # (8)
+   4,         # (9)
+   500,       # (10)
+   400        # (11)
+   )
> 
> 
> rjesenje <- solveLP(c, b, A, maximum = FALSE, dir)
> rjesenje

## 
## 
## Results of Linear Programming / Linear Optimization
## 
## Objective function (Minimum): 610000 
## 
## Iterations in phase 1: 5
## Iterations in phase 2: 1
## Solution
##   opt
## 1 100
## 2   0
## 3   0
## 4  50
## 5   0
## 6   8
## 7   0
## 
## Basic Variables
##          opt
## 1        100
## 4         50
## 6          8
## S 2        2
## S 3       50
## S 5       20
## S 6        6
## S 7       50
## S 8  4390000
## S 10     400
## S 11     350
## 
## Constraints
##    actual dir  bvec    free  dual   dual.reg
## 1     100  >= 1e+02       0  2000     400.00
## 2       4  >= 2e+00       2     0       2.00
## 3     100  >= 5e+01      50     0      50.00
## 4      50  >= 5e+01       0  5000     350.00
## 5      40  >= 2e+01      20     0      20.00
## 6       8  >= 2e+00       6     0       6.00
## 7     100  >= 5e+01      50     0      50.00
## 8  610000  <= 5e+06 4390000     0 4390000.00
## 9       4  >= 4e+00       0 40000     109.75
## 10    100  <= 5e+02     400     0     400.00
## 11     50  <= 4e+02     350     0     350.00
## 
## All Variables (including slack variables)
##          opt   cvec  min.c  max.c   marg marg.reg
## 1        100   2000  -4000    Inf     NA       NA
## 2          0 150000     99     77 110000   2.0000
## 3          0 100000     99     77 100000  43.9000
## 4         50   5000 -10000    Inf     NA       NA
## 5          0 120000     99     77 120000  20.0000
## 6          8  20000 -40000  75000     NA       NA
## 7          0 110000     99     77 110000  39.9091
## S 1        0      0  -2000    Inf   2000 400.0000
## S 2        2      0     NA 110000      0       NA
## S 3       50      0     NA 100000      0       NA
## S 4        0      0  -5000    Inf   5000 350.0000
## S 5       20      0     NA 120000      0       NA
## S 6        6      0 -20000  55000      0       NA
## S 7       50      0     NA 110000      0       NA
## S 8  4390000      0     NA      1      0       NA
## S 9        0      0 -40000    Inf  40000 109.7500
## S 10     400      0     NA   2000      0       NA
## S 11     350      0     NA   5000      0       NA

Pitanja za ponavljanje

1) Opći cilj Simplex metode pri maksimizaciji je…
A) pronaći proizvoljno izvedivo rješenje
B) minimizirati broj ograničenja
C) kretati se po rubovima izvedivog poligona do vrha s najvećim \(Z\)
D) izbjegavati jednadžbe u kanonskom zapisu

2) U kanonskom zapisu, ograničenje tipa \(a^\top x \le b\) pretvara se u jednadžbu dodavanjem:
A) surplus varijable
B) slack varijable
C) artificijelne varijable
D) binarne varijable

3) U kanonskom zapisu, ograničenje \(a^\top x \ge b\) pretvara se u jednadžbu dodavanjem:
A) slack varijable i često artificijelne varijable
B) samo slack varijable
C) samo artificijelne varijable
D) ničega — već je jednadžba

4) U početnoj Simplex tablici bazične varijable su one čiji stupci…
A) imaju sve nule
B) imaju 1 u jednom retku i 0 u ostalima
C) imaju najveći koeficijent u cilju
D) imaju negativan koeficijent u Z-retku

5) Kod maksimizacije, varijabla koja ulazi u bazu po „Most Negative Rule” je ona s…
A) najmanjim pozitivnim koeficijentom u Z-retku
B) najnegativnijim koeficijentom u Z-retku
C) najvećim pivot elementom u stupcu
D) najmanjim omjerom RHS/\(a_{ij}\)

6) Pravilo za izbor izlazne varijable (pivot red) u standardnom Simplexu je…
A) najveći omjer \(b_i/a_{ij}\) za \(a_{ij}>0\)
B) najmanji omjer \(b_i/a_{ij}\) za \(a_{ij}>0\)
C) najmanji koeficijent u Z-retku
D) prvi redak po abecedi

7) Ako su svi koeficijenti u stupcu ulazne varijable \(\le 0\), tada je problem…
A) degeneriran
B) nedopustiv
C) neograničen (unbounded)
D) višestruko optimalan

8) Blandovo pravilo služi primarno za…
A) bržu konvergenciju
B) sprječavanje ciklusa/degeneracijskih petlji
C) maksimizaciju rasta \(Z\) po iteraciji
D) odabir najvećeg pivota

9) Reducirani trošak \(\bar{c}_j = c_j - c_B^\top B^{-1}A_j\) u maksimizaciji — uvjet optimalnosti je ispunjen ako za sve nebazične \(j\):
A) \(\bar{c}_j \ge 0\)
B) \(\bar{c}_j < 0\)
C) \(\bar{c}_j = 0\)
D) \(\bar{c}_j \le 0\)

10) Koji pristup je prikladan kad početno bazično rješenje nije očito zbog „=” i „\(\ge\)” ograničenja?
A) Two-phase (dvostupanjski) Simplex ili Big-M
B) Revidirani Simplex C) Nasumični pivot
D) Grafička metoda

11) Revidirani Simplex (u odnosu na tablični/„tableau”)…
A) uvijek ažurira cijelu tablicu
B) radi s baznom matricom \(B\) i \(B^{-1}\), bez pune tablice
C) ne koristi reducirane troškove
D) nije pogodan za velike probleme

12) Dualni Simplex je najprikladniji kada…
A) je Z-redak „uređen” (dualno dopustivo), ali neki RHS su negativni (primarno nedopustivo)
B) imamo stroge „<” nejednakosti
C) želimo nasumični izbor pivota
D) tražimo grafičko rješenje

13) Steepest-edge pravilo za izbor ulazne varijable maksimizira približno…
A) \(|c_j|\)
B) \(|c_j|/\|a_j\|\)
C) \(\sum \min b_i/a_{ij}\)
D) broj iteracija

14) Devex metoda je…
A) točna implementacija steepest-edge
B) heuristika koja aproksimira steepest-edge s težinama
C) isto što i Blandovo pravilo
D) tehnika LU dekompozicije

15) Kod stroge nejednakosti „\(<\)” u LP-u tipično…
A) pišemo je kao „\(\le\)” bez ikakve napomene
B) zamjenjujemo ju „\(\ge\)”
C) modeliramo kao \(\le\) s vrlo malom razlikom (npr. \(b - 10^{-9}\))
D) nije moguće riješiti LP

16) U kojoj situaciji je preporučljivo koristiti LU dekompoziciju ili PFI u Simplexu?
A) U malim problemima s 2 varijable
B) U velikim problemima radi stabilnosti i efikasnosti umjesto eksplicitnog \(B^{-1}\)
C) Samo kad je \(B\) singularna
D) Kad postoji jedino rješenje

17) U početnoj tablici:

BV	x₁	x₂	s₁	s₂	RHS
s₁	2	1	1	0	100
s₂	1	2	0	1	80
Z	-3	-5	0	0	0

Po „Most Negative Rule” ulazi:
A) \(x_1\)
B) \(x_2\)
C) \(s_1\)
D) \(s_2\)

18) Za tablicu iz 17) i odabranu ulaznu varijablu, tko izlazi po pravilu minimalnog omjera?
(Računajte \(\text{RHS}/a_{ij}\) za pozitivne \(a_{ij}\) u tom stupcu.)
A) \(s_1\)
B) \(s_2\)
C) niti jedna — neograničeno
D) veza — treba Bland

19) U Fazi I dvostupanjskog Simplexa cilj je…
A) maksimizirati \(Z\)
B) pronaći LU dekompoziciju \(B\)
C) minimizirati broj slack varijabli
D) minimizirati zbroj artificijelnih varijabli \(W\) i ukloniti ih iz baze

20) Koja je izjava o optimalnosti u Simplexu za maksimizaciju najtočnija?
A) Ako postoji barem jedan negativan koeficijent u Z-retku, rješenje je optimalno
B) Ako su svi reducirani troškovi \(\bar{c}_j \le 0\), rješenje je optimalno
C) Ako je barem jedan RHS negativan, rješenje je optimalno
D) Optimalnost se ne može testirati bez grafičkog prikaza

Provjera odgovora

C; 2) B; 3) A; 4) B; 5) B; 6) B; 7) C; 8) B; 9) D; 10) A;
B; 12) A; 13) B; 14) B; 15) C; 16) B; 17) B; 18) A; 19) D; 20) B.

Korištena literatura

Anderson, D. R., Sweeney, D. J., Williams, T. A., Camm, J. D., & Cochran, J. J. (2012). Quantitative methods for business (12th ed.). South-Western Cengage Learning.

Dantzig, G. B. (1951). Maximization of a linear function of variables subject to linear inequalities. In T. C. Koopmans (Ed.), Activity analysis of production and allocation (pp. 339–347). Wiley. (Original work published 1947)

Dantzig, G. B. (1963). Linear programming and extensions. Princeton University Press.

Bland, R. G. (1977). New finite pivoting rules for linear programming. Mathematics of Operations Research, 2(2), 103–107.

Goldfarb, D., & Reid, J. (1977). A practicable steepest‐edge simplex algorithm. Mathematical Programming, 12, 361–371. (provjeriti točan raspon stranica ako je potrebno)

Forrest, J. J., & Goldfarb, D. (1992). Steepest‐edge simplex algorithms for linear programming. Mathematical Programming, 57(3), 341–374. (provjeriti točan raspon stranica ako je potrebno)

Harris, P. M. J. (1973). Pivot selection methods of the Devex LP code. Mathematical Programming, 5(1), 1–28. (provjeriti točan raspon stranica ako je potrebno)

Lemke, C. E. (1954). The dual method of solving the linear programming problem. Naval Research Logistics Quarterly, 1(1), 36–47.

Karmarkar, N. (1984). A new polynomial-time algorithm for linear programming. Combinatorica, 4(4), 373–395.

Paketi / dokumentacija (R)

lpSolve — CRAN dokumentacija/udžbenički primjeri.
linprog — CRAN dokumentacija (revidirani/simplex pristupi).
ROI.plugin.lpsolve — ROI sučelje za lpSolve.
ROI.plugin.highs / HiGHS — moderni solver (revidirani dualni Simplex, interior-point).

Simplex metoda

Katarina Kostelić

Intuicija Simplex metode

Provedba Simplex metode

Priprema

Simplex tablica

Odabir ulazne varijable

Kad odabrati pojedini pristup?

Odabir izlazne varijable

Kad odabrati pojedini pristup?

Update tablice

Standardni pristup (ručna verzija) – Gauss-Jordanova eliminacija

Računalni pristupi (optimizirani)

Revidirani Simplex (Revised Simplex)

Product Form of the Inverse (PFI)

LU dekompozicija bazne matrice

Paketi dostupni u R-u

Primjeri ručnog rješavanja i upotrebe paketa

1. primjer

2. primjer

3. primjer

4. primjer

4. primjer a)

4. primjer b)

5. primjer

6. primjer

7. primjer

8. primjer

9.primjer

10. primjer

11. primjer

Pitanja za ponavljanje

Provjera odgovora

Korištena literatura