R语言与优化模型(一):规划问题和运输问题
作者:鲁伟,热爱数据,坚信数据技术和代码改变世界。R语言和Python的忠实拥趸,为成为一名未来的数据科学家而奋斗终生。
个人公众号:数据科学家养成记 (微信ID:louwill12)
线性规划、整数规划和运输问题
虽说目前数学建模领域大家用的软件工具都是Matlab或者Lingo,但在个人偏好的驱使下还是决定用自己擅长的R语言来实现数学建模算法。就首先从优化模型中简单的线性规划、整数规划以及多目标规划开始。R语言在针对各类优化模型时都能快速方便的求解,对运输问题、生产计划问题、产销问题和旅行商问题等都有专门的R包来解决。
线性规划与整数规划的区别主要在于对决策变量的取值约束有所不同。线性规划的决策变量为正实数,而整数规划则要求决策变量为正整数。在R语言中,有众多相关的R包可以解决这两类问题,例如stat包中的optim、optimize函数,这里给大家推荐Rglpk包,Rglpk包用法简单,核心函数调用方便,对解决大型的线性规划和整数规划问题十分好用。
Rglpk包的核心函数为Rglpk_solve_LP
函数
,其用法如下:
Rglpk_solve_LP(obj,mat,dir,rhs,types=NULL,max=FALSE,bounds=NULL, verbose=FALSE)
其中obj为目标函数系数,mat为约束矩阵,dir为约束符号,rhs为约束向量,types为变量类型,可选类型有“B”“I”“C”,分别代表0-1变量、正整数和正实数,max为逻辑参数,当其取TRUE 时求目标函数最大值,反之为最小值,bounds为决策变量的额外约束,verbose表示是否输出中间过程的的控制参数,默认为FALSE。
看一个求解线性规划的例子:
library(Rglpk)
obj<-c(10,15,12)#目标函数系数
mat<-matrix(c(5,-5,2,3,6,1,1,15,1),nrow=3)#约束矩阵
dir<-rep("<=",3)#约束符号
rhs<-c(9,15,5)#约束向量
Rglpk_solve_LP(obj,mat,dir,rhs,max=TRUE)#求解
$optimum
[1] 42
$solution
[1] 0.200000 2.666667 0.000000
求解过程清晰明了,比起专业的线性规划软件Lingo有过之而无不及。
再看一例整数规划(0-1规划):
obj<-c(4,3,2)
mat<-matrix(c(2,4,0,-5,1,1,3,3,1),nrow=3)
dir<-c("<=",">=",">=")
rhs<-c(4,3,1)
types<-c("B","B","B")#指定变量类型,模型为0-1规划
Rglpk_solve_LP(obj,mat,dir,rhs,types,max=FALSE)
$optimum
[1] 2
$solution
[1] 0 0 1
整数规划(0-1规划)和线性规划并无本质区别,所以在Rglpk中其实现函数一样,只是依靠types参数来控制两个模型的区别。
运输问题为常见的线性规划问题,但Rglpk包并没有普适性。R针对运输问题、生产计划等问题有专门的R包lpSovle,其核心函数 lp.transport 用法如下:
lp.transport(costs.mat,direction="min",row.signs,row.rhs,
col.signs,col.rhs,presolve=0,compute.sense=0,integers=1:(nc*nr))
其中cost.mat为运费矩阵,direction参数决定取最大值还是最小值,row.sign(产量约束符号)和row.rhs(产量约束向量)构成产量约束条件;col.sign(销量约束符号)和col.rhs(销量约束向量)构成销量约束条件,其余参数可忽略。
例:使用lpSovle包求解6个生产点和8个销售点的最小费用运输问题。
library(lpSolve)
costs<-matrix(c(6,4,5,7,2,5,2,9,2,6,3,5,6,5,1,7,9,
2,7,3,9,3,5,2,4,8,7,9,7,8,2,5,4,2,2,1,5,8,3,7,6,4,
9,2,3,1,5,3),nrow=6)#运费矩阵
row.signs<-rep("<=",6)#产量约束符号
row.rhs<-c(60,55,51,43,41,52)#产量约束向量
col.signs<-rep("=",8)#销量约束符号
col.rhs<-c(35,37,22,32,41,32,43,38)#销量约束向量
res<-lp.transport(costs,"min",row.signs,row.rhs,col.signs,
col.rhs)
res
Success: the objective function is 664
res$solution
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0 19 0 0 41 0 0 0
[2,] 0 0 0 32 0 0 0 1
[3,] 0 12 0 0 0 0 39 0
[4,] 0 0 0 0 0 6 0 37
[5,] 35 6 0 0 0 0 0 0
[6,] 0 0 22 0 0 26 4 0
由R输出结果可知6个生产点8个销售点的最小运费为664,相应的运送方案为A1→B2:19个单位,A1→B5:41个单位,A2→B4:32个单位,A2→B8:1个单位,A3→B2:12个单位,A3→B7:39个单位,A4→B8:37个单位,A5→B1:35个单位,A5→B2:6个单位,A6→B3:22单位,A6→B6:26个单位,A6→B7:4个单位。
非线性规划与多目标规划
与线性规划不同的是,非线性规划要求目标函数或约束条件中含有非线性函数。相应的求解这类问题就要用到非线性规划的方法。约束条件或者目标函数的放宽使得规划模型更具普适性,但也增加了问题求解的难度。对于简单的非线性规划问题,R语言中stat包即可求解。在这里我们给大家介绍R语言中求解非线性规划更为专业的Rdonlp2包。
Rdonlp2在求解非线性规划问题上功能十分强大,用户可自行安装并查阅帮助文档:
library(Rdonlp2)
Rdonlp2包的核心函数为 donlp2,可以快速求解非线性规划的最值。其用法如下:
donlp2(par,fn,
par.upper=rep(+Inf,length(par)),
par.lower=rep(+Inf,length(par)),
A=NULL,
lin.upper=rep(+Inf,length(par)),
lin.lower=rep(-Inf,length(par)),
nlin=list(),
nlin.upper=rep(+Inf,length(nlin)),
nlin.lower=rep(-Inf,length(nlin)),
control=donlp2.control(),
control.fun=function(lst){return(TRUE)},
env=.GlobalEnv,
name="Rdonlp2"
)
对该函数众多参数进行解释:
par : 初始迭代向量
fn : 连续型函数
par.upper、par.lower : 自变量的上下限
A:线性约束矩阵
lin.upper、lin.lower : 线性约束条件的上下界
nlin : 非线性约束条件函数列表
nlin.upper、nlin.lower : 非线性约束条件的上下界
其余参数可忽略。
先看一例donlp2函数求解非线性规划的例子:
用Rdonlp2包编写非线性规划计算命令:
library(Rdonlp2)
p=c(0,0) #迭代初始值
par.l=c(-100,100);par.u=c(100,100)#自变量定义域约束
fn=function(x){
exp(x[1])*(4*x[1]^2+2*x[2]^2+4*x[1]*x[2]+2*x[2]+1)
} #目标函数
A=matrix(c(1,1,3,-1),2,byrow=TRUE)
lin.l=c(2,1);lin.u=c(+Inf,3) #线性约束
nlcon1=function(x){
-x[1]*x[2]
} #非线性约束
donlp2(p,fn,par.u,par.l,A,lin.l=lin.l,lin.u=lin.u,
nlin=list(nlcon1))
$message
[1] "KT-conditions satisfied, no further correction computed"
$par
[1] 33.66667 100.00000
$gradf
[1] 1.625065e+19 2.243635e+17
相应运行结果包括中间过程、算法参数、运行时间等数据,这里就不全部列出。
再看一例利用Rdonlp2包求解二元rastrigin函数最小值的问题。该函数因为有众多极值点而仅有一个最值点对各类优化算法有一定的欺骗性,但用其来测试各类优化算法效果显著。该二元rastrigin函数为:
先作出该函数的三维图像:
a<-5
x<-seq(-a,a,0.01)
y<-seq(-a,a,0.01)
f<-function(x,y){
x^2-10*sin(2*pi*x)+y^2-10*cos(2*pi*y)+20
}
z<-outer(x,y,f)
image(x,y,z,col=heat.colors(24))
library(rgl)
zorder<-rank(z)
persp3d(x,y,z,col=rainbow(as.integer(max(zorder)))[zorder])
fn<-function(x){
f<-sum(x^2*cos(2*pi*x)+10)
}
par<-c(-500,500)
ret<-donlp2(par,fn)
ret$f
[1] 20
ret$par
[1] -2.654965e-06 2.654965e-06
由输出结果可知该函数最小值为20 。所以Rdonlp2求解非线性规划也是非常方便的。当然了,还有很多更为一般的非线性规划问题是Rdonlp2包也不能解决的,这时候可能会向遗传算法、模拟退火这些启发式算法求助啦。这里且不作讨论。
在许多实际问题中,衡量一个方案的好坏标准往往不止一个,例如制造一枚导弹,既要射程远又要最省燃料,还得精度高。这一类问题就不是单一的目标规划问题了,我们称之为多目标规划问题。
求解多目标规划问题通常有主要目标法、分层序列法和线性加权求和法等方法。其中主要目标法通过确定一个主要目标,将多目标优化问题转化为线性或非线性规划问题;分层序列法是将目标中多个目标按照其重要程度排一个次序而逐步求解的过程;线性加权法则是对各目标赋予一个权数,加权求和得到一个新的目标函数而进行求解的过程。
R语言中多目标规划问题的求解通常用goalprog包,其核心函数为llgp:
llgp(coefficient,targets,achievements,maxiter=1000,verbose=FALSE)
其中coefficient为约束系数矩阵,targets为系数矩阵约束向量,achievements是包含objective、priority、p、n的目标函数,其中objective表示第几对偏差变量,priority表示该偏差变量的优先级,p和n分别为正负偏差变量的权系数。
看一道多目标规划问题计算实例:(该题来自钱颂迪《运筹学》教材4.4节例5)
llgp函数计算过程如下:
library(goalprog)
coefficients<-matrix(c(1,1,5,1,1,0,3,1),4)
targets<-c(10,4,56,12)
achievements<-data.frame(objective=1:4,priority=c(1,1,3,4),
p=c(2,3,0,1),n=c(0,0,1,0))
soln<-llgp(coefficients,targets,achievements)
soln$converged
[1] TRUE
soln$out
Decision variables
X
X1 4.000000e+00
X2 6.000000e+00
Summary of objectives
Objective Over Under Target
G1 1.000000e+01 0.000000e+00 0.000000e+00 1.000000e+01
G2 4.000000e+00 0.000000e+00 0.000000e+00 4.000000e+00
G3 3.800000e+01 0.000000e+00 1.800000e+01 5.600000e+01
G4 1.000000e+01 0.000000e+00 2.000000e+00 1.200000e+01
Achievement function
A
P1 0.000000e+00
P2 0.000000e+00
P3 1.800000e+01
P4 0.000000e+00
由输出的结果我们可以看到x1、x2最优值分别为4和6。
由该例我们可以看出,R语言求解多目标规划问题只需按照函数格式将规划的目标函数、约束条件套入格式即可,过程非常简单。
参考文献:
魏太云.R软件与最优化[C]//中国r语言会议.2008.
往期推送
R语言交互式绘制杭州市地图:leafletCN包简介
kaggle:NBA球员投篮数据分析与可视化
R语言数据分析练手小项目:杭州二手房数据分析
如何像一个机器学习老司机一样跟别人解释SVM算法?
人工神经网络算法及其简易R实现
RStudio|用R Markdown生成你的R语言数据分析报告
如何使用reshape/reshape2使劲揉你的数据
R语言爬虫|15行代码教你抓取拉勾网招聘信息
用数据分析告诉你数据分析师能挣多少钱
R语言典型相关分析:NBA球员身体素质与统计数据关联性
为什么R语言是当今最值得学习的数据科学语言
R语言兵器谱:数据科学家的十八般武艺
R语言爬虫系列1|HTML基础与R语言解析
R语言爬虫系列2|XML&XPath表达式与R爬虫应用
R语言爬虫系列3|HTTP协议
R语言爬虫系列4|AJAX与动态网页介绍
R语言爬虫系列5|正则表达式与字符串处理函数
R语言爬虫系列6|动态数据抓取范例
公众号后台回复关键字即可学习
回复 爬虫 爬虫三大案例实战
回复 Python 1小时破冰入门回复 数据挖掘 R语言入门及数据挖掘
回复 人工智能 三个月入门人工智能
回复 数据分析师 数据分析师成长之路
回复 机器学习 机器学习的商业应用
回复 数据科学 数据科学实战
回复 常用算法 常用数据挖掘算法
R语言与优化模型(一):规划问题和运输问题相关推荐
- 055B ENMTools教程-基于R语言对MaxEnt模型优化-MaxEnt调参教程--更新日期2021-9
055B-1 视频附带资料下载和密码:软件-数据-文献下载-持续更新 055B-2 ENMTools软件下载安装 055B-3 R软件和工具包安装 055B-4 生物气候因子的精度说明与选择方法(理论 ...
- 基于R语言、MaxEnt模型融合技术的物种分布模拟、参数优化方法、结果分析制图与论文写作
详情链接 :基于R语言.MaxEnt模型融合技术的物种分布模拟.参数优化方法.结果分析制图与论文写作 内容介绍: 第一章 .理论篇 以问题导入的方式,深入掌握原理基础 : 什么是MaxEnt模型? ...
- R语言、MaxEnt模型融合技术的物种分布模拟、参数优化方法、结果分析制图与论文写作
基于R语言.MaxEnt模型融合技术的物种分布模拟.参数优化方法.结果分析制图与论文写作技术应用 第一章.理论篇以问题导入的方式,深入掌握原理基础 什么是MaxEnt模型? MaxEnt模型的原理是什 ...
- ENMTools教程-基于R语言对MaxEnt模型优化-MaxEnt调参教程介绍
MaxEnt3.4.4软件下载网盘: http://lucky-boy.ys168.com (如有侵权请联系删除) 055B-1 视频附带资料:软件-数据-文献下载-持续更新 055B-2 ENMTo ...
- 基于R语言的代理模型(高斯过程、贝叶斯优化、敏感性分析、异方差性等)高级技术应用
基于R语言的代理模型(高斯过程.贝叶斯优化.敏感性分析.异方差性等)高级技术应用 直播时间:10月30日-10月31日.11月6日-7日(4天+1周辅导练习) (上午9:30-12:00 下午14: ...
- r语言解释回归模型的假设_模型假设-解释
r语言解释回归模型的假设 Ever heard of model assumptions? What are they? And why are they important? A model is ...
- r语言 计算模型的rmse_直播丨R语言与作物模型高级应用实战技术应用
随着基于过程的作物生长模型(Process-based Crop Growth Simulation Model)的发展,R语言在作物生长模型和数据分析.挖掘和可视化中发挥着越来越重要的作用.想要成为 ...
- R语言ARIMA集成模型预测时间序列分析
全文链接:http://tecdat.cn/?p=18493 本文我们使用4个时间序列模型对每周的温度序列建模.第一个是通过auto.arima获得的,然后两个是SARIMA模型,最后一个是Buys- ...
- R语言如何计算回归模型参数的置信区间?
R语言如何计算回归模型参数的置信区间? 目录 R语言如何计算回归模型参数的置信区间? R语言是解决什么问题的? R语言如何计算回归模型参数的置信区间? 安利一个R语言的优秀博主及其CSDN专栏: R语 ...
- 基于R语言一元线性回归模型实例及代码
基于R语言一元线性回归模型实例及代码 题目描述 数据特征及可视化 建立模型与初步评价 (自己写lm()代码) 显著性检验 整体显著性检验 数学理论 系数显著性检验 代码实现系统显著性检验 回归诊断 异 ...
最新文章
- 【idea】Springboot整合jpa
- arduino串口绘图_一起打造一款光驱迷你绘图仪
- 7个HTML5移动开发框架,初学HTML5必看
- 版本分支管理标准 - Trunk Based Development 主干开发模型
- 深度阅读之《Concurrency in Go》
- 肾有多好人就有多年轻 男女通用的补肾秘方
- CAN和CANOpen的关系
- python创建文件的编码格式
- pytorch 生成随机数Tensor的方法 torch.rand torch.randn torch.normal torch.linespace
- 【iCore3双核心板】iCore3双核心板使用说明(图文)
- java iecapt.exe_IECapt生成网页快照IECapt.exe下载 CutyCapt
- 计算机操作系统--思维导图
- mouseenter和mouseleave与mouseover和mouseout的区别
- 【计算理论】计算理论总结 ( 上下文无关文法 ) ★★
- js 获取当天23点59分59秒 时间戳
- 怎么查看无线路由器连接的设备连接服务器,路由器怎么看几个人连接
- XXL之整合SpringBoot
- c++编译报错:ld returned 1 exit status
- 【热门收藏】iOS开发人员必看的精品资料(100个)
- 2021J - Circular Billiard Table
热门文章
- 66岁比尔盖茨突然宣布离婚!27年前与下属恋爱修成正果,现在“无法共同成长”,分割8000亿财产...
- Google 放话:要教会我家宝宝开发Android App!
- 数字孪生应用白皮书_工信部发布数字孪生应用白皮书:特斯联入选智慧城市建设标杆案例...
- windows操作系统_windows下用深度系统安装器安装深度操作系统实现双系统分别运行...
- vscode settings.json配置
- 数据结构和算法9——哈希表
- C#基础知识回顾-- 反射(1)
- Environment.CommanLine返回的文件路径使用注意
- 本地项目部署到服务器 启动 报错 数据表不存原因 解决
- 数据库优化之简单理解