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

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

---

1. 课件下载

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 非线性规划

• 最优投资组合

$$\mathop {Max}\limits_\theta U = \theta\mu + (1-\theta)r_f - \frac{1}{2}A{\theta^2}\sigma^2\\ \Rightarrow \theta=\frac{\mu -r_f}{A\sigma^2}$$

mu=0.1; r_f = 0.02; A =10; sigma=0.2
theta= seq(0,1, by=0.05)
U = theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2
data=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()

• 企业最优决策和要素需求
  $$\mathop {Max}\limits_{K,L} \Pi = {\left( {KL} \right)^{0.5}} - 0.1K - 5L \Rightarrow K=500,~~L=10$$

r=0.1; w=5; alpha=0.5; Cost=100; N=12
K=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$Y
dim(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)                            # 最优化结果

• 消费者最优问题

$$\begin{array}{l} \mathop {Max}\limits_{{c_1},{c_2}} {\rm{~}}u\left( {{c_1},{c_2}} \right) = \alpha \ln {c_1} + \left( {1 - \alpha } \right)\ln {c_2}\\ s.t.{\rm{~~~}}G\left( {{c_1},{c_2}} \right) = w - {p_1}{c_1} - {p_2}{c_2} = 0 \end{array}$$

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$U
dim(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/p2
data.frame(c1_check,c2_check)