使用 R 的 lpSolve
與 lpSolveAPI
套件解決線性規劃問題,並提供實際的應用範例。
在日常生活中有許多的問題都可以使用線性規劃來尋求最好的解決方案,甚至在商業管理領域中,它也被大量應用在降低成本、提升產值與營收的策略上,而一般人沒有修過這門課程,不知道原來有那麼好用的工具,這裡我們將以實際而且常見的範例,介紹如何使用 R 解決線性規劃的問題。
生產褲子與外套問題
線性規劃(linear programming)是數學上的一種最佳化問題,以下我們直接透過範例來了解什麼是線性規劃問題。
假設今天有一個成衣的經銷商,要求工廠生產褲子與外套,這兩種產品都是使用棉與聚酯纖維來製作的,工廠內現有的原料中,棉有 750 m2,而聚酯纖維則有 1000 m2,生產褲子與外套所需要的原料用量如下:
褲子 | 外套 | |
---|---|---|
棉 | 1 m2 | 1.5 m2 |
聚酯纖維 | 2 m2 | 1 m2 |
每件褲子的價格是 1100 元,而外套則是 800 元,現在的問題就是工廠應該生產多少件的褲子,以及多少件的外套,才能讓總銷售金額達到最高呢?
像這種在多種選擇之間,尋找最佳組合的問題,就是一個典型的線性規劃問題。
數學條件式
若要解決這種問題,首先要把問題以數學的條件式寫出來。
假設我們要生產的褲子與外套的件數分別為 x
與 y
:
y = 生產的外套件數
而由於原料有限,所以在生產的時候,使用的棉與聚酯纖維不可以超過現有的存量,因此就可以寫出兩條限制式:
2x + y ≤ 1000
另外生產的件數不可以是負的,所以再加上兩條限制式:
y ≥ 0
而我們要做的事情就是在這些限制條件下,找到最好的 x
與 y
,讓銷售金額(目標函數)最高:
使用 R 解線性規劃問題
若要使用 R 來解這種簡單的線性規劃問題,最簡單的方式就是使用 lpSolve
套件,使用前先安裝並載入之:
# 安裝 lpSolve 套件 install.packages("lpSolve") # 載入 lpSolve 套件 library(lpSolve)
接著將上面寫好的數學條件式,依據 lpSolve
的格式輸入進來:
# 銷售金額(目標函數) # 1100x + 800y f.obj <- c(1100, 800) # 條件限制 # x + 1.5y <= 750 # 2x + y <= 1000 # x >= 0 # y >= 0 f.con <- matrix(c(1, 1.5, 2, 1, 1, 0, 0, 1), nrow = 4, byrow = TRUE) f.dir <- c("<=", "<=", ">=", ">=") f.rhs <- c(750, 1000, 0, 0)
然後把這些條件放進 lp
函數,這樣它就會自動找出最佳的解了:
# 解線性規劃問題 result <- lp("max", f.obj, f.con, f.dir, f.rhs) result
Success: the objective function is 612500
這個輸出的值就是目標函數的最大值,也就是說最大的銷售金額是 612500 元。
接著將解出來的 x
與 y
顯示來:
# 輸出 x 與 y 的解 result$solution
[1] 375 250
所以我們若要讓銷售金額達到最大值,就要生產 375 件褲子,以及 250 件外套。
畫圖
這是將線性規劃問題的各條限制式,畫在座標軸上的圖:
圖上的四條線就是四個限制條件,而我們就是要在這四條線所包住的區域中,找到能夠使銷售金額達到最大值的組合,最後找出來的解就是 (375, 250) 這個點。
畫這張圖的 R 程式碼如下:
# 畫圖 label.points <- data.frame(x = c(375, 500, 0, 0), y = c(250, 0, 0, 500)) ggplot(data.frame(x=c(-100,1000), y = c(-100,1000)),aes(x,y)) + stat_function(fun=function(x) (750 - x)/1.5 , geom="line", aes(col='x + 1.5y = 750'), size = 2) + stat_function(fun=function(x) 1000 - 2*x, geom="line", aes(col='2x + y = 1000'), size = 2) + geom_vline(xintercept = 0, aes(col='x = 0'), size = 2) + geom_hline(yintercept = 0, aes(col='y = 0'), size = 2) + geom_polygon(data = label.points, alpha = 0.2) + geom_point(data = label.points, aes(x, y), size = 5) + annotate('text', x = 375, y = 250, label="(375, 250)", size=5, hjust = -0.1, vjust = -0.5) + annotate('text', x = 500, y = 0, label="(500, 0)", size=5, hjust = -0.1, vjust = -0.5) + annotate('text', x = 0, y = 500, label="(0, 500)", size=5, hjust = -0.1, vjust = -0.5) + ylim(-100, 1200)
lpSolveAPI 套件
除了 lpSolve
之外,也可以使用 lpSolveAPI
套件,它的功能更多,但是語法也比較複雜,以下是解這個線性規劃問題的範例,首先安裝並載入 lpSolveAPI
套件:
# 安裝 lpSolveAPI 套件 install.packages("lpSolveAPI") # 載入 lpSolveAPI 套件 library(lpSolveAPI)
一開始要先建立一個線性規劃模型的物件:
# 建立線性規劃模型 # 包含兩個變數、兩個條件式 my.lp <- make.lp(2, 2)
接著指定限制條件與目標函數,而 lpSolveAPI
的模型預設就有限制變數要大於或等於零,所以我們只需要輸入其餘兩條條件就夠了:
# 條件限制 # x + 1.5y <= 750 # 2x + y <= 1000 set.column(my.lp, 1, c(1, 2)) set.column(my.lp, 2, c(1.5, 1)) set.constr.type(my.lp, c("<=", "<=")) set.rhs(my.lp, c(750, 1000)) # 銷售金額(目標函數) # 1100x + 800y lp.control(my.lp, sense="max") set.objfn(my.lp, c(1100, 800))
查看一下建立好的模型:
my.lp
Model name: C1 C2 Maximize 1100 800 R1 1 1.5 <= 750 R2 2 1 <= 1000 Kind Std Std Type Real Real Upper Inf Inf Lower 0 0
呼叫 solve
來解這個解線性規劃問題:
# 解線性規劃問題 solve(my.lp)
[1] 0
解完之後,輸出目標函數的最大值:
# 目標函數的最大值 get.objective(my.lp)
[1] 612500
輸出解出來的 x 與 y 值:
# 解出來的 x 與 y get.variables(my.lp)
[1] 375 250