用R语言做线性规划求解

简介

本文将对《数据、模型与决策:基于Excel的建模和商务应用》这本书的第10章规划求解的题目使用R语言来进行线性规划求解。

用到了两个用于做线性规划求解的软件包,分别是:Rglpk、Rsymphony。

例题一:生产计划问题

在本例中,将使用Rglpk包来进行规划求解。

Rglpk简介

安装Rglpk

install.packages("Rglpk")

Rglpk用法

Rglpk_solve_LP((obj, mat, dir, rhs, bounds = NULL, types = NULL, max = FALSE,control = list(), ...))

参数详解

参数 作用
obj 规划目标系数
mat 约束向量矩阵
dir 约束方向向量,有’>’、’<’、’=’构成
rhs 约束值
bounds 上下限的约束,默认0到INF
type 限定目标变量的类型,’B’指的是0-1规划,’C’代表连续,’I’代表整数,默认是’C’
control 包含四个参数verbose、presolve、tmlimit、canonicalizestatus。

题目

一个工厂有车床、刨床、钻床和铣床四种设备。生产A、B、C、D、E五中产品。每种设备每天生产时间为8小时,每年工作日为250天。各种设备台数,全年能力(可用工时),每种设备生产一件分别占用这四种设备的工时(单位:小时),五种产品可以获得的利润(单位:元/件),如下表所示:

生产计划问题的数据
设备类型 产品A 产品B 产品C 产品D 产品E 设备能力
车床 0.23 0.44 0.17 0.08 0.36 24000
刨床 0.13 0 0.2 0.37 0.19 22000
钻床 0 0.25 0.34 0 0.18 16000
铣床 0.55 0.72 0 0.61 0 12000
利润 123 94 105 132 118

现在我们要确定这五种产品的生产数量,使得占用的设备工时不超过各种设备的能力,同时使利润最大。

数学模型

