• Load Parameters & Libraries
  • Seed Initialization
  • Simulation
    • Building Surface Area
    • Single Vegetable Yield
    • CO2 & Water Usage
    • Plant Costs
    • Total Costs
  • Yields
    • All
    • Tomatoes
    • Green Beans
    • Kale
  • CO2
  • H20
  • Costs by Plant
  • Yield per Plant Costs

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

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*y*b, 30.22, 266.5,148.18) 
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

Single Vegetable Yield

#tomatoes with  modifier of 3 to 5 for vertical ag per lit
tomatoes=rnorm(i*y*b, 11440,2000)*(surface_area*runif(i*y*b, 3,5)+community_area)
#green beans with  modifier of 3 to 5 for vertical ag per lit
greenbeans=rnorm(i*y*b, 4000, 1000*.2)*(surface_area*runif(i*y*b, 3,5)+community_area)
#broccoli with  modifier of 3 to 5 for vertical ag per lit
broccoli=rnorm(i*y*b, 2950,1000)*(surface_area*runif(i*y*b, 3,5)+community_area)
#kale
kale=rnorm(i*y*b, 12940, 2000)*(surface_area*runif(i*y*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*y*b,.37,.1)*.45*all
h20=runif(p*i*y*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

Plant Costs

#total plant costs with surface area converted back to square feet
tom=surface_area*.5/6*runif(i*y*b,3.77, 4.00)/0.00002295684 #.5 plants per square feet, 6 per package,$3.77 base cost, in square feet
gb=surface_area*4/4*runif(i*y*b,3.77, 4.00)/0.00002295684
broc=surface_area*1/6*runif(i*y*b,5.98, 6.50)/0.00002295684
k=surface_area*1/6*runif(i*y*b,3.77, 4.00)/0.00002295684
cost=array(rbind(tom,gb,broc,k), c(4,b*y,i))
yieldperdollar=all/as.vector(cost)

Total Costs

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

Yields

All

scale=1000000
all=array(all, c(4,b*y,i))/scale
p1=colMeans(all[1,,])
p2=colMeans(all[2,,])
p3=colMeans(all[3,,])
p4=colMeans(all[4,,])
Yield=c(p1,p2,p3,p4)
Iteration=rep(seq(1:30),4)
Plant=c(rep("Tomatoes",30),rep("Green Beans",30),rep("Broccoli",30),rep("Kale",30))

mydf=data.frame(cbind(Yield, Iteration, Plant))
mydf$Yield=as.numeric(mydf$Yield)
mydf$Iteration=as.numeric(mydf$Iteration)

ggplot(mydf, aes(x=Iteration, y=Yield, group=Plant, color=Plant)) + 
  geom_line()+ylab('Yield in Million Pounds')+ggtitle('Yield in Million Pounds')

Tomatoes

myf=function(x) sd(x)/sqrt(length(x))
tomatoes=matrix(tomatoes, ncol=i)
Totals=as.data.frame(colMeans(tomatoes)/(y*scale))
colnames(Totals)="Total"
myse=apply(tomatoes/(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 Tomato Yield in Million lbs per Year per Building 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')

Green Beans

myf=function(x) sd(x)/sqrt(length(x))
scale=1000000
greenbeans=matrix(greenbeans, ncol=i)
Totals=as.data.frame(colMeans(greenbeans)/(y*scale))
colnames(Totals)="Total"
myse=apply(greenbeans/(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 Green Bean Yield in Million lbs per Year per Building 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')

Kale

myf=function(x) sd(x)/sqrt(length(x))
scale=1000000
kale=matrix(kale, ncol=i)
Totals=as.data.frame(colMeans(kale)/(y*scale))
colnames(Totals)="Total"
myse=apply(kale/(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 Kale Yield in Million lbs per Year per Building 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')

## Broccoli

myf=function(x) sd(x)/sqrt(length(x))
scale=1000000
broccoli=matrix(broccoli, ncol=i)
Totals=as.data.frame(colMeans(broccoli)/(y*scale))
colnames(Totals)="Total"
myse=apply(broccoli/(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 Broccoli Yield in Million lbs per Year per Building 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)/p
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))
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')

H20

#https://www.edengreen.com/blog-collection/water-crisis-drought
h20=matrix(h20, ncol=i)/p
Totals=as.data.frame(colMeans(h20)/(y*scale))
colnames(Totals)="Total"
myse=apply(h20/(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))
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')

Costs by Plant

p1=colMeans(cost[1,,])
p2=colMeans(cost[2,,])
p3=colMeans(cost[3,,])
p4=colMeans(cost[4,,])
Cost=c(p1,p2,p3,p4)
Iteration=rep(seq(1:30),4)
Plant=c(rep("Tomatoes",30),rep("Green Beans",30),rep("Broccoli",30),rep("Kale",30))
mydf=data.frame(cbind(Cost, Iteration, Plant))
mydf$Cost=as.numeric(mydf$Cost)
mydf$Iteration=as.numeric(mydf$Iteration)
ggplot(mydf, aes(x=Iteration, y=Cost, group=Plant, color=Plant)) + 
  geom_line()+ylab('Total Plant Costs')+ggtitle('Total Plant Costs in $')+
   scale_y_continuous(labels = function(x) paste0('$',x))

# Costs by Yield

tc=array(tc,c(4,b*y,i))
p1=colMeans(tc[1,,])
p2=colMeans(tc[2,,])
p3=colMeans(tc[3,,])
p4=colMeans(tc[4,,])
tc=c(p1,p2,p3,p4)
Iteration=rep(seq(1:30),4)
Plant=c(rep("Tomatoes",30),rep("Green Beans",30),rep("Broccoli",30),rep("Kale",30))
mydf=data.frame(cbind(tc, Iteration, Plant))
mydf$tc=as.numeric(mydf$tc)
mydf$Iteration=as.numeric(mydf$Iteration)
ggplot(mydf, aes(x=Iteration, y=tc, group=Plant, color=Plant)) + 
  geom_line()+ylab('Total Costs, Millions')+ggtitle('Total Costs in Millions')+
   scale_y_continuous(labels = function(x) paste0('$',x))

Yield per Plant Costs

yieldperdollar=array(yieldperdollar, c(4,b*i,30))
p1=colMeans(yieldperdollar[1,,])
p2=colMeans(yieldperdollar[2,,])
p3=colMeans(yieldperdollar[3,,])
p4=colMeans(yieldperdollar[4,,])
Cost=c(p1,p2,p3,p4)
Iteration=rep(seq(1:30),4)
Plant=c(rep("Tomatoes",30),rep("Green Beans",30),rep("Broccoli",30),rep("Kale",30))
mydf=data.frame(cbind(Cost, Iteration, Plant))
mydf$Cost=as.numeric(mydf$Cost)
mydf$Iteration=as.numeric(mydf$Iteration)
ggplot(mydf, aes(x=Iteration, y=Cost, group=Plant, color=Plant)) + 
  geom_line()+ylab('Yield in Pounds per Dollar of Plant Costs')+ggtitle('Yield in Pounds per $ in Plant Costs')