@fanxy 2016-07-08T07:57:30.000000Z 字数 3162 阅读 2350

# 第四讲 最优化问题基础与应用

樊潇彦 复旦大学经济学院 经济数学

Ch04.pdf868.1kB

# 2. 例题程序附录

# 准备工作setwd("D:\\...") install.packages("Rdonlp2", repos="http://R-Forge.R-project.org")install.packages("lpSolve")install.packages(c("ggplot2","plotly"))library(lpSolve)library(Rdonlp2)library(ggplot2)library(plotly)

## 2.1 线性规划

# 生产计划production = lp(objective.in=c(24, 18),            const.mat=matrix(c(1, 1, 2, 3, 3, 2, 1, 0, 0,1), nrow=5, byrow=T),            const.rhs=c(150, 240, 300,0,0),            const.dir=c("<=", "<=", "<=",">=",">="),             direction="max")production$solution # 规划结果production$objval                          # 目标值# 运输问题transportation = lp(objective.in=c(4,5,6,5,2,4),            const.mat=rbind(matrix(c(1, 1, 1, 0, 0, 0,                                      0, 0, 0, 1, 1, 1,                                     1, 0, 0, 1, 0, 0,                                     0, 1, 0, 0, 1, 0,                                     0, 0, 1, 0, 0, 1), nrow=5,byrow=T),                            diag(6)),            const.rhs=c(12, 8, 8,6,6,rep(0,6)),            const.dir=c("<=", "<=",     ">=",">=",">=",                         ">=",">=",">=", ">=",">=",">="),             direction="min")transportation$solution # 规划结果transportation$objval                      # 目标值

## 2.2 非线性规划

• 最优投资组合

mu=0.1; r_f = 0.02; A =10; sigma=0.2theta= seq(0,1, by=0.05)U = theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2data=data.frame(theta, U)ggplot(data) +  geom_line(aes(x=theta, y=U)) +  geom_vline(xintercept=(mu-r_f)/A/sigma^2, color="red") +  labs(x="风险资产占比", y="效用") +  theme_bw()
• 企业最优决策和要素需求
r=0.1; w=5; alpha=0.5; Cost=100; N=12K=seq(0.01,Cost/r, length.out=N)L=seq(0.01,Cost/w, length.out=N)# 用ggplot2包做二维图firmdata <- transform(expand.grid(K=K,L=L), Y=(K*L)^alpha)ggplot(firmdata) + stat_contour(aes(x=K, y=L, z=Y, colour=Y)) +  geom_abline(slope=-r/w, intercept=Cost/w, color="red") +  labs(xlab="K", ylab="L") +  theme_bw()# 用plotly包作三维图Y=firmdata$Ydim(Y)=c(N,N) plot_ly(z=Y, type="contour") %>% # 两维产出等高线 layout(scene = list( xaxis = list(title = "K"),  yaxis = list(title = "L"),  zaxis = list(title = "Y")))plot_ly(z=Y, type="surface") %>% # 三维产出曲面 layout(scene = list( xaxis = list(title = "K"),  yaxis = list(title = "L"),  zaxis = list(title = "Y")))# 求解par.l= c(0,0); par.u = c(100/r,100/w) # 参数的下限和上限fn=function(x){ # 目标函数 f=(x[1]*x[2])^alpha-r*x[1]-w*x[2] return(-f) # 默认为最小化，加负号改为最大化}A=matrix(c(r,w),nrow=1) # 线性约束lin.l=c(Cost); lin.u=c(Cost) # 线性约束上限和下限p0=c(1,1) # 参数（K和L）初始值ret=donlp2(p0,fn,par.l=par.l,par.u=par.u,A,lin.l=lin.l,lin.u=lin.u) # 解最优化问题K=ret$par[1]; L=ret$par[2]data.frame(K,L) # 最优化结果 • 消费者最优问题 alpha=0.4; p1=10; p2=5; inc=100; N=12# 作图c1=seq(0.01, inc/p1, length.out=N)c2=seq(0.01, inc/p2, length.out=N)data <- transform(expand.grid(c1=c1,c2=c2),  U=alpha*log(c1)+(1-alpha)*log(c2))U=data$Udim(U)=c(N,N) plot_ly(z=U, type="contour") %>%  layout(scene = list(           xaxis = list(title = "C1"),            yaxis = list(title = "C2"),            zaxis = list(title = "U")))# 求解par.l= c(0,0); par.u = c(inc/p1,inc/p2)    u=function(x){                               u=alpha*log(x[1])+(1-alpha)*log(x[2])  return(-u)                                 }A=matrix(c(p1,p2),1)                       lin.l=c(inc); lin.u=c(inc)   p0=c(1,1)                                   ret=donlp2(p0,fn=u,par.l=par.l,par.u=par.u,A,lin.l=lin.l,lin.u=lin.u) c1=ret$par[1]; c2=ret$par[2]data.frame(c1,c2) # 结果检查 c1_check=alpha*inc/p1                      c2_check=(1-alpha)*inc/p2data.frame(c1_check,c2_check) 

• 私有
• 公开
• 删除