Wstęp
Wyobraźmy sobie, że jesteśmy właścicielami fabryki, w której sprzedawane są dwa rodzaje mebli: krzesła i fotele. Mamy określoną ilość materiału i ograniczoną przestrzeń w magazynie. Ile krzeseł i ile foteli powinna produkować nasza fabryka, tak aby zyski były maksymalne? Aby odpowiedzieć na to pytanie wykorzystamy zagadnienie programowania liniowego. Programowanie liniowe to metoda rozwiązania problemu masymalizacji lub minimalizacji liniowej funkcji celu \(f (x_1, x_2, ..., x_n)\) przy pewnej ilości ograniczeń nałożonych na zmienne \(x_1, x_2, ... x_n\). Zagadnienie to zostało omówione na przykładzie przedstawionego zadania optymalizacyjnego.
Zadanie:
Fabryka „Wygodnie” produkuje dwa rodzaje mebli: krzesła i fotele. Do produkcji wykorzystuje drewno, materiał i metal (gwoździe, śrubki itp.). Stan zasobów w magazynie przedstawiony jest w Tabeli 1.
Tabela 1
| Tworzywo | Zasoby_kg |
|---|---|
| Drewno | 7500 |
| Materiał | 3000 |
| Metal | 500 |
Ilość zasobów potrzebnych do produkcji oraz zysk jednostkowy wytworzonego produktu został zaprezentowany w Tabeli 2.
Tabela 2
| x | Krzesło | Fotel |
|---|---|---|
| Drewno [w kg] | 3.0 | 1.0 |
| Materiał [w kg] | 0.5 | 1.0 |
| Metal [w kg] | 0.1 | 0.1 |
| Zysk jednostkowy produktu [zł za sztukę] | 200.0 | 100.0 |
Ze względu na pojemność magazynu firma może wyprodukować tylko 3000 krzeseł lub 4000 foteli. Ile krzeseł, a ile foteli powinna produkować firma, aby zmaksymalizować swój zysk?
Rozwiązanie:
Aby rozwiązać powyższe zadanie najpierw trzeba przedstawić problem za pomocą wyrażeń matematycznych:
1. Definicja zmiennych
Naszym celem jest znalezieniem optymalnych wielkości produkcji krzeseł i foteli, w związku z czym definujemy zmienne następująco:
x- liczba krzeseł
y- liczba foteli
2. Układ ograniczeń
Poniższe nierównoście przedstawiają opisane w zadaniu ograniczenia nałożone na wielkości produkcji:
Ograniczenia ze względu na dostępność zasobów:
\(3x+y \leq 7 500\)
\(0,5x+y \leq 3 000\)
\(0,1x+0,1y \leq 500\)
Ograniczenia ze względu na wielkość magazynu:
\(x/3000+y/4000 \leq 1\)
Ograniczenia brzegowe:
\(x \leq 0\)
\(y \leq 0\)
3. Funkcja celu
W tym przypadku funkcją celu jest funkcja zysku określona wzorem:
\(Z=200x+100y ~~-> max\)
Rozwiązanie numeryczne
Rozwiązanie przy użyciu pakietu lpSolve
Zacznijmy od rozwiązania zadania z wykorzystaniem pakietu lpSolve.
# do C przypisano współczynniki funkcji celu -> C
C <- c(200, 100)
# macierz ograniczeń
A <- matrix(c(3, 1,
0.5, 1,
0.1, 0.1,
4, 3), nrow=4, byrow=TRUE)
# prawe strony ogarniczeń
B <- c(7500, 3000, 500, 12000)
# znaki ograniczeń
constranints_direction <- c("<=", "<=", "<=", "<=")
# znalezienie optymalnego rozweiązania
optimum <- lp(direction="max",
objective.in = C,
const.mat = A,
const.dir = constranints_direction,
const.rhs = B,
all.int = T)
# Status: 0 = sukces, 2 = brak możliwości rozwiązania
print(optimum$status)## [1] 0#Wyświetl optymalne wartości dla x, y
best_sol <- optimum$solution
names(best_sol) <- c("x", "y")
print(best_sol)## x y
## 2100 1200# Sprawdź wartość funkcji celu w optymalnym punkcie
print(paste("Całkowity zysk = ", optimum$objval, sep=""))## [1] "Całkowity zysk = 540000"Rozwiązanie przy użyciu pakietu lpSolveAPI
Spróbujmy ponownie rozwiązać ten sam problem, tym razem korzystając z pakietu lpSolveAPI.
# Ustaw 4 ograniczenia i 3 zmienne decyzyjne
lprec <- make.lp(nrow = 4, ncol = 2)
# Ustaw typ problemu, który próbujemy rozwiązać
lp.control (lprec, sense = "max")## $anti.degen
## [1] "fixedvars" "stalling"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] 1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "maximize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"# Ustaw typ zmiennych decyzyjnych
set.type(lprec, 1:2, type=c("integer"))
# Ustaw wektor współczynników funkcji celu
set.objfn (lprec, C)
# Dodaj ograniczenia
add.constraint(lprec, A[1, ], "<=", B[1])
add.constraint(lprec, A[2, ], "<=", B[2])
add.constraint(lprec, A[3, ], "<=", B[3])
add.constraint(lprec, A[4, ], "<=", B[4])
# Rozwiąż problem
solve (lprec)## [1] 0# Uzyskaj wartości zmiennych decyzyjnych
a <- get.variables (lprec)
a## [1] 2100 1200# Uzyskaj wartość funkcji celu
get.objective (lprec)## [1] 540000Wizualizacja
Rozwiązanie zadania numerycznie to nie wszystko. Wiadomo, że ludzie są z natury wzrokowcami, więc musimy zadbać także o aspekt wizualny rozwiązania. W końcu może się zdarzyć, że szukając potencjalych inwestorów, którzy pomogą nam rozwijać naszą działalność, będziemy musieli przekonać ich, że wykorzystujemy optymalnie nasze zasoby. Przejdźmy więc do prezentacji naszego rozwiązania w sposób czytelny i zrozumiały nawet dla “matematycznych ignorantów”.
Opcja 1: Wykres przy użyciu grafiki wbudowanej
Zacznijmy od grafiki wbudowanej. Korzyścią tego typu rozwiązania jest brak konieczności wyszukiwania i wgrywania specjalnych pakietów graficznych. Jeśli zależy nam na szybkim i prostym wykresie możemy skorzystać z polecenia plot().
funkcja_liniowa <- function(x,a,b){
y <- a*x + b
return(y) }
x <- seq(0,8000,by=0.1)
y <- funkcja_liniowa(x,a = -0.5,b =3000)
y1 <- funkcja_liniowa(x,a = -4/3,b =4000)
y2 <- funkcja_liniowa(x,a = -3,b =7500)
y3 <- funkcja_liniowa(x,a = -1,b =5000)
# punkty końcowe na osi poziomej
xx = c(1, x, 8000)
yy = c(0, y, 0)
yy1 = c(0, y1, 0)
yy2 = c(0, y2, 0)
yy3 = c(0, y3, 0)
plot(1,type="n",xlim=c(0,7000),ylim=c(0,8000),xlab="liczba krzeseł",ylab="liczba foteli",xaxs="i",yaxs="i",
main="Wizulizacja rozwiązania programowania liniowego")
polygon(xx,yy,col=rgb(0, 0, 1, alpha = 0.2), border=rgb(0, 0, 1, alpha = 0.5))
polygon(xx,yy1,col=rgb(1, 0, 0.5, alpha = 0.2), border = rgb(1, 0, 0.5, alpha = 0.5))
polygon(xx,yy2,col=rgb(0., 0, 0.200, alpha = 0.2), border = rgb(0, 0, 0.200, alpha = 0.5))
polygon(xx,yy3,col=rgb(1, 0, 0, alpha = 0.2), border = rgb(1, 0, 0, alpha = 0.5))
abline(2000,-2,lty=4)
abline(4000,-2,lty=4)
abline(5400,-2,lty=4)
points(a[1],a[2],pch=19)
text(4/3*a[1],a[2],paste("(" , a[1], ",", a[2],")"), col = "gray0")
text(1/2*a[1],a[2],"zbiór rozwiązań", col = "gray0")
Opcja 2: Wykres przy użyciu pakietu ggplot2
Jeśli zależy nam, na czymś trochę bardziej skomplikowanym możemy skorzystać z pakietu ggplot2. Poniżej przedstawiono przykład wizualizacji, którą możemy w prosty i szybki sposób modyfikować i wykorzystywać przy nowych zadaniach.
# definiujemy funkcje dotyczące organiczeń
f1 = function(x) -1/2*x+3000
f2 = function(x) -4/3*x+4000
f3 = function(x) -3*x+7500
f4 = function(x) -x+5000
# definujemy funkcję celu przy dwóch ustalonych stałych (stałe można dobrać na końcu)
f = function(x) -2*x + 5400
ff = function(x) -2*x + 2800
# definiujemy zestaw kolorów, które wykorzystamy na wykresie
colo <- c("purple1", "royalblue4", "purple4", "purple3", "gray0")
# tworzymy wykres
ggplot(data.frame(x=c(0, 7000)), aes(x)) +
#zaciemniamy pola odpowiadające nierównościom
stat_function(fun=f1, geom="area",fill= colo[1],alpha=0.3, na.rm = TRUE) +
stat_function(fun=f2, geom="area",fill= colo[2],alpha=0.3, na.rm = TRUE) +
stat_function(fun=f3, geom="area",fill= colo[3],alpha=0.3, na.rm = TRUE) +
stat_function(fun=f4, geom="area",fill= colo[4],alpha=0.3, na.rm = TRUE) +
stat_function(fun=f1, geom="line",col= colo[1],alpha=1, na.rm = TRUE) +
stat_function(fun = f2, geom = "line", col = colo[2], alpha = 1, na.rm = TRUE) + #dodajemy linie
stat_function(fun = f3, geom = "line", col = colo[3], alpha = 1, na.rm = TRUE) +
stat_function(fun = f4, geom = "line", col = colo[4], alpha = 1, na.rm = TRUE) +
stat_function(fun = f, geom = "line", col = colo[5], alpha = 1, linetype= "dashed", na.rm = TRUE)+
stat_function(fun = ff, geom = "line", col = colo[5], alpha = 1, linetype= "dashed", na.rm = TRUE)+
geom_segment(mapping=aes(x=0,y=0,xend=6500, yend=0), arrow=arrow(), size=1) +
geom_segment(mapping=aes(x=0,y=0,xend=0, yend=8000), arrow=arrow(), size=1) +
ylab("Liczba foteli") + xlab ("Liczba krzeseł") +
ggtitle("Wizualizacja rozwiązania programowania liniowego") +
ylim(0,8000) +
xlim(0,6500) +
geom_text(aes(x=5/3*a[1],y=3/2*a[2],label=paste("Punkt optymalny", "(" , a[1], ",", a[2],")")), col=colo[5]) +
geom_point(aes(x=a[1],y=a[2]), col= colo[5]) +
geom_text(aes(x=1/2*a[1],y=a[2],label="zbiór rozwiązań"), col= colo[5]) +
theme( #określamy wygląd wykresu
panel.background = element_rect(fill = "white", color = "#cdcdcd", size = 2),
panel.grid.major = element_line(colour = "#cdcdcd", linetype = "dashed",size = 1),
plot.title = element_text(face= "bold", size = 15),
axis.title = element_text(size = 10)
)
Kilka słów od autorek
Gdy napotykamy problem maksymalizacji/minimalizacji funkcji liniowej przy nałożonych ograniczeniach, nie musimy popadać w rozpacz, bowiem istnieje prosta i skuteczna metoda poradzenia sobie z tym problemem. Mamy nadzieję, że po przeczytaniu tej pracy programowanie liniowe już nie będzie dla Ciebie zagadką. Jeśli do tego dysponujesz w miarę sprawnym komputerem i potrafisz “przeklikać”" kilka linijek w R, Twoje zadanie okaże się prostsze niż myślisz.
Jak można zauważyć, każda z przedstawionych metod doprowadziła nas do tego samego wyniku, możemy więc jednoznacznie stwierdzić, że rozwiązaniem naszego zadania jest punkt (2100, 1200). Obie metody rozwiązania numerycznego, a także obie wizualizacje okazały się więc równie skuteczne. To od Ciebie zależy, którą z nich wybierzesz… Tylko wybierz jedną, bo robienie kilku wersji jest jak widać bez sensu :)
| Źródła: |