R语言求解线性规划和非线性规划
非线性规划
1. 非线性规划我们讨论过线性规划,其目标函数和约束条件都是自变量的线性函数。
如果目标函数是非线性函数或至少有一个约束条件是非线性等式(不等式),则这一类数学规划就称为非线性规划。
在科学管理和其他领域中，很多实际问题可以归结为线性规划，但还有另一些问题属于非线性规划。
由于非线性规划含有深刻的背景和丰富的内容，已发展为运筹学的重要分支，并且在最优设计,管理科学,风险管理,系统控制,求解均衡模型,以及数据拟合等领域得到越来越广泛的应用。
非线性规划的研究始于三十年代末,是由W.卡鲁什首次进行的，40年代后期进入系统研究，1951年.库恩和.塔克提出带约束条件非线性规划最优化的判别条件,从而奠定了非线性规划的理论基础,后来在理论研究和实用算法方面都有很大的发展。
非线性规划求解方法可分为无约束问题和带约束问题来讨论，前者实际上就是多元函数的极值问题，是后一问题的基础。
无约束问题的求解方法有最陡下降法、共轭梯度法、变尺度法和鲍威尔直接法等。
关于带约束非线性规划的情况比较复杂，因为在迭代过程中除了要使目标函数下降外，还要考虑近似解的可行性。
总的原则是设法将约束问题化为无约束问题；把非线性问题化为线性问题从而使复杂问题简单化。
求解方法有可行方向法、约束集法、制约函数法、简约梯度法、约束变尺度法、二次规划法等。
虽然这些方法都有较好的效果，但是尚未找到可以用于解决所有非线性规划的统一算法。
非线性规划举例[库存管理问题] 考虑首都名酒专卖商店关于啤酒库存的年管理策略。
假设该商店啤酒的年销售量为A 箱,每箱啤酒的平均库存成本为H 元,每次订货成本都为F 元。
如果补货方式是可以在瞬间完成的,那么为了降低年库存管理费用,商店必须决定每年需要定多少次货,以及每次订货量。
我们以Q 表示每次定货数量,那么年定货次数可以为QA,年订货成本为Q A F ⨯。
由于平均库存量为2Q,所以,年持有成本为2Q H ⨯,年库存成本可以表示为:Q HQ A F Q C ⨯+⨯=2)( 将它表示为数学规划问题:min Q H Q A F Q C ⋅+⋅=2)( ..t s 0≥Q其中Q 为决策变量,因为目标函数是非线性的,约束条件是非负约束,所以这是带约束条件的非线性规划问题。
线性和非线性最优化理论、方法、软件及应用
线性和非线性最优化理论、方法、软件及应用最优化在航空航天、生命科学、水利科学、地球科学、工程技术等自然科学领域和经济金融等社会科学领域有着广泛和重要的应用, 它的研究和发展一直得到广泛的关注. 最优化的研究包含理论、方法和应用.最优化理论主要研究问题解的最优性条件、灵敏度分析、解的存在性和一般复杂性等.而最优化方法研究包括构造新算法、证明解的收敛性、算法的比较和复杂性等.最优化的应用研究则包括算法的实现、算法的程序、软件包及商业化、在实际问题的应用. 这里简介一下线性和非线性最优化理论、方法及应用研究的发展状况.1. 线性最优化线性最优化, 又称线性规划, 是运筹学中应用最广泛的一个分支.这是因为自然科学和社会科学中许多问题都可以近似地化成线性规划问题. 线性规划理论和算法的研究及发展共经历了三个高潮, 每个高潮都引起了社会的极大关注. 线性规划研究的第一高潮是著名的单纯形法的研究. 这一方法是Dantzig在1947年提出的,它以成熟的算法理论和完善的算法及软件统治线性规划达三十多年. 随着60年代发展起来的计算复杂性理论的研究, 单纯形法在七十年代末受到了挑战. 1979年前苏联数学家Khachiyan提出了第一个理论上优于单纯形法的所谓多项式时间算法--椭球法, 曾成为轰动一时的新闻, 并掀起了研究线性规划的第二个高潮. 但遗憾的是广泛的数值试验表明, 椭球算法的计算比单纯形方法差.1984年Karmarkar提出了求解线性规划的另一个多项式时间算法. 这个算法从理论和数值上都优于椭球法,因而引起学术界的极大关注, 并由此掀起了研究线性规划的第三个高潮. 从那以后, 许多学者致力于改进和完善这一算法,得到了许多改进算法.这些算法运用不同的思想方法均获得通过可行区域内部的迭代点列,因此统称为解线性规划问题的内点算法. 目前内点算法正以不可抗拒的趋势将超越和替代单纯形法.线性规划的软件, 特别是由单纯形法所形成的软件比较成熟和完善.这些软件不仅可以解一般线性规划问题, 而且可以解整数线性规划问题、进行灵敏度分析, 同时可以解具有稀疏结构的大规模问题.CPLEX是Bi xby基于单纯形法研制的解线性和整数规划的软件, CPLEX的网址是/. 此外,这个软件也可以用来解凸二次规划问题, 且特别适合解大规模问题. PROC LP是SAS软件公司研制的SAS商业软件中OR模块的一个程序.这个程序是根据两阶段单纯形法研制的,可以用来解线性和整数规划问题并可进行灵敏度分析, 是一个比较完善的程序.用户可以根据需要选择不同的参数来满足不同的要求。
r语言 对指定范围的数据拟合曲线
文章标题：深入探讨R语言对指定范围的数据拟合曲线在数据分析和统计学中，拟合曲线是一种对观测数据进行数学建模的方法。
通过拟合曲线，我们可以揭示数据中的规律和趋势，并进一步进行预测和分析。
R语言作为一种强大的统计分析工具，提供了丰富的函数和包，可以帮助我们对指定范围的数据进行拟合曲线分析。
在本文中，我们将深入探讨R语言对指定范围的数据拟合曲线的方法和应用，以及个人对这一主题的理解和观点。
一、拟合曲线的基本概念拟合曲线是一种通过数学函数来拟合数据的方法，旨在找到最符合数据趋势的函数模型。
在R语言中，拟合曲线通常使用lm()函数进行线性拟合，或者使用非线性拟合的函数，如nls()或者自定义的函数。
通过拟合曲线，我们可以得到数据的趋势线，了解数据的分布规律和变化情况。
二、在R语言中进行数据拟合曲线在R语言中进行数据拟合曲线分析通常需要以下步骤：导入数据并进行数据预处理；选择适当的拟合函数进行数据拟合；对拟合结果进行评估和分析。
在导入数据方面，我们可以使用read.csv()或者read.table()等函数导入数据。
在进行数据拟合前，通常需要对数据进行清洗和处理，去除异常值或缺失值，确保数据的完整性和准确性。
接下来，根据数据的特点和需求，选择适当的拟合函数进行数据拟合。
对于线性拟合，我们可以使用lm()函数进行简单的线性回归分析；对于非线性拟合，我们可以使用nls()函数或者自定义的函数进行拟合。
对拟合结果进行评估和分析，可以通过summary()函数查看拟合参数的具体数值，或者绘制拟合曲线和原始数据的散点图进行可视化比较。
三、在指定范围的数据拟合曲线在实际应用中，我们经常需要对指定范围的数据进行拟合曲线分析。
我们可能只关注数据集中的某一部分数据，或者对数据的某一时间段或空间区域进行拟合。
在R语言中，可以通过子集操作或者筛选条件的方式，对指定范围的数据进行拟合曲线分析。
可以使用subset()函数或者逻辑条件来选择感兴趣的数据，并进行相应的拟合分析。
R语言实验三
实验三数组的运算、求解方程（组）和函数极值、数值积分【实验类型】验证性【实验学时】2 学时【实验目的】1、掌握向量的四则运算和内积运算、矩阵的行列式和逆等相关运算；2、掌握线性和非线性方程（组）的求解方法，函数极值的求解方法；3、了解 R 中数值积分的求解方法。
【实验内容】1、向量与矩阵的常见运算；2、求解线性和非线性方程（组）；3、求函数的极值，计算函数的积分。
【实验方法或步骤】第一部分、课件例题：1.向量的运算x<-c(-1,0,2)y<-c(3,8,2)v<-2*x+y+1vx*yx/yy^xexp(x)sqrt(y)x1<-c(100,200); x2<-1:6; x1+x22.x<-1:5y<-2*1:5x%*%ycrossprod(x,y)x%o%ytcrossprod(x,y)outer(x,y)3.矩阵的运算A<-matrix(1:9,nrow=3,byrow=T);AA+1 #A的每个元素都加上1B<-matrix(1:9,nrow=3); BC<-matrix(c(1,2,2,3,3,4,4,6,8),nrow=3); C D<-2*C+A/B; D #对应元素进行四则运算x<-1:9A+x #矩阵按列与向量相加E<-A%*%B; E #矩阵的乘法y<-1:3A%*%y #矩阵与向量相乘crossprod(A,B) #A的转置乘以Btcrossprod(A,B) #A乘以B的转置4.矩阵的运算A<-matrix(c(1:8,0),nrow=3);At(A) #转置det(A) #求矩阵行列式的值diag(A) #提取对角线上的元素A[lower.tri(A)==T]<-0;A #构造A对应的上三角矩阵qr.A<-qr(A);qr.A #将矩阵A分解成正交阵Q与上三角阵R的乘积，该结果为一列表Q<-qr.Q(qr.A);Q;R<-qr.R(qr.A);R #显示分解后对应的正交阵Q与上三角阵Rdet(Q);det(R);Q%*%R #A=Q*Rqr.X(qr.A) #显示分解前的矩阵5.解线性方程组A<-matrix(c(1:8,0),nrow=3,byrow=TRUE)b<-c(1,1,1)x<-solve(A,b); x #解线性方程组Ax=bB<-solve(A); B #求矩阵A的逆矩阵BA%*%B #结果为单位阵6.非线性方程求根f<-function(x) x^3-x-1 #建立函数uniroot(f,c(1,2)) #输出列表中f.root为近似解处的函数值，iter为迭代次数，estim.prec为精度的估计值uniroot(f,lower=1,upper=2) #与上述结果相同polyroot(c(-1,-1,0,1)) #专门用来求多项式的根，其中c(-1,-1,0,1)表示对应多项式从零次幂项到高次幂项的系数7.求解非线性方程组（1）自编函数： (Newtons.R)Newtons<-function (funs, x, ep=1e-5, it_max=100){index<-0; k<-1while (k<=it_max){ #it_max 表示最大迭代次数x1 <- x; obj <- funs(x);x <- x - solve(obj$J, obj$f); #Newton 法的迭代公式norm <- sqrt((x-x1) %*% (x-x1))if (norm<ep){ index<-1; break #index=1 表示求解成功}; k<-k+1 }obj <- funs(x);list(root=x, it=k, index=index, FunVal= obj$f)} # 输出列表（2）调用求解非线性方程组的自编函数funs<-function(x){ f<-c(x[1]^2+x[2]^2-5, (x[1]+1)*x[2]-(3*x[1]+1)) # 定义函数组J<-matrix(c(2*x[1], 2*x[2], x[2]-3, x[1]+1), nrow=2,byrow=T) # 函数组的 Jacobi 矩阵list(f=f, J=J)} # 返回值为列表 : 函数值 f 和 Jacobi 矩阵 Jsource("F:/wenjian_daima/Newtons.R") # 调用求解非线性方程组的自编函数Newtons(funs, x=c(0,1))8.一元函数极值f<-function(x) x^3-2*x-5 # 定义函数optimize(f,lower=0,upper=2) # 返回值 : 极小值点和目标函数f<-function(x,a) (x-a)^2 # 定义含有参数的函数optimize(f,interval=c(0,1),a=1/3) # 在函数中输入附加参数9.多元函数极值（1）obj <-function (x){ # 定义函数F<-c(10*(x[2]-x[1]^2),1-x[1]) # 视为向量sum (F^2) } # 向量对应分量平方后求和nlm(obj,c(-1.2,1))（2）fn<-function(x){ # 定义目标函数F<-c(10*(x[2]-x[1]^2), 1-x[1])t(F)%*%F } # 向量的内积gr <- function(x){ # 定义梯度函数F<-c(10*(x[2]-x[1]^2), 1-x[1])J<-matrix(c(-20*x[1],10,-1,0),2,2,byrow=T) #Jacobi 矩阵2*t(J)%*%F } # 梯度optim(c(-1.2,1), fn, gr, method="BFGS")最优点 (par) 、最优函数值 (value)10.梯形求积分公式（1）求积分程序： (trape.R)trape<-function(fun, a, b, tol=1e-6){ # 精度为 10 -6N <- 1; h <- b-a ; T <- h/2 * (fun(a) + fun(b)) # 梯形面积 repeat{h <- h/2; x<-a+(2*1:N-1)*h; I <-T/2 + h*sum(fun(x)) if(abs(I-T) < tol) break; N <- 2 * N; T = I }; I}（2）source("F:/wenjian_daima/trape.R") # 调用函数f<-function(x) exp(-x^2)trape(f,-1,1)（3）常用求积分函数f<-function(x)exp(-x^2) # 定义函数integrate(f,0,1)integrate(f,0,10)integrate(f,0,100)integrate(f,0,10000) # 当积分上限很大时，结果出现问题integrate(f,0,Inf) # 积分上限为无穷大ft<-function(t) exp(-(t/(1-t))^2)/(1-t)^2 # 对上述积分的被积函数 e 2 作变量代换 t=x/(1+x) 后的函数integrate(ft,0,1) # 与上述计算结果相同，且精度较高第二部分、教材例题：1.随机抽样（1）等可能的不放回的随机抽样:> sample(x, n) 其中x为要抽取的向量, n为样本容量（2）等可能的有放回的随机抽样:> sample(x, n, replace=TRUE)其中选项replace=TRUE表示有放回的, 此选项省略或replace=FALSE表示抽样是不放回的sample(c("H", "T"), 10, replace=T)sample(1:6, 10, replace=T)（3）不等可能的随机抽样:> sample(x, n, replace=TRUE, prob=y)其中选项prob=y用于指定x中元素出现的概率, 向量y与x等长度sample(c("成功", "失败"), 10, replace=T, prob=c(0.9,0.1))sample(c(1,0), 10, replace=T, prob=c(0.9,0.1))2.排列组合与概率的计算1/prod(52:49)1/choose(52,4)3.概率分布qnorm(0.025) #显著性水平为5%的正态分布的双侧临界值qnorm(0.975)1 - pchisq(3.84, 1) #计算假设检验的p值2*pt(-2.43, df = 13) #容量为14的双边t检验的p值4.limite.central( )的定义limite.central <- function (r=runif, distpar=c(0,1), m=.5,s=1/sqrt(12),n=c(1,3,10,30), N=1000) {for (i in n) {if (length(distpar)==2){x <- matrix(r(i*N, distpar[1],distpar[2]),nc=i)}else {x <- matrix(r(i*N, distpar), nc=i)}x <- (apply(x, 1, sum) - i*m )/(sqrt(i)*s)hist(x,col="light blue",probability=T,main=paste("n=",i), ylim=c(0,max(.4, density(x)$y)))lines(density(x), col="red", lwd=3)curve(dnorm(x), col="blue", lwd=3, lty=3, add=T)if( N>100 ) {rug(sample(x,100))}else {rug(x)}}}5.直方图x=runif(100,min=0,max=1)hist(x)6.二项分布B（10，0.1）op <- par(mfrow=c(2,2))limite.central(rbinom,distpar=c(10,0.1),m=1,s=0.9)par(op)7.泊松分布: pios（1）op <- par(mfrow=c(2,2))limite.central(rpois, distpar=1, m=1, s=1, n=c(3, 10, 30 ,50)) par(op)8.均匀分布：unif(0,1)op <- par(mfrow=c(2,2))limite.central( )par(op)9.指数分布：exp(1)op <- par(mfrow=c(2,2))limite.central(rexp, distpar=1, m=1, s=1)par(op)10.混合正态分布的渐近正态性mixn <- function (n, a=-1, b=1){rnorm(n, sample(c(a,b),n,replace=T))}limite.central(r=mixn, distpar=c(-3,3),m=0, s=sqrt(10), n=c(1,2,3,10)) par(op)11.混合正态分布的渐近正态性op <- par(mfrow=c(2,2))mixn <- function (n, a=-1, b=1){rnorm(n, sample(c(a,b),n,replace=T))}limite.central(r=mixn, distpar=c(-3,3),m=0,s=sqrt(10),n=c(1,2,3,10)) par(op)第三部分、课后习题：3.1a=sample(1:100,5)asum(a)3.2（1）抽到10、J、Q、K、A的事件记为A，概率为P(A)=1(5220)其中在R中计算得：> 1/choose(52,20)[1] 7.936846e-15（2）抽到的是同花顺P(B)=(41)(91) (525)在R中计算得：> (choose(4,1)*choose(9,1))/choose(52,5) [1] 1.385e-053.3#(1)x<-rnorm(1000,mean=100,sd=100)hist(x)#(2)y<-sample(x,500)hist(y)#(3)mean(x)mean(y)var(x)var(y)3.4x<-rnorm(1000,mean=0,sd=1) y=cumsum(x)plot(y,type = "l")plot(y,type = "p")3.5x<-rnorm(100,mean=0,sd=1) qnorm(.025)qnorm(.975)t.test(x)由R结果知：理论值为[-1.96,1.96]，实际值为：[-0.07929,0.33001]3.6op <- par(mfrow=c(2,2))limite.central(rbeta, distpar=c(0.5 ,0.5),n=c(30,200,500,1000))par(op)3.7N=seq(-4,4,length=1000)f<-function(x){dnorm(x)/sum(dnorm(x))}n=f(N)result=sample(n,replace=T,size = 1000)standdata=rnorm(1000)op<-par(mfrow=c(1,2)) #1行2列数组按列（mfcol）或行（mfrow）各自绘图hist(result,probability = T)lines(density(result),col="red",lwd=3)hist(standdata,probability = T)lines(density(standdata),col="red",lwd=3) par(op)。
实验二利用Lingo求解整数规划及非线性规划问题
例 3 用Lingo软件求解非线性规划问题
min z x1 12 x2 22
x2 x1 1,
x1
x2
2,
x1
0,
x2
0.
Lingo 程序: min= x1-1 ^2+ x2-2 ^2;
x2-x1=1;
x1+x2<=2;
注意: Lingo 默认变量的取值从0到正无穷大, 变量定界函数可以改变默认状态. @free x : 取消对变量x的限制 即x可取任意实数值
例 4 求函数 zx22y22 的最小值.
例 4 求函数 zx22y22 的最小值.
解: 编写Lingo 程序如下:
min= x+2 ^2+ y-2 ^2; @free x ; 求得结果: x=-2, y=2
二、Lingo 循环编程语句
1 集合的定义 包括如下参数: 1 集合的名称.
sets: endsets
44
minZ
aijxij
i1 j1
4
xij
1
j 1,2,3,4
s.t.
i 1 4
xij
1
i 1,2,3,4
j1
xij 0或1 i, j 1,2,3,4
LINGO程序如下：
MODEL: SETS: person/A,B,C,D/; task/1..4/; assign person,task :a,x; ENDSETS DATA: a=1100,800,1000,700,
77
63
67
丁
55
76
62
62
甲, 乙, 丙, 丁 四名队员各自游什么姿势 , 才最有可能取得好成绩
线性规划与非线性规划
一、线性规划
模型 设在甲车床上加工工件1、2、3的数量分别为 在乙车床上加工工件1、2、3的数量分别为
问题二：某厂每日8小时的产量不低于1800件.为了进行质量控制，计划聘请两种不同水平的检验员.一级检验员的标准为：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15件/小时，正确率95%，计时工资3元/小时.检验员每错检一次，工厂要损失2元.为使总检验费用最省，该工厂应聘一级、二级检验员各几名？
01
模型 设需要一级、二级检验员的人数分别为 人， 应付检验员工资为 因检验员错检而造成的损失为
02
注；当前MATLAB只支持第一种形式。
或矩阵形式 其中 是决策变量， 是约束矩阵， ，
二、非线性规划
1、二次规划 标准形式： MATLAB调用格式： (1) x=quadprog(H,C,A1,b1); (2)x=quadprog(H,C,A1,b1,A2,b2,v1,v2); (3)[x,fval,exitflag,output]= quadprog(H,C,A1,b1, A2,b2 ,v1,v2,x0,options);
见MATLAB程序（xianxingguihua4）
例4：问题二的解答 改写为
结果： 即只需聘用9个一级检验员。 注：本问题应还有一个约束条件：x1、x2取整数，故它属于一个整数线性规划问题，这里当成一个线性规划求解，求得最优解刚好是整数x1=9，x2=0，故它就是该整数规划的最优解.若用线性规划解法求得的最优解不是整数，将其取整后不一定是相应整数规划的最优解，这样的整数规划应用专门的方法求解.
2、状态窗口(LINDO Solver Status)
当前状态：已达最优解 迭代次数：18次 约束不满足的“量”(不是“约束个数”)：0 当前的目标值：94 最好的整数解：94 整数规划的界：93.5 分枝数：1 所用时间：0.00秒（太快了，还不到0.005秒） 刷新本界面的间隔:1(秒)
第5章 非线性规划
(水力约束) (水力摩阻系数约束)
KD GC
L (热力约束)
(粘温关系约束)
(工艺要求约束) (管道强度约束)
在目标函数中，f1(TR)、f2(Pd)一般为非线性函数，约束条 件中亦存在不少非线性函数，显然是一个NLP问题。
非线性规划的基本概念和定理
例3：最小二乘问题：该问题大量存在于工业生产和科学 实验的数据处理中。例如原油的粘度可以表示为：
凹函数的几何意义：
对 于 一 元 函 数 f(x) ， 若
函数曲线上任意两点之 间的连线永远不在曲线
的 上 方 ， 则 f(x) 为 凹 函
数(参见右图) 。
非线性规划的基本概念和定理 f(X)
f [X 1 (1 ) X 2 ]
对于二元函数 f(x1,x2)， 若函数曲面上任意两点 之间的连线永远不在曲 面的上方，则f(x1,x2)为 凹函数(参见右图)。
1、一元函数：
①必要条件：f(x)在x*处取得极值的必要条件是f'(x*)＝0；
②充分条件：若f"(x*)＜0，则x*为极大点； 若f"(x*)＞0，则x*为极小点。 2、多元函数： ①必要条件： f(X)在D域内存在极值点X*的必要条件为 * f ( X ) 0 (即f(X)在X*处的所有一阶偏导数等于0)。
非线性规划的基本概念和定理
根据定义，线性函数既是凸函数，又是凹函数。 凸函数的几何意义： 对 于 一 元 函 在曲线的下方， 则 f(x) 为 凸 函 数 ( 参 见 右
图) 。
非线性规划的基本概念和定理 f(X)
f ( X 1 ) (1 ) f ( X 2 )
§5.1 非线性规划的基本概念和定理
一、什么是非线性规划?
非线性规划
非线性规划（nonlinear programming）1.非线性规划概念非线性规划是具有非线性约束条件或目标函数的数学规划，是运筹学的一个重要分支。
非线性规划研究一个n元实函数在一组等式或不等式的约束条件下的极值问题，且目标函数和约束条件至少有一个是未知量的非线性函数。
目标函数和约束条件都是线性函数的情形则属于线性规划。
2.非线性规划发展史公元前500年古希腊在讨论建筑美学中就已发现了长方形长与宽的最佳比例为0.618,称为黄金分割比。
其倒数至今在优选法中仍得到广泛应用。
在微积分出现以前，已有许多学者开始研究用数学方法解决最优化问题。
例如阿基米德证明：给定周长，圆所包围的面积为最大。
这就是欧洲古代城堡几乎都建成圆形的原因。
但是最优化方法真正形成为科学方法则在17世纪以后。
17世纪，I.牛顿和G.W.莱布尼茨在他们所创建的微积分中，提出求解具有多个自变量的实值函数的最大值和最小值的方法。
以后又进一步讨论具有未知函数的函数极值，从而形成变分法。
这一时期的最优化方法可以称为古典最优化方法。
最优化方法不同类型的最优化问题可以有不同的最优化方法，即使同一类型的问题也可有多种最优化方法。
反之,某些最优化方法可适用于不同类型的模型。
最优化问题的求解方法一般可以分成解析法、直接法、数值计算法和其他方法。
(1)解析法：这种方法只适用于目标函数和约束条件有明显的解析表达式的情况。
求解方法是：先求出最优的必要条件，得到一组方程或不等式，再求解这组方程或不等式，一般是用求导数的方法或变分法求出必要条件，通过必要条件将问题简化，因此也称间接法。
(2)直接法：当目标函数较为复杂或者不能用变量显函数描述时，无法用解析法求必要条件。
此时可采用直接搜索的方法经过若干次迭代搜索到最优点。
这种方法常常根据经验或通过试验得到所需结果。
对于一维搜索(单变量极值问题)，主要用消去法或多项式插值法；对于多维搜索问题(多变量极值问题)主要应用爬山法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第七章线性规划与非线性规划
例1m a x z=10x
1+5x2
s.t.5x1+2x2<=8
3x1+4x2=9
x1+x2>=1
x1,x2>=0
首先可化为标准形式：min - z = -10x1 -5x2
s.t. 5x1+2x1<=8
-x1-x2<=-1
3x1+4x2=9
x1,x2>=0
library(Rglpk)
obj<-c(-10,-5)
mat<-matrix(c(5,2,-1,-1,3,4),3,2,T)
dir<-c("<=","<=","==")
rhs<-c(8,-1,9)
Rglpk_solve_LP(obj,mat,dir,rhs)
#直接求解
library(Rglpk)
obj<-c(10,5)
mat<-matrix(c(5,2,1,1,3,4),3,2,T)
dir<-c("<=",">=","==")
rhs<-c(8,1,9)
Rglpk_solve_LP(obj,mat,dir,rhs,max=T)
非线性规划求解（Rdonlp2）
例2 有如下的条件约束最优化问题：
22min(sin cos )
1001001001002133
2sin cos 3z x y y x x y x y x y xy x y =+-<<⎧⎪-<<⎪⎪≤+⎨≤-≤⎪⎪=⎪≤⎩
library (Rdonlp2) p = c(10,10) #迭代初始值
#对求解问题进行描述
22min(sin cos )z x y y x =+
fn = function (x){
x[1]^2*sin(x[2])+x[2]^2*cos(x[1])
}
#对x,y 值域描述
100100100100
x y -<<-<< ## par.l 和par.u 分别为约束的左边和右边
par.l = c(-100,-100); par.u = c(100,100) ## 目标值域
#对线性约束进行描述
2133
x y x y ≤+≤-≤ A = matrix(c(1,1,3,-1),2,byrow=TRUE ) ##线性约束系数
lin.l = c(2,1); lin.u = c(+Inf ,3) ## 分别为约束的左边和右边
#对非线性约束进行描述
2sin cos 3
xy x y =≤ nlcon1 = function (x){
x[1]*x[2] ##公式 x*y
}
nlcon2 = function (x){
sin(x[1])*cos(x[2]) ##公式 sin(x)*cos(y)
}
## 两个非线性约束的左右边
## x*y=2 等价于 2<=x*y<=2
nlin.l = c(2,-Inf ) ; nlin.u = c(2,0.6) #将参数输入donlp2函数中进行求解
## 输入参数第一行： x,y 值域及目标函数
## 输入参数第二行： 线性约束条件
## 输入参数第三，四行： 非线性约束条件
ret = donlp2(p, fn, par.u=par.u, par.l=par.l,
A, lin.l=lin.l,lin.u=lin.u,
nlin=list(nlcon1,nlcon2),
nlin.u=nlin.u, nlin.l=nlin.l)
## 输出结果
ret$par
ret$par
例3 解下列二次规划
[][]111222111()26122x x f x x x x x -⎡⎤⎡⎤⎡⎤=-⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦
二次规划的优化问题，这是一种特殊形式的非线性约束优化问题。
二次规划在许多领域都有运用，比如投资组合优化、求解支持向量机(SVM)分类问题等。
想要用quadprog 包求解二次规划，我们需要同时转化我们的目标函数和约束条件为矩阵形式。
quadprog 包默认是求解最小化问题，目标函数二次，约束一次。
所以，我们的约束条件默认的形式也就是AX>=bvec 。
通常我们需要把一些原来是求极大值的问题或者<=约束通过乘以负号来转化 library(quadprog)
Dmat <-matrix(c(1,-1,1,2),2,2,T)
Dmat
dvec <- c(2,6)
A<- matrix(-c(1,1,-1,2,2,1),3,2,T)
A
bvec <- c(-2,-2,-3)
Amat <- t(A)
sol <- solve.QP(Dmat, dvec, Amat, bvec)
sol
参数Dmat 表示海赛矩阵
参数dvet 表示一阶向量，和Dmat 的维数要相对应。
参数Amat 表示约束矩阵，默认的约束都是>=。
参数bvet 表示右边值，由向量，和Amat 的维数要相对应。
参数 meq 表示从哪一行开始Amat 矩阵中的约束是需要被当作等式约束的。
例4 假设以决策变量x1、x2、x3分别表示甲、乙、丙、丁4种肥料的用量，得线性规划模型
1234124134141234min 0.040.150.10.125..0.030.30.15320.050.20.1240.140.0742,,,0z x x x x s t x x x x x x x x x x x x =+++⎧⎪++≥⎪⎪++=⎨⎪+≤⎪⎪≥⎩
library(Rglpk)
obj<-c(0.04,0.15,0.1,0.125)
mat<-matrix(c(0.03,0.3,0,0.15,0.05,0,0.2,0.1,0.14,0,0,0.07),3,4,T)
mat
dir<-c(">=","==","<=")
rhs<-c(32,24,42)
Rglpk_solve_LP(obj,mat,dir,rhs)。