MIS 381N SCO HW7
Yuwen WANG
April 25, 2016
LJIWY
Problem 1
1.1 Linear programming
# create lp model
path.1 = make.lp(0,10*10)
path.2 = make.lp(0,10*10)
# Set up the cost matrix
lmatrix = matrix(0, nrow = 10, ncol = 10)
lmatrix[1,3] = lmatrix[3,6] = lmatrix[6,10] = lmatrix[9,10] = 3
lmatrix[1,4] = lmatrix[3,7] = lmatrix[2,4] = lmatrix[4,8] = lmatrix[5,7] = lmatrix[5,8] = lmatrix[8,10] = 2
lmatrix[2,3] = lmatrix[2,5] = lmatrix[4,6] = lmatrix[7,10] = 4
lmatrix[4,9] = 5
# Set objective function
set.objfn(path.1,as.vector(t(lmatrix)))
set.objfn(path.2,as.vector(t(lmatrix)))
# Set objective direction
lp.control(path.1,sense='min')## $anti.degen
## [1] "fixedvars" "stalling"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] -1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "minimize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"lp.control(path.2,sense='min')## $anti.degen
## [1] "fixedvars" "stalling"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] -1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "minimize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"# Add constraints
## Force coef for non-exist paths to be 0
rhs1 = c(1,rep(0,8),-1) # Enter from node 1
rhs2 = c(0,1,rep(0,7),-1) # Enter from node 2
for (n in 1:10){
coef=c(lmatrix[n,1:10]/lmatrix[n,1:10],-lmatrix[1:10,n]/lmatrix[1:10,n])
ind=c((n-1)*10+c(1:10),(c(1:10)-1)*10+n)
nz=is.finite(coef)
add.constraint(path.1,coef[nz], "=",rhs1[n],ind[nz])
}
for (n in 1:10){
coef=c(lmatrix[n,1:10]/lmatrix[n,1:10],-lmatrix[1:10,n]/lmatrix[1:10,n])
ind=c((n-1)*10+c(1:10),(c(1:10)-1)*10+n)
nz=is.finite(coef)
add.constraint(path.2,coef[nz], "=",rhs2[n],ind[nz])
}
## Binary choice variables
set.type(path.1,c(1:100),'binary')
set.type(path.2,c(1:100),'binary')
# Solve lp
solve(path.1)## [1] 0solve(path.2)## [1] 0print(get.objective(path.1))## [1] 6print(matrix(get.variables(path.1),10,10,byrow = T))## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 0 0 1 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 1 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0 1
## [9,] 0 0 0 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0print(get.objective(path.2))## [1] 6print(matrix(get.variables(path.2),10,10,byrow = T))## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 1 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 1 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0 1
## [9,] 0 0 0 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0The shortest path from node 1 to node 10 is 1–> 4–> 8 –>10. The total length is 6.
The shortest path from node 1 to node 10 is 2–> 4–> 8 –>10. The total length is 6.
1.2 Dynamic Programming
Alternatively, we can walk backwards in time. The network can be divided into 4 stages, where nodes 1~2 belongs to Stage 1, 3~5 belongs to Stage 2, 6~9 belongs to Stage 3, and 10 belongs to Stage 4.
\((1)\ State:\ Stage\ i,\ i = 1,2,3,4\)
\((2)\ Action:\ all\ paths\ available\)
\((3)\ B.E.:\)
\(L(4)=0\)
\(L(i)=min\{p_{i(i+1)}\}+L(i+1)\)
where \(L(i)\) is the shortest path from Stage i to Stage 4, \(p_{ij}\) is the length of any available path from Stage i to Stage j.
Problem 2
2.1 Linear Programming
2.1.1 LP Formulation
Let \(N_1i\) stand for not sell at year \(i\), for \(i=1,2,...,6\).
Let \(S_i\) stand for sell at year \(i\), for \(i=1,2,...,6\).
\(min\ 20000-\sum S_iRv_i+\sum N_iOC_i\)
\(s.t.\)
\((1)\ N_i+S_i \leq 1,\ i=1,2,...,6\)
\((2)\ \sum S_i=1,\ i=1,2,...,6\)
\((3)\ S_i + \sum_{t=i+1}^6 N_t \leq 1,\ i=1,2,...,5\)
2.1.2 LP Solution
cost = make.lp(0,12)
# Set objective function
OC = c(6,10,16,24,32,44) * 100
RV = -c(14,12,8,6,4,2) * 1000
obj.coef = c(OC,RV)
set.objfn(cost,obj.coef)
# set objective direction
lp.control(cost,sense='min')## $anti.degen
## [1] "fixedvars" "stalling"
##
## $basis.crash
## [1] "none"
##
## $bb.depthlimit
## [1] -50
##
## $bb.floorfirst
## [1] "automatic"
##
## $bb.rule
## [1] "pseudononint" "greedy" "dynamic" "rcostfixing"
##
## $break.at.first
## [1] FALSE
##
## $break.at.value
## [1] -1e+30
##
## $epsilon
## epsb epsd epsel epsint epsperturb epspivot
## 1e-10 1e-09 1e-12 1e-07 1e-05 2e-07
##
## $improve
## [1] "dualfeas" "thetagap"
##
## $infinite
## [1] 1e+30
##
## $maxpivot
## [1] 250
##
## $mip.gap
## absolute relative
## 1e-11 1e-11
##
## $negrange
## [1] -1e+06
##
## $obj.in.basis
## [1] TRUE
##
## $pivoting
## [1] "devex" "adaptive"
##
## $presolve
## [1] "none"
##
## $scalelimit
## [1] 5
##
## $scaling
## [1] "geometric" "equilibrate" "integers"
##
## $sense
## [1] "minimize"
##
## $simplextype
## [1] "dual" "primal"
##
## $timeout
## [1] 0
##
## $verbose
## [1] "neutral"# Add constraints
## (1) Ni + Si <= 1
for (i in 1:6){
add.constraint(cost, xt = rep(1,2),
type = '<=', rhs = 1,
indices = c(i,i+6))
}
## (2) sum(Si) = 1
add.constraint(cost, xt = rep(1,6),
type = '=', rhs = 1,
indices = c(7:12))
## (3) Si + N(i+1) + ... <= 1
for (i in 1:5){
add.constraint(cost, xt = rep(1,8-i),
type = '<=', rhs = 1,
indices = c(i:6,i+6))
}
solve(cost)## [1] 0print(get.objective(cost)+20000)## [1] 6000print(get.variables(cost))## [1] 0 0 0 0 0 0 1 0 0 0 0 0The solution with minimal net cost is selling the car at year 1.
2.2 Dynamic Programming
2.2.1 DP Formulation
Alternatively, we can walk backwards in time.
\((1)\ State:\ year\ t,\ t=0,1,2,...,6\)
\((2)\ Action:\ sell(S),\ not\ sell(N)\)
\((3)\ B.E.:\)
\(\ v(6)=min(c(S,6),c(N,6))+20000\)
\(\ v(i)=min(c(S,i),c(N,i)+v(i+1))+20000,\ i=1,2,...,5\)
where \(c(A,i)\) is the cost of action(S or N) at state i, \(v(i)\) is the total cost since state i.
2.2.2 DP Solution
cN = c(6,10,16,24,32,44)*100
cS = -c(14,12,8,6,4,2)*1000
# place-holders for value and choice for each state
v = rep(0,6)
choice = rep(0,6)
for (i in 6:1){
if (i==6){
v[i] = min(cN[i], cS[i]) + 20000
if (cN[i]>cS[i]){
choice[i] = 'S'
}
} else {
v[i] = min(cN[i] + v[i+1], cS[i]) + 20000
if (cN[i] + v[i+1] > cS[i]){
choice[i] = 'S'
choice[(i+1):6] = 'N'
}
}
# print(choice[i])
# print(v[i])
}
print(choice)## [1] "S" "N" "N" "N" "N" "N"
Problem 3
DP Formulation
\((1)\ State:\ (C,t),\ C=\{B,C,I\},\ t=1,2,3,4\)
\((2)\ Action:\ (to)\ B,C,I\)
\((3)\ B.E.:\)
\(v(I_4)=0\)
\(v(C_t)=\max_{C_t} \{p(C_t)-tc(C_t,C_{t+1}\}+v(C_{t+1}),\ t=1,2,3\)
where \(p(C_t)\) is profit at City C at time t, \(tc(C_t,C_{t+1})\) is the travel cost, \(v(C_t)\) is total profit since time t.
DP Solution
# profit and travel cost: 1-I, 2-B, 3-C
p = c(120,160,170)
tc = matrix(c(0,50,20,
50,0,70,
20,70,0), nrow = 3, byrow = T)
# Place-holders for city and total profit for each day
city = rep(NA,4)
v = rep(0,4)
for (t in 4:1){
if (t==4){
v[t] = 0
city[t] = 1
} else {
v_t = rep(0,3)
for (Ct in 1:3){
v_t[Ct] = p[Ct]-tc[Ct,city[t+1]]
}
v[t] = max(v_t) + v[t+1]
city[t] = which.max(v_t)
}
}
print(v)## [1] 490 320 150 0print(city)## [1] 3 3 3 1Therefore, the maximal profit is 490, which is achieved by traveling to Chicago on the first day, staying there until traveling to Indianapolis on the last day.
Problem 4
DP Solution
r3 = function(d3){
if (d3==0){
return(0)
} else {
return (4*d3+5)
}
}
r2 = function(d2){
if (d2==0){
return(0)
} else {
return (3*d2+7)
}
}
r1 = function(d1){
if (d1==0){
return(0)
} else {
return (7*d1+2)
}
}
v = matrix(NA,nrow = 3,ncol = 7)
for (t in 3:1){
if (t==3){
for (s in 0:6){
v[t,s+1]=r3(s)
}
} else {
for (s in 0:6){
v_temp = rep(0,7)
for (d in 0:s){
if (t==2){
v_temp[d+1] = r2(d) + v[t+1,s+1-d]
} else if (t==1){
v_temp = r1(d) + v[t+1,s+1-d]
}
}
#print (v_temp)
v[t,s+1]=max(v_temp)
#print(c(s,v_temp))
}
}
}
colnames(v) = paste('s=',c(0:6),sep = '')
rownames(v) = paste('t=',c(1:3),sep='')
print(v)## s=0 s=1 s=2 s=3 s=4 s=5 s=6
## t=1 0 9 16 23 30 37 44
## t=2 0 10 19 23 27 31 35
## t=3 0 9 13 17 21 25 29print(v['t=1','s=6'])## [1] 44Therefore, the optimal strategy yields $44,000 by investing all $6,000 in NO.1.