$$ max(z) = 123x1+94x2+105x3+132x4+118x5
\left{\begin{array}{rl} \quad 0.23x
1+0.44x2+0.17x3+0.08x4+0.36x5 \leq 24000
\quad 0.13x1 \qquad\qquad +0.20x3+0.37x4+0.19x5 \leq 22000
\quad 0.25x2+0.34x3 \qquad\qquad +0.18x5 \leq 16000
\quad 0.55x
1+0.72x2 \qquad\qquad +0.61x4 \qquad\qquad \leq 12000
\quad x1,x2,x3,x4,x_5 \geq0 ,且为正整数
\end{array} \right.
$$

代码实现

library(Rglpk)

obj <- c(123, 94, 105, 132, 118)
row1 <- c(0.23, 0.44, 0.17, 0.08, 0.36)
row2 <- c(0.13, 0.00, 0.20, 0.37, 0.19)
row3 <- c(0.00, 0.25, 0.34, 0.00, 0.18)
row4 <- c(0.55, 0.72, 0.00, 0.61, 0.00)
mat <- as.matrix(rbind(row1,row2,row3,row4))
mat
##      [,1] [,2] [,3] [,4] [,5]
## row1 0.23 0.44 0.17 0.08 0.36
## row2 0.13 0.00 0.20 0.37 0.19
## row3 0.00 0.25 0.34 0.00 0.18
## row4 0.55 0.72 0.00 0.61 0.00
dir <- c("<=", "<=", "<=", "<=")
rhs <- c(24000, 22000, 16000, 12000)
max <- TRUE

rst <- Rglpk_solve_LP(obj, mat, dir, rhs, max = max, types = "I")
rst
## $optimum
## [1] 10872517
## 
## $solution
## [1]     0     0 18771 19672 53431
## 
## $status
## [1] 0
## 
## $solution_dual
## [1] NA
## 
## $auxiliary
## $auxiliary$primal
## [1] 23999.99 21184.73 15999.72 11999.92
## 
## $auxiliary$dual
## [1] NA

得到结果是:

$$ \left{\begin{array}{lr} x1 = 0
x
2 = 0
x3 = 18771
x
4 = 19672
x_5 = 53431
\end{array} \right.
$$

利润最大值为:10872517。

四台设备的使用工时分别为:23999.99, 21184.73, 15999.72, 11999.92

例题二:配料问题

本例中将会尝试使用Rsymphony包来进行规划求解。

Rsymphony简介

安装

install.package("Rsymphony")

使用方法

Rsymphony,混合整数线性规划SYMPHONY 求解器,其中主函数有:

Rsymphony_solve_LP(obj, mat, dir, rhs, bounds = NULL, 
                   types = NULL, max = FALSE, verbosity = -2, 
                   time_limit = -1, node_limit = -1, gap_limit = -1, 
                   first_feasible = FALSE,
                   write_lp = FALSE, write_mps = FALSE)

参数详解

主要参数 作用
obj 规划目标系数
mat 约束向量矩阵
dir 约束方向向量,有’>’、’<’、’=’构成
rhs 约束值
bounds 上下限的约束,默认0到INF
type 限定目标变量的类型,’B’指的是0-1规划,’C’代表连
max 逻辑值,T为求最大值,F为最小值

题目

化肥厂用四种原料A、B、C、D混合成复合肥料,这四种原料的单价以及复合肥料(M)所要求的氮、磷、钾的最低百分量如下表:

复合肥料配料问题数据

A B C D M
氮N 0.3 0.15 0 0.15 0.15
磷P 0.1 0 0.25 0.15 0.15
钾K 0 0.2 0.15 0.15 0.1
单价(元/吨) 2200 1800 2400 2700

要求配1000吨复合肥料,并假定在配置过程中物料没有损耗。

求使得总成本最低的配料方案。

数学模型

$$ min(z) = 2200x1 + 1800x2 + 2400x3 + 2700x4
\left{\begin{array}{lr} 0.30x1 + 0.15x2 + 0.00x3 + 0.15x4 \geq 150
0.10x1 + 0.00x2 + 0.25x3 + 0.15x4 \geq 150
0.00x1 + 0.20x2 + 0.15x3 + 0.15x4 \geq 100
x1 + x2 + x3 + x4 = 1000
x1,x2,x3,x4 \geq 0
\end{array} \right.
$$

书中$min(z)$那条公式的系数有误。

代码实现

library(Rsymphony)

obj <- c(2200, 1800, 2400, 2700)
row1 <- c(0.3, 0.15, 0, 0.15)
row2 <- c(0.1, 0.00, 0.25, 0.15)
row3 <- c(0, 0.2, 0.15, 0.15)
row4 <- c(1,1,1,1)

mat <- as.matrix(rbind(row1, row2, row3,row4))

dir <- c(">=", ">=", ">=","==")
rhs <- c(0.15*1000, 0.15*1000, 0.10*1000, 1000)
max <- F

Rsymphony_solve_LP(obj, mat, dir, rhs, max = F)
## $solution
## [1] 375 125 375 125
## 
## $objval
## [1] 2287500
## 
## $status
## TM_OPTIMAL_SOLUTION_FOUND 
##                         0
# 再用例子1中的Rglpk包验证一下答案
Rglpk_solve_LP(obj, mat, dir, rhs, max = F)
## $optimum
## [1] 2287500
## 
## $solution
## [1] 375 125 375 125
## 
## $status
## [1] 0
## 
## $solution_dual
## [1] 0 0 0 0
## 
## $auxiliary
## $auxiliary$primal
## [1]  150  150  100 1000
## 
## $auxiliary$dual
## [1]  7833.333  8750.000  8250.000 -1025.000

可见两个包计出来的结果都是一样的,也跟书中的结果一致。

运输问题有一个特点:只要供应量和需求量都是整数,那么最优解中的决策变量一定是整数,无需将决策定义成整数变量。

例题三:背包问题

题目

一艘货船最大装载重量为5000千克,现有A, B, C, D, E, F六种货物待装运,每种货物单件的价值和重量如下表:

| 物货名称 | A | B | C | D | E | F | | :-: | -: | -: | -: | -: | -: | -: | | 重量(公斤/件) | 320 | 420 | 530 | 550 | 590 | 640 | | 价值(万元/件) | 2.75 | 3.22 | 4.55 | 4.73 | 5.01 | 5.5 |

数学模型

$$ max(z) = 2.75x1 + 3.22x2 + 4.55x3 + 4.73x4 + 5.01x5 + 5.5x6
\left{\begin{array}{lr} 320x1 + 420x2 + 530x3 + 550x4 + 590x5 + 640x6 \leq 5000
\end{array} \right.
$$

代码实现

obj <- c(2.75, 3.22, 4.55, 4.73, 5.01, 5.5)
row1 <- c(320, 420, 530, 550, 590, 640)

mat <- as.matrix(rbind(row1))

dir <- c("<=")
rhs <- c(5000)
max <- T
types <- c("I")

Rsymphony_solve_LP(obj, mat, dir, rhs, max = max, types = types)
## $solution
## [1] 2 0 2 6 0 0
## 
## $objval
## [1] 42.98
## 
## $status
## TM_OPTIMAL_SOLUTION_FOUND 
##                         0

例题四:物流配送问题

某种产品从两个生产地A1、A2运往三个需求地B1、B2、B3。各生产低端的生产量、各需求地的需求量、每个生产地到每个需求地每吨产品的运输价格如下表:

物流配送问题的数据

| 运价(元/吨) | B1 | B2 | B3 |
| :-: | -: | -: | -: |
| A1 | 12 | 13 | 21 |
| A2 | 14 | 17 | 8 |

总运费最低的配送方案。

变量设置

| 运量 | B1 | B2 | B3 | 生产量(吨) | | :-: | :-: | :-: | :-: | -: | | A1 | x11 | x12 | x13 | 520 | | A2 | x21 | x22 | x23 | 480 | | 需求量(吨) | 200 | 400 | 400 | 1000 |

两个产地的总产量和三个需求地的总需求量是相等的,这个问题称为供求平衡的物流配送问题

数学模型

$$ min(z) = 12x{11} + 13x{12} + 21x{13} + 14x{21} + 17x{22} + 8x{23}
\left{\begin{array}{rl} x{11} + x{12} + x{13} = 520
x
{21} + x{22} + x{23} = 480
x{11} + x{21} = 200
x{12} + x{22} = 400
x{13} + x{23} = 400
\end{array} \right.
$$

代码实现

obj <- c(12, 13, 21, 14, 17, 8)
row1 <- c(1, 1, 1, 0, 0, 0)
row2 <- c(0, 0, 0, 1, 1, 1)
row3 <- c(1, 0, 0, 1, 0, 0)
row4 <- c(0, 1, 0, 0, 1, 0)
row5 <- c(0, 0, 1, 0, 0, 1)

mat <- as.matrix(rbind(row1,row2,row3,row4,row5))

dir <- c("==", "==", "==", "==", "==")
rhs <- c(520, 480, 200, 400, 400)
max <- T
types <- c("I")

Rsymphony_solve_LP(obj, mat, dir, rhs, max = max, types = types)
## $solution
## [1] 120   0 400  80 400   0
## 
## $objval
## [1] 17760
## 
## $status
## TM_OPTIMAL_SOLUTION_FOUND 
##                         0

例题五:[0, 1]选择问题

题目

一家控股公司要在下属的五家子公司中选择三家准备上市。这五家子公司的资产、负债和税后利润如下表:

公司选择问题的数据

| 子公司 | A | B | C | D | E | | :-: | -: | -: | -: | -: | -: | | 资产(亿元) | 3.48 | 5.62 | 7.33 | 6.27 | 2.14 | | 负债(亿元) | 1.28 | 2.53 | 1.02 | 3.55 | 0.53 | | 税后利润(万元) | 5400 | 2300 | 4600 | 3300 | 980 |

要求所选的三家子公司的总资产不低于10亿元,负债不超过5亿元,并使得新组建的公司税后利润最大。

设有5个c(0, 1)的变量,x1,x2,x3,x4,x5,取值如下:

$$ x_i=\left{\begin{array}{rl} 0 & 不选择的子公司
1 & 选择的子公司
\end{array} \right.
$$

数学模型

$$ min(z) = 5400x{1} + 2300x{2} + 4600x{3} + 3300x{4} + 980x{5}
\left{\begin{array}{rl} & 3.48x
1 + 5.62x2 + 7.33x3 + 6.27x4 + 2.14x5 \geq 10
& 1.28x1 + 2.53x2 + 1.02x3 + 3.55x4 + 0.53x5 \leq 5
& x
1 + x2 + x3 + x4 + x5 = 3
& x1, x2, x3, x4, x_5 为 0或1 \end{array} \right.
$$

代码实现

obj <- c(5400, 2300, 4600, 3300, 980)
row1 <- c(3.48, 5.62, 7.33, 6.27, 2.14)
row2 <- c(1.28, 2.53, 1.02, 3.55, 0.53)
nx <- length(obj)
mat.tmp <- matrix(rep.int(0, nx^2), nrow = nx)
for (i in 1:nx) {
  for (j in 1:nx) {
    if (i == j) {
      mat.tmp[i, j] <- 1
    }
  }
} 

mat <- as.matrix(rbind(row1,row2, mat.tmp))

dir <- c(">=", "<=", rep.int("<=", 5))
rhs <- c(10, 5, 1, 1, 1, 1, 1)
max <- TRUE
types <- c("I")

Rsymphony_solve_LP(obj, mat, dir, rhs, max = max, types = types)
## $solution
## [1] 1 1 1 0 0
## 
## $objval
## [1] 12300
## 
## $status
## TM_OPTIMAL_SOLUTION_FOUND 
##                         0

输出的结果跟书中一致。

例题六:指派问题

指派问题简介

指派问题:有n项任务由n个人完成,每项任务交给一个人,每人都有一项任务。由i个人完成j项任务的成本(或效益)为$c_{ij}$。求使总成本最小(或总效益最大)的分配方案。

$$ x_{ij}\left{\begin{array}{rl} 0 & 第i个人不是从事第j项任务
1 & 第i个人被指派完成第j项任务
\end{array} \right.
$$

这也是一个c(0, 1)规划问题,这个问题的0-1线性规划模型如下:

$$ max(z)=\sum{i=1}^{n}\sum{j=1}^{n}{c{ij}x{ij}}
\sum{i=1}^n{x{ij}}=1 , j = 1,2,\dots,n
\sum{j=1}^n{x{ij}}=1 , i = 1,2,\dots,n
$$

这类问题称为指派问题(assignment problem)。

题目

市政府有四项市政建设工程招标,经过初选,四家建设公司最后参与这四项工程竞标。每家公司对每个工程的报价如下:

各公司的报价

公司 工程A 工程B 工程C 工程D
甲公司 920 480 650 340
乙公司 870 510 700 350
丙公司 880 500 720 400
丁公司 930 490 680 410

市政府规定,每项工程只能有一家公司中标,每个公司只能承担一项工程。

求总价最低的决标方案。

数学模型

$$ max(z)=920x{11} + 480x{12} + 650x{13} + 340x{14}
+ 870x{21} + 510x{22} + 700x{23} + 350x{24}
+ 880x{31} + 500x{32} + 720x{33} + 400x{34}
+ 930x{41} + 490x{42} + 680x{43} + 410x{44}
\left{\begin{array}{rl} x{11} + x{12} + x{13} + x{14} = 1
x{21} + x{22} + x{23} + x{24} = 1
x{31} + x{32} + x{33} + x{34} = 1
x{41} + x{42} + x{43} + x{44} = 1
x{11} + x{21} + x{31} + x{41} = 1
x{12} + x{22} + x{32} + x{42} = 1
x{13} + x{23} + x{33} + x{43} = 1
x{14} + x{24} + x{34} + x{44} = 1
\end{array} \right.
$$

前四条约束和后四条约束分别是每家公司只能承接一项工程每项工程只能由一家公司来承接

代码实现

obj <- c(920, 480, 650, 340, 
         870, 510, 700, 350,
         880, 500, 720, 400,
         930, 490, 680, 410
         )
nx <- length(obj) 
mat.length <- nx*8
mat.tmp <- matrix(rep.int(0, mat.length), ncol = 16)

# 生成前四条约束矩阵
j = 1:4
for(i in 1:4) {
    mat.tmp[i, j] = 1
    j = j + 4
    if (j[4] > 16) {
      j = 1:4
    }
}
# 生成后四条约束矩阵
j = c(1,5,9,13)
for(i in 5:8) {
    mat.tmp[i, j] = 1
    j = j + 1
    if (j[4] > 16) {
      j = c(1,5,9,13)
    }
}

mat <- mat.tmp
mat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,]    1    1    1    1    0    0    0    0    0     0     0     0     0
## [2,]    0    0    0    0    1    1    1    1    0     0     0     0     0
## [3,]    0    0    0    0    0    0    0    0    1     1     1     1     0
## [4,]    0    0    0    0    0    0    0    0    0     0     0     0     1
## [5,]    1    0    0    0    1    0    0    0    1     0     0     0     1
## [6,]    0    1    0    0    0    1    0    0    0     1     0     0     0
## [7,]    0    0    1    0    0    0    1    0    0     0     1     0     0
## [8,]    0    0    0    1    0    0    0    1    0     0     0     1     0
##      [,14] [,15] [,16]
## [1,]     0     0     0
## [2,]     0     0     0
## [3,]     0     0     0
## [4,]     1     1     1
## [5,]     0     0     0
## [6,]     1     0     0
## [7,]     0     1     0
## [8,]     0     0     1
dir <- rep.int("==", 8)
rhs <- rep.int(1, 8)
max <- FALSE

rst <- Rsymphony_solve_LP(obj, mat, dir, rhs, max = max)
matrix(rst$solution, ncol = 4)
##      [,1] [,2] [,3] [,4]
## [1,]    0    0    1    0
## [2,]    0    0    0    1
## [3,]    1    0    0    0
## [4,]    0    1    0    0
rst
## $solution
##  [1] 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0
## 
## $objval
## [1] 2370
## 
## $status
## TM_OPTIMAL_SOLUTION_FOUND 
##                         0

输出结果与书中结果一致,总报价为2370。书中数学模型中的第一行max z应该为min z,题目要求求最小值,应属笔误,因为结果是对的。

指派问题也和运输问题有类似的性质,只要约束条件的右边常数都是等于1,最优解中的决策变量自然就等于0或1,不比将决策变量定义成0-1变量,如例子六中并没有像例子五那样控制变量小于等于1。

线性规划问题总结

线性规划两条基本定理:

  • 如果线性规划的可行域不是空集,则线性规划的可行域是一个凸集。
  • 线性规划如果有最优解,最优解至少在可行域的一个极点上。

线性规划问题的基本概念:

  • 决策变量
  • 线性规划的解
  • 目标函数
  • 约束条件
  • 右边常数
  • 变量的符号约束(决策变量>=0或者<=0或没有符号限制)
  • 可行解
  • 可行域
  • 目标函数等值线
  • 目标函数的梯度方向
  • 最优解

求解线性规划的算法——单纯形法

线性规模模型的目标函数必须是变量的线性函数,约束条件必须是变量的现行等式或不等式。

下面这样的就不是线性规划问题:

$$ min(z)=3x1^2+2x1x2
\left{\begin{array}{rl} 2x
1+\frac{x2+x3}{x1} \leq 8
|x
1+x2|+4x3 \leq 9
x1,x2,x_3 \geq 0
\end{array} \right.
$$

这是因为:

  1. 目标函数中出现变量的平方以及两个变量的乘积
  2. 第一个约束条件的分母中出现变量
  3. 第二个约束条件中出现非线性函数

另外一些(例子一):

  • 影子价格越大,那么该变量的变化对结果影响越大。
  • 如果设备能力(例子一)有剩余(松弛值>0)该影子价格一定等于0。
  • 最优生产计划中安排生产的产品他们的机会成本都等于利润(即差额成本=0)
  • 机会成本大于利润的产品,不安排生产(即产量为0)

线性规划模型是一种典型的单目标决策模型

Excel打开规划求解

文件 -> 选项 -> 加载项 -> 转到 -> 规划求解加载项 -> 确定

也可以顺带把分析工具库也选上。

参考资料

  1. http://blog.csdn.net/qq_27755195/article/details/53759566
  2. 蒋绍忠.数据、模型与决策:基于Excel的建模和商务应用[M].北京:北京大学出版社,2010:355-386.