Load Parameters
#plants = tomatoes, green beans, kale
b=20223 #b=buildings: will replace this with the actual number of building
y=10 #y=years: running for a 10-year lifespan
i=30 #i=iterations: number of iterations
Storage of Results / Seed Initialization
#matrices to hold results
yield=matrix(rep(0,i*b*y), c(i,b*y))
co2=matrix(rep(0,i*b*y), c(i,b*y))
netincome=matrix(rep(0,i*b*y), c(i,b*y))
consumercost=matrix(rep(0,i*b*y), c(i,b*y))
water_use=matrix(rep(0,i*b*y), c(i,b*y))
soil_use=matrix(rep(0,i*b*y), c(i,b*y))
composting_material=matrix(rep(0,i*b*y), c(i,b*y))
energy_use=matrix(rep(0,i*b*y), c(i,b*y))
#seed for replication of results
set.seed(1234)
Height
height=rtri(i*y*b, 20,400, 145) #Triangular
width=runif(i*y*b, 1,3)/10*height #.1 to .3 x height
length=runif(i*y*b, 1,3)/10*height #.1 to .3 x length
fraction=runif(i*y*b,25,75)/100 #25% to 75% of building surface area utilized
fraction2=runif(i*y*b,10,20)/100 #10% to 20% of conservation area utilized
#surface area is LW + 2LH + 2WH..multiply by fraction and convert to acres with constant
surface_area=fraction*(length*width+2*length*height+2*width*height)*0.00002295684
#add to the surface area some fraction that accounts for community land
community_area=fraction2*14311 #fixed conservation acreage
#pounds yield of 35 to 40K lbs per acre with modifier of 3 to 5 for vertical ag per lit
yield=runif(i*y*b, 35000,40000)*(surface_area*runif(i*y*b, 3,5)+community_area) #https://cropcirclefarms.com/farm-yield-calculator.html
#CO2 is about 8 lbs per pound yielded
co2=rnorm(i*y*b,2,.2)*yield
#compare CO2 with typical farming
#https://pubmed.ncbi.nlm.nih.gov/30045579/ ~2 lbs compared to ~8 lbs
#https://www.agritechtomorrow.com/article/2019/11/vertical-farming-methods-compared-to-traditional-vertical-farming/11822
#Cost is between 3 and 10 bucks per pound yielded...Get the numbers $1 to $7 per pound
consumercost=runif(i*y*b,2.6,4)*yield
#https://practicalfarmers.org/research/enterprise-budget-for-cherry-tomatoes-2/
#2.6 to $4.0 per pound to buy it
netincome=runif(i*y*b, 1.79, 2.42)*yield
#https://practicalfarmers.org/research/enterprise-budget-for-cherry-tomatoes-2/
Yield
myf=function(x) sd(x)/sqrt(length(x))
scale=1000000
yield=matrix(yield, ncol=i)
Totals=as.data.frame(colMeans(yield)/(y*scale))
colnames(Totals)="Total"
myse=apply(yield/(y*scale),2,myf)
Totals$Lower=as.data.frame(rep(mean(Totals$Total)-1.96*mean(myse), i))
Totals$Upper=as.data.frame(rep(mean(Totals$Total)+1.96*mean(myse), i))
colnames(Totals)=c("Total", "Lower", "Upper")
plot(Totals$Total~seq(1,i), main="Mean Yield in Million lbs per Year by Iteration", ylab='Yield in Million lbs',type='l', xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
CO2
co2=matrix(co2, ncol=i)
Totals=as.data.frame(colMeans(co2)/(y*scale))
colnames(Totals)="Total"
myse=apply(co2/(y*scale),2,myf)
Totals$Lower=as.data.frame(rep(mean(Totals$Total)-1.96*mean(myse), i))
Totals$Upper=as.data.frame(rep(mean(Totals$Total)+1.96*mean(myse), i))
colnames(Totals)=c("Total", "Lower", "Upper")
plot(Totals$Total~seq(1,i), main="Mean CO2 in Million lbs per Year by Iteration", ylab='Yield in Million lbs',type='l', xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
Consumer Cost
consumercost=matrix(consumercost, ncol=i)
Totals=as.data.frame(colMeans(consumercost)/(y*scale))
colnames(Totals)="Total"
myse=apply(consumercost/(y*scale),2,myf)
Totals$Lower=as.data.frame(rep(mean(Totals$Total)-1.96*mean(myse), i))
Totals$Upper=as.data.frame(rep(mean(Totals$Total)+1.96*mean(myse), i))
colnames(Totals)=c("Total", "Lower", "Upper")
plot(Totals$Total~seq(1,i), main="Mean Consumer Cost in Million $ per Year by Iteration", ylab='Yield in Million lbs',type='l', xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
Net Income
netincome=matrix(netincome, ncol=i)
Totals=as.data.frame(colMeans(netincome)/(y*scale))
colnames(Totals)="Total"
myse=apply(netincome/(y*scale),2,myf)
Totals$Lower=as.data.frame(rep(mean(Totals$Total)-1.96*mean(myse), i))
Totals$Upper=as.data.frame(rep(mean(Totals$Total)+1.96*mean(myse), i))
colnames(Totals)=c("Total", "Lower", "Upper")
plot(Totals$Total~seq(1,i), main="Mean Consumer Cost in Million $ per Year by Iteration", ylab='Yield in Million lbs',type='l', xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
