[关闭]
@fanxy 2016-07-08T07:57:30.000000Z 字数 3162 阅读 2321

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

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


1. 课件下载

Ch04.pdf868.1kB

2. 例题程序附录

  1. # 准备工作
  2. setwd("D:\\...")
  3. install.packages("Rdonlp2", repos="http://R-Forge.R-project.org")
  4. install.packages("lpSolve")
  5. install.packages(c("ggplot2","plotly"))
  6. library(lpSolve)
  7. library(Rdonlp2)
  8. library(ggplot2)
  9. library(plotly)

2.1 线性规划

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

2.2 非线性规划

  1. mu=0.1; r_f = 0.02; A =10; sigma=0.2
  2. theta= seq(0,1, by=0.05)
  3. U = theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2
  4. data=data.frame(theta, U)
  5. ggplot(data) +
  6. geom_line(aes(x=theta, y=U)) +
  7. geom_vline(xintercept=(mu-r_f)/A/sigma^2, color="red") +
  8. labs(x="风险资产占比", y="效用") +
  9. theme_bw()
  1. r=0.1; w=5; alpha=0.5; Cost=100; N=12
  2. K=seq(0.01,Cost/r, length.out=N)
  3. L=seq(0.01,Cost/w, length.out=N)
  4. # 用ggplot2包做二维图
  5. firmdata <- transform(expand.grid(K=K,L=L), Y=(K*L)^alpha)
  6. ggplot(firmdata) +
  7. stat_contour(aes(x=K, y=L, z=Y, colour=Y)) +
  8. geom_abline(slope=-r/w, intercept=Cost/w, color="red") +
  9. labs(xlab="K", ylab="L") +
  10. theme_bw()
  11. # 用plotly包作三维图
  12. Y=firmdata$Y
  13. dim(Y)=c(N,N)
  14. plot_ly(z=Y, type="contour") %>% # 两维产出等高线
  15. layout(scene = list(
  16. xaxis = list(title = "K"),
  17. yaxis = list(title = "L"),
  18. zaxis = list(title = "Y")))
  19. plot_ly(z=Y, type="surface") %>% # 三维产出曲面
  20. layout(scene = list(
  21. xaxis = list(title = "K"),
  22. yaxis = list(title = "L"),
  23. zaxis = list(title = "Y")))
  24. # 求解
  25. par.l= c(0,0); par.u = c(100/r,100/w) # 参数的下限和上限
  26. fn=function(x){ # 目标函数
  27. f=(x[1]*x[2])^alpha-r*x[1]-w*x[2]
  28. return(-f) # 默认为最小化,加负号改为最大化
  29. }
  30. A=matrix(c(r,w),nrow=1) # 线性约束
  31. lin.l=c(Cost); lin.u=c(Cost) # 线性约束上限和下限
  32. p0=c(1,1) # 参数(K和L)初始值
  33. ret=donlp2(p0,fn,par.l=par.l,par.u=par.u,A,lin.l=lin.l,lin.u=lin.u) # 解最优化问题
  34. K=ret$par[1]; L=ret$par[2]
  35. data.frame(K,L) # 最优化结果

  1. alpha=0.4; p1=10; p2=5; inc=100; N=12
  2. # 作图
  3. c1=seq(0.01, inc/p1, length.out=N)
  4. c2=seq(0.01, inc/p2, length.out=N)
  5. data <- transform(expand.grid(c1=c1,c2=c2),
  6. U=alpha*log(c1)+(1-alpha)*log(c2))
  7. U=data$U
  8. dim(U)=c(N,N)
  9. plot_ly(z=U, type="contour") %>%
  10. layout(scene = list(
  11. xaxis = list(title = "C1"),
  12. yaxis = list(title = "C2"),
  13. zaxis = list(title = "U")))
  14. # 求解
  15. par.l= c(0,0); par.u = c(inc/p1,inc/p2)
  16. u=function(x){
  17. u=alpha*log(x[1])+(1-alpha)*log(x[2])
  18. return(-u)
  19. }
  20. A=matrix(c(p1,p2),1)
  21. lin.l=c(inc); lin.u=c(inc)
  22. p0=c(1,1)
  23. ret=donlp2(p0,fn=u,par.l=par.l,par.u=par.u,A,lin.l=lin.l,lin.u=lin.u)
  24. c1=ret$par[1]; c2=ret$par[2]
  25. data.frame(c1,c2)
  26. # 结果检查
  27. c1_check=alpha*inc/p1
  28. c2_check=(1-alpha)*inc/p2
  29. data.frame(c1_check,c2_check)
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注