Load Parameters & Libraries

library(EnvStats)
## 
## Attaching package: 'EnvStats'
## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm
## The following object is masked from 'package:base':
## 
##     print.default
library(ggplot2)

p=4 #tomatoes, green beans, kale, broccoli
b=131536 #b=buildings in Greater Springfield:
y=1 #y=years, year to year
i=30  #i=iterations:  number of iterations
myscale=1000000

Seed Initialization

#seed for replication of results
set.seed(1234)

Simulation

Building Surface Area

#weighted average min mean max = 266.5, 30.22, 148.18 in GSA
  height=rtri(i*b, 30.22, 266.5,148.18) 
  width=runif(i*b, 1,3)/10*height #.1 to .3 x height
  length=runif(i*b, 1,3)/10*height #.1 to .3 x length
fraction=runif(i*b,25,75)/100 #25% to 75% of building surface area utilized
fraction2=runif(i*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
acreage=0.00002295684 #convert from square feet to acres
surface_area=fraction*(length*width+2*length*height+2*width*height)*acreage
#add to the surface area some fraction that accounts for community land
community_area=fraction2*14311 #fixed conservation acreage

Single Vegetable Yield

#Yield is 3 to 5 times greater for vertical agriculture.

#tomatoes with  modifier of 3 to 5 for vertical ag per lit
tomatoes=rnorm(i*b, 11440,2288)*(surface_area*runif(i*b, 3,5)+community_area)
#green beans with  modifier of 3 to 5 for vertical ag per lit
greenbeans=rnorm(i*b, 4000, 800)*(surface_area*runif(i*b, 3,5)+community_area)
#broccoli with  modifier of 3 to 5 for vertical ag per lit
broccoli=rnorm(i*b, 2950,590)*(surface_area*runif(i*b, 3,5)+community_area)
#kale
kale=rnorm(i*b, 12940, 2588)*(surface_area*runif(i*b, 3,5)+community_area)
#all yields
all=as.vector(rbind(tomatoes, greenbeans,broccoli,kale))

#https://cropcirclefarms.com/farm-yield-calculator.html
#https://nevegetable.org/cultural-practices/table-15-approximate-yields

#CO2 is about .37 per kg CO2 eqivalent which is .45 * yield
#https://www.sciencedirect.com/science/article/pii/S0959652616303584?casa_token=2irdiKKNSCYAAAAA:t03z922lYsZqGjP--XRrnXdLi8zxX1YmNipGcuYx4lwF4-ew0VWPSor2sEwpgrfyrLVFBx3J#tbl1

CO2 & Water Usage

co2=rnorm(p*i*b,.31,.1)*.45*all
h20=runif(p*i*b, 0.5, 1)*all

#compare CO2 with typical farming
#https://www.agritechtomorrow.com/article/2019/11/vertical-farming-methods-compared-to-traditional-vertical-farming/11822
#https://www.edengreen.com/blog-collection/water-crisis-drought
#https://agritecture.medium.com/calculating-the-carbon-footprint-of-vertical-farming-and-traditional-farming-e36ff5245d5f

Plant Costs

#total plant costs with surface area converted back to square feet
tom=(surface_area+community_area)*.5/6*runif(i*b,3.77, 1.2*3.77)/acreage #.5 plants per square feet, 6 per package,$3.77 base cost, in square feet
gb=(surface_area+community_area)*4/4*runif(i*b,3.77, 1.2*3.77)/acreage
broc=(surface_area+community_area)*1/6*runif(i*b,5.98, 1.2*5.98)/acreage
k=(surface_area+community_area)*1/6*runif(i*b,3.77, 1.2*3.77)/acreage

Total Costs

tc=all*rtri(i*b*y*p, 3, 4, 3.07)

Yields and Costs

Tomatoes

myf=function(x) sd(x)/sqrt(length(x))
tomatoes=matrix(tomatoes, ncol=i)
Totals=as.data.frame(colMeans(tomatoes))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(myscale)
myse=apply(tomatoes/(myscale),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 Tomato Yield in Million lbs", 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')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: ', round(tmp[1],2))
text(5,tmp[1]+.02, mytext)

Tomato Cost per Building

myf=function(x) sd(x)/sqrt(length(x))
tom=matrix(tom, ncol=i)
Totals=data.frame(colMeans(tom))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(b)
myse=apply(tom/(b),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 Tomato Cost per Building in Dollars", ylab='Dollars',type='l',  xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: $', round(tmp[1],2))
text(5,tmp[1]+.2, mytext)

Green Beans

myf=function(x) sd(x)/sqrt(length(x))
myscale=1000000
greenbeans=matrix(greenbeans, ncol=i)
Totals=data.frame(colMeans(greenbeans))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(myscale)
myse=apply(greenbeans/(myscale),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 Green Bean Yield in Million lbs", 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')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: ', round(tmp[1],2))
text(5,tmp[1]+.002, mytext)

## Green Beans Cost

myf=function(x) sd(x)/sqrt(length(x))
gb=matrix(gb, ncol=i)
Totals=data.frame(colMeans(gb))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(b)
myse=apply(gb/(b),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 Green Beans Costs per Building ", ylab='Dollars',type='l',  xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: $', round(tmp[1],2))
text(5,tmp[1]+1, mytext)

## Kale

myf=function(x) sd(x)/sqrt(length(x))
myscale=1000000
kale=matrix(kale, ncol=i)
Totals=data.frame(colMeans(kale))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(myscale)
myse=apply(kale/(myscale),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 Kale Yield in Million lbs", 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')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: ', round(tmp[1],2))
text(5,tmp[1]+.02, mytext)

Kale Cost

myf=function(x) sd(x)/sqrt(length(x))
k=matrix(k, ncol=i)
Totals=data.frame(colMeans(k))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(b)
myse=apply(k/(b),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 Kale Costs per Building", ylab='Dollars',type='l',  xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: $', round(tmp[1],2))
text(5,tmp[1]+.2, mytext)

Broccoli

myf=function(x) sd(x)/sqrt(length(x))
myscale=1000000
broccoli=matrix(broccoli, ncol=i)
Totals=data.frame(colMeans(broccoli))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(myscale)
myse=apply(broccoli/(myscale),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 Broccoli Yield in Million lbs", 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')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: ', round(tmp[1],2))
text(5,tmp[1]+.002, mytext)

Broccoli Costs

myf=function(x) sd(x)/sqrt(length(x))
broc=matrix(broc, ncol=i)
Totals=data.frame(colMeans(broc))
colnames(Totals)="Total"
Totals$Total=Totals$Total/(b)
myse=apply(broc/(b),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 Broccoli Costs per Building", ylab='Dollars',type='l',  xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: $', round(tmp[1],2))
text(5,tmp[1]+.5, mytext)

CO2

co2=matrix(co2, ncol=i)/p
Totals=data.frame(colMeans(co2)/(myscale))
colnames(Totals)="Total"
myse=apply(co2/(myscale),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))
Totals$Iteration=seq(1:30)
colnames(Totals)=c("Total", "Lower", "Upper", "Iteration")
plot(Totals$Total~seq(1,i), main="Mean CO2 in Million lbs by Iteration, All Plants", ylab='CO2 in Million Pounds',type='l',  xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: ', round(tmp[1],4))
text(5,tmp[1]+.001, mytext)

H20

#https://www.edengreen.com/blog-collection/water-crisis-drought
h20=matrix(h20, ncol=i)/p
Totals=data.frame(colMeans(h20)/(myscale))
colnames(Totals)="Total"
myse=apply(h20/(myscale),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))
Totals$Iteration=seq(1:30)
colnames(Totals)=c("Total", "Lower", "Upper", "Iteration")
plot(Totals$Total~seq(1,i), main="Mean H20 in Gallons by Iteration, All Plants", ylab='H20 in Gallons',type='l',  xlab='Iteration', ylim=c(min(Totals$Lower),max(Totals$Upper)))
lines(Totals$Lower, col='red')
lines(Totals$Upper, col='blue')
tmp=rep(mean(Totals$Total), i)
lines(x=seq(1:30),y=tmp, col='black')
mytext=paste('Mean: ', round(tmp[1],4))
text(5,tmp[1]+.001, mytext)