Brian GPT 02

Initial, Primary Monte Carlo Method Simulation Code

Libraries

require(dplyr)

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

require(triangle)

## Loading required package: triangle

require(ggplot2)

## Loading required package: ggplot2

require(stats)

Initial Parameters

iterations <- 1200
months <- 60
set.seed(321)

Cars per Month (Deterministic, Not in Simulation)

This value is indexed only by month.

# for (i in 1:months) <--- I'm not sure the loop changes anything in the next four lines?
# I even tried adding brackets {} and it didn't do anything different
ice_car_growth_rate <- 1.005^(1:months)  #0.5% per month
ice_truck_growth_rate <- 1.003^(1:months)  # 0.3% per month
bev_car_growth_rate <- 1.02^(1:months)  # 2% per month
bev_truck_growth_rate <- 1.015^(1:months)  # 1.5% per month


# Initial number of vehicles (estimates for the start of the simulation)
initial_ice_cars <- 546000
initial_ice_trucks <- 234000
initial_bev_cars <- 14000
initial_bev_trucks <- 6000

# Set initial values
ice_cars_vector <- initial_ice_cars*ice_car_growth_rate
ice_trucks_vector <- initial_ice_trucks*ice_truck_growth_rate
bev_cars_vector <- initial_bev_cars*bev_car_growth_rate
bev_trucks_vector <- initial_bev_trucks*bev_truck_growth_rate

Simulation

# Monthly Miles
set.seed(1)
miles <- rtriangle(iterations*months, 4785, 18858, 13476)/12

# Average battery size in kWh for BEVs
set.seed(2)
bev_battery_size_kwh <- rtriangle(iterations*months, 40, 118, 82)  

#Average battery size in kwh for BEV trucks
set.seed(3)
bev_truck_battery_size_kwh <- rtriangle(iterations*months, 92, 212, 135) 

# CO2 emission factors for ICE cars (sedans + small crossover SUVs) 
set.seed(4)
truck_co2_per_10k_miles <- rnorm(iterations*months, mean = (4.2 + 5.8) / 2, sd = (5.8 - 4.2) / 4)  

# CO2 emission factors for ICE trucks (Full sized pickups and larger SUVs)
set.seed(5)
car_co2_per_10k_miles <- rnorm(iterations*months, mean = (2.8 + 5.5) / 2, sd = (5.5 - 2.8) / 4)  # Normal distribution between 2.8 and 5.5 tons

# CO2 emissions per kWh for battery production using triangular distribution
set.seed(6)  # Set seed for reproducibility
bev_battery_co2_per_kwh <- rtriangle(iterations*months, 0.05, 0.2, 0.14)  # tons of CO2 per kWh

# CO2 for battery production
bev_battery_co2 <- bev_battery_co2_per_kwh * bev_battery_size_kwh  # Total CO2 for battery production
bev_truck_battery_co2 <- bev_battery_co2_per_kwh * bev_truck_battery_size_kwh #Total CO2 for battery production (trucks)

# Energy consumption values
set.seed(10) #SET SEEDS FOR ALL DIST
ice_car_consumption <- rtriangle(iterations*months, 3e6, 5e6, 4e6) / 3600000 * miles
set.seed(11)
ice_truck_consumption <- rtriangle(iterations*months, 6e6, 8e6, 7e6) / 3600000 * miles
set.seed(12)
# OLD: bev_car_consumption <- runif(iterations*months, 720000, 1260000) / 3600000 * miles
bev_car_consumption <- rnorm(iterations * months, mean = 990000, sd = 156174) / 3600000 * miles
set.seed(13)
bev_truck_consumption <- rnorm(iterations*months, 1800000, 120000) / 3600000 * miles
set.seed(14)
bev_car_grid_demand <- 1/rtriangle(iterations*months, 0.6, 0.89, 0.8)
set.seed(15)
bev_truck_grid_demand <- 1/rtriangle(iterations*months, 0.6, 0.89, 0.8)
bev_car_energy <- bev_car_consumption * (bev_car_grid_demand) 
bev_truck_energy <- bev_truck_consumption * bev_truck_grid_demand 
  
# CO2 emissions based on miles drive
ice_car_co2 <- (miles / 10000) * car_co2_per_10k_miles
ice_truck_co2 <- (miles / 10000) * truck_co2_per_10k_miles

#Grid CO2 production as a grid emission factor
grid_emissions_factor <- 0.000233  # Tons of CO₂ per kWh

# Calculate incremental vehicle populations
bev_cars_increment <- c(bev_cars_vector[1], diff(bev_cars_vector))  # Growth in BEV cars each month
bev_trucks_increment <- c(bev_trucks_vector[1], diff(bev_trucks_vector))  # Growth in BEV trucks each month

# Matrix for monthly BEV battery production emissions
bev_car_battery_production_matrix <- matrix(
  bev_battery_co2_per_kwh * bev_battery_size_kwh,  # CO2 per vehicle
  nrow = months,
  ncol = iterations
) * bev_cars_increment  # Scale by incremental population growth

bev_truck_battery_production_matrix <- matrix(
  bev_battery_co2_per_kwh * bev_truck_battery_size_kwh,  # CO2 per vehicle
  nrow = months,
  ncol = iterations
) * bev_trucks_increment  # Scale by incremental population growth

# Cumulative battery emissions for BEVs
bev_car_cumulative_battery_co2 <- apply(bev_car_battery_production_matrix, 2, cumsum)
bev_truck_cumulative_battery_co2 <- apply(bev_truck_battery_production_matrix, 2, cumsum)

# Grid-Based CO2 (unchanged from your existing logic)
bev_car_grid_co2 <- matrix(bev_car_energy * bev_cars_vector * grid_emissions_factor, nrow = months, ncol = iterations)
bev_truck_grid_co2 <- matrix(bev_truck_energy * bev_trucks_vector * grid_emissions_factor, nrow = months, ncol = iterations)

# Combine battery production and grid-based emissions
bev_car_co2_matrix <- bev_car_cumulative_battery_co2 + apply(bev_car_grid_co2, 2, cumsum)
bev_truck_co2_matrix <- bev_truck_cumulative_battery_co2 + apply(bev_truck_grid_co2, 2, cumsum)


#Cumulative Miles 

mat_miles <- matrix(miles, nrow = months)
cumulative_miles <- apply(mat_miles, 2, cumsum)  # Cumulative miles for each iteration
cumulative_miles_mean <- rowMeans(cumulative_miles)

Results

ice_car_energy_matrix=matrix(ice_car_consumption * ice_cars_vector, nrow=months)
ice_truck_energy_matrix=matrix(ice_truck_consumption * ice_trucks_vector, nrow=months)
bev_car_energy_matrix=matrix(bev_car_energy * bev_cars_vector, nrow=months)
bev_truck_energy_matrix=matrix(bev_truck_energy * bev_trucks_vector, nrow=months)
ice_car_co2_matrix=matrix(ice_car_co2 * ice_cars_vector, nrow=months)
ice_truck_co2_matrix=matrix(ice_truck_co2 * ice_trucks_vector, nrow=months)
bev_car_co2_matrix=matrix(cumsum(bev_car_co2_matrix), nrow=months)
bev_truck_co2_matrix=matrix(cumsum(bev_truck_co2_matrix), nrow=months)

Column Means, SDs, and CIs

# Calculate mean, standard error, and 95% confidence intervals
myf <- function(x) {
  means = rowMeans(x)  # Row means
  se = apply(x, 1, sd) / sqrt(iterations)  # Row standard error
  lower = means - qt(0.975, iterations - 1) * se  # Lower bound
  upper = means + qt(0.975, iterations - 1) * se  # Upper bound
  mymat = cbind(lower, means, upper)  # Combine into matrix
  return(mymat)
  
}

# ICE Car Energy Scatter
ice_energy_means=myf(ice_car_energy_matrix)
plot(ice_energy_means)

#Start Making Visual Plots #ICE truck Energy Scatter

ice_truck_energy_means=myf(ice_truck_energy_matrix)
plot(ice_truck_energy_means)

#BEV Car Energy Scatter

bev_car_energy_means=myf(bev_car_energy_matrix)
plot(bev_car_energy_means)

#BEV Truck Energy Scatter

bev_truck_energy_means=myf(bev_truck_energy_matrix)
plot(bev_truck_energy_means)

#ICE CO2 Scatter

ice_car_co2_means=myf(ice_car_co2_matrix)
plot(ice_car_co2_means)

ICE Truck CO2 Scatter

ice_truck_co2_means=myf(ice_truck_co2_matrix)
plot(ice_truck_co2_means)

#BEV Car CO2 Scatter

bev_car_co2_means=myf(bev_car_co2_matrix)
plot(bev_car_co2_means)

BEV Truck CO2 Scatter

bev_truck_co2_means=myf(bev_truck_co2_matrix)
plot(bev_truck_co2_means)

#Run CI function, create dataframe, and plot BEV energy over miles

# Apply the `myf` function to the BEV car energy matrix
bev_car_energy_ci <- myf(bev_car_energy_matrix)

# Create a data frame for plotting
bev_energy_plot_data <- data.frame(
  Miles = cumulative_miles_mean,  # Use cumulative miles mean or other suitable x-axis values
  Mean = bev_car_energy_ci[, "means"],
  Lower_CI = bev_car_energy_ci[, "lower"],
  Upper_CI = bev_car_energy_ci[, "upper"]
)

# Plot using ggplot2
library(ggplot2)

ggplot(bev_energy_plot_data, aes(x = Miles)) +
  geom_line(aes(y = Mean), color = "green") +  # Mean line
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "green", alpha = 0.2) +  # Confidence intervals
  labs(
    title = "BEV Car Energy Consumption with Confidence Intervals",
    x = "Cumulative Miles Driven",
    y = "Energy Consumption (kWh)"
  ) +
  theme_minimal()

# Combine ICE cars and trucks energy matrices
ice_total_energy_matrix <- ice_car_energy_matrix + ice_truck_energy_matrix

# Apply the `myf` function to the combined ICE energy matrix
ice_energy_ci <- myf(ice_total_energy_matrix)

# Create a data frame for plotting
ice_energy_plot_data <- data.frame(
  Miles = cumulative_miles_mean,  # Use cumulative miles mean for x-axis
  Mean = ice_energy_ci[, "means"],
  Lower_CI = ice_energy_ci[, "lower"],
  Upper_CI = ice_energy_ci[, "upper"]
)

# Plot using ggplot2
library(ggplot2)

ggplot(ice_energy_plot_data, aes(x = Miles)) +
  geom_line(aes(y = Mean), color = "blue") +  # Mean line for ICE vehicles
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "blue", alpha = 0.2) +  # Confidence intervals
  labs(
    title = "ICE Vehicle Energy Consumption with Confidence Intervals",
    x = "Cumulative Miles Driven",
    y = "Energy Consumption (kWh)"
  ) +
  theme_minimal()

# Combine BEV cars and trucks CO2 matrices
bev_total_co2_matrix <- bev_car_co2_matrix + bev_truck_co2_matrix

# Apply the `myf` function to the combined BEV CO2 matrix
bev_total_co2_ci <- myf(bev_total_co2_matrix)

# Create a data frame for plotting
bev_total_co2_plot_data <- data.frame(
  Miles = cumulative_miles_mean,  # Use cumulative miles mean for x-axis
  Mean = bev_total_co2_ci[, "means"],  # Mean of total BEV CO2
  Lower_CI = bev_total_co2_ci[, "lower"],  # Lower confidence interval
  Upper_CI = bev_total_co2_ci[, "upper"]  # Upper confidence interval
)

# Plot BEV CO2 emissions with confidence intervals over cumulative miles
ggplot(bev_total_co2_plot_data, aes(x = Miles)) +
  geom_line(aes(y = Mean), color = "red", size = 1) +  # Mean line for BEV CO2 emissions
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "red", alpha = 0.2) +  # Confidence intervals
  labs(
    title = "Total BEV CO2 Emissions with Confidence Intervals",
    x = "Cumulative Miles Driven",
    y = "CO2 Emissions (Tons)"
  ) +
  theme_minimal() +
  scale_x_continuous(limits = c(0, 60000)) +  # Set x-axis range to 0 - 60000 miles
  scale_y_continuous(limits = c(1.7e+10, 2.5e+10),
    name = "CO2 Emissions (Tons)"
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Combine ICE cars and trucks CO2 matrices
ice_total_co2_matrix <- ice_car_co2_matrix + ice_truck_co2_matrix

# Apply the `myf` function to the combined ICE CO2 matrix
ice_co2_ci <- myf(ice_total_co2_matrix)

# Create a data frame for plotting
ice_co2_plot_data <- data.frame(
  Miles = cumulative_miles_mean,  # Use cumulative miles mean for x-axis
  Mean = ice_co2_ci[, "means"],
  Lower_CI = ice_co2_ci[, "lower"],
  Upper_CI = ice_co2_ci[, "upper"]
)

# Plot using ggplot2
library(ggplot2)

ggplot(ice_co2_plot_data, aes(x = Miles)) +
  geom_line(aes(y = Mean), color = "blue") +  # Mean line for ICE CO2 emissions
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill = "blue", alpha = 0.2) +  # Confidence intervals
  labs(
    title = "ICE CO2 Emissions with Confidence Intervals",
    x = "Cumulative Miles Driven",
    y = "CO2 Emissions (Tons)"
  ) +
  theme_minimal() +
  scale_x_continuous(limits = c(0, max(cumulative_miles_mean)))  # Ensure x-axis shows the full range

# Combine BEV and ICE total CO2 data into a data frame
co2_comparison_data <- data.frame(
  Miles = cumulative_miles_mean,  # Shared x-axis for cumulative miles
  BEV_CO2_Mean = bev_total_co2_ci[, "means"],  # BEV CO2 means
  BEV_CO2_Lower_CI = bev_total_co2_ci[, "lower"],
  BEV_CO2_Upper_CI = bev_total_co2_ci[, "upper"],
  ICE_CO2_Mean = ice_co2_ci[, "means"],  # ICE CO2 means
  ICE_CO2_Lower_CI = ice_co2_ci[, "lower"],
  ICE_CO2_Upper_CI = ice_co2_ci[, "upper"]
)

# Plot BEV and ICE total CO2 on the same chart
ggplot(co2_comparison_data, aes(x = Miles)) +
  # BEV CO2 line and ribbon
  geom_line(aes(y = BEV_CO2_Mean, color = "BEV"), size = 1) +
  geom_ribbon(aes(ymin = BEV_CO2_Lower_CI, ymax = BEV_CO2_Upper_CI, fill = "BEV"), alpha = 0.2) +
  # ICE CO2 line and ribbon
  geom_line(aes(y = ICE_CO2_Mean, color = "ICE"), size = 1) +
  geom_ribbon(aes(ymin = ICE_CO2_Lower_CI, ymax = ICE_CO2_Upper_CI, fill = "ICE"), alpha = 0.2) +
  # Labels and legend
  labs(
    title = "Total CO2 Emissions: BEV vs ICE",
    x = "Miles Driven",
    y = "CO2 Emissions (Tons)",
    color = "Vehicle Type",
    fill = "Vehicle Type"
  ) +
  scale_color_manual(values = c("BEV" = "green", "ICE" = "blue")) +
  scale_fill_manual(values = c("BEV" = "green", "ICE" = "blue")) +
  theme_minimal()

# Combine BEV and ICE total CO2 data into a data frame
co2_comparison_data <- data.frame(
  Miles = cumulative_miles_mean,  # Shared x-axis for cumulative miles
  BEV_CO2_Mean = bev_total_co2_ci[, "means"],  # BEV CO2 means
  BEV_CO2_Lower_CI = bev_total_co2_ci[, "lower"],
  BEV_CO2_Upper_CI = bev_total_co2_ci[, "upper"],
  ICE_CO2_Mean = ice_co2_ci[, "means"],  # ICE CO2 means
  ICE_CO2_Lower_CI = ice_co2_ci[, "lower"],
  ICE_CO2_Upper_CI = ice_co2_ci[, "upper"]
)

# Plot BEV and ICE total CO2 on the same chart with log scale
ggplot(co2_comparison_data, aes(x = Miles)) +
  # BEV CO2 line and ribbon
  geom_line(aes(y = BEV_CO2_Mean, color = "BEV"), size = 1) +
  geom_ribbon(aes(ymin = BEV_CO2_Lower_CI, ymax = BEV_CO2_Upper_CI, fill = "BEV"), alpha = 0.2) +
  # ICE CO2 line and ribbon
  geom_line(aes(y = ICE_CO2_Mean, color = "ICE"), size = 1) +
  geom_ribbon(aes(ymin = ICE_CO2_Lower_CI, ymax = ICE_CO2_Upper_CI, fill = "ICE"), alpha = 0.2) +
  # Labels and legend
  labs(
    title = "Total CO2 Emissions (Log Scale): BEV vs ICE",
    x = "Miles Driven",
    y = "CO2 Emissions (Tons, Log Scale)",
    color = "Vehicle Type",
    fill = "Vehicle Type"
  ) +
  scale_y_log10() +  # Apply log scale to y-axis
  scale_color_manual(values = c("BEV" = "green", "ICE" = "blue")) +
  scale_fill_manual(values = c("BEV" = "green", "ICE" = "blue")) +
  theme_minimal()

# Combine BEV and ICE total CO2 data into a data frame
energy_comparison_data <- data.frame(
  Miles = cumulative_miles_mean,  # Shared x-axis for cumulative miles
  BEV_energy_Mean = bev_car_energy_ci[, "means"],  # BEV CO2 means
  BEV_energy_Lower_CI = bev_car_energy_ci[, "lower"],
  BEV_energy_Upper_CI = bev_car_energy_ci[, "upper"],
  ICE_energy_Mean = ice_energy_ci[, "means"],  # ICE CO2 means
  ICE_energy_Lower_CI = ice_energy_ci[, "lower"],
  ICE_energy_Upper_CI = ice_energy_ci[, "upper"]
)

# Plot BEV and ICE total CO2 on the same chart with log scale
ggplot(energy_comparison_data, aes(x = Miles)) +
  # BEV CO2 line and ribbon
  geom_line(aes(y = BEV_energy_Mean, color = "BEV"), size = 1) +
  geom_ribbon(aes(ymin = BEV_energy_Lower_CI, ymax = BEV_energy_Upper_CI, fill = "BEV"), alpha = 0.2) +
  # ICE CO2 line and ribbon
  geom_line(aes(y = ICE_energy_Mean, color = "ICE"), size = 1) +
  geom_ribbon(aes(ymin = ICE_energy_Lower_CI, ymax = ICE_energy_Upper_CI, fill = "ICE"), alpha = 0.2) +
  # Labels and legend
  labs(
    title = "Total Energy Consumption (Log Scale): BEV vs ICE",
    x = "Miles Driven",
    y = "CO2 Emissions (Tons, Log Scale)",
    color = "Vehicle Type",
    fill = "Vehicle Type"
  ) +
  scale_y_log10() +  # Apply log scale to y-axis
  scale_color_manual(values = c("BEV" = "red", "ICE" = "purple")) +
  scale_fill_manual(values = c("BEV" = "red", "ICE" = "purple")) +
  theme_minimal()

# Create a histogram of monthly miles
hist(miles, 
     main = "Histogram of Monthly Miles Driven", 
     xlab = "Miles Driven (Monthly)", 
     col = "skyblue", 
     breaks = 30,  # Number of bins
     border = "white")

# Create side-by-side histograms
par(mfrow = c(1, 2))  # Set up a 1x2 plotting layout

# Histogram for BEV battery size
hist(bev_battery_size_kwh,
     main = "Histogram of BEV Battery Sizes",
     xlab = "Battery Size (kWh)",
     col = "lightgreen",
     breaks = 30,  # Number of bins
     border = "white")

# Histogram for BEV Truck battery size
hist(bev_truck_battery_size_kwh,
     main = "Histogram of BEV Truck Battery Sizes",
     xlab = "Battery Size (kWh)",
     col = "lightblue",
     breaks = 30,  # Number of bins
     border = "white")

# Reset plotting layout to default
par(mfrow = c(1, 1))

# Create side-by-side histograms
par(mfrow = c(1, 2))  # Set up a 1x2 plotting layout

# Histogram for Truck CO2 emissions per 10k miles
hist(truck_co2_per_10k_miles,
     main = "ICE Truck CO2 Emissions (per 10k miles)",
     xlab = "CO2 Emissions (Tons)",
     col = "lightblue",
     breaks = 30,  # Number of bins
     border = "white")

# Histogram for Car CO2 emissions per 10k miles
hist(car_co2_per_10k_miles,
     main = "ICE Car CO2 Emissions (per 10k miles)",
     xlab = "CO2 Emissions (Tons)",
     col = "lightgreen",
     breaks = 30,  # Number of bins
     border = "white")

# Reset plotting layout to default
par(mfrow = c(1, 1))

# Histogram for BEV battery CO2 emissions per kWh
hist(bev_battery_co2_per_kwh,
     main = "BEV Battery CO2 Emissions per kWh",
     xlab = "CO2 Emissions (Tons per kWh)",
     col = "purple",
     breaks = 30,  # Number of bins
     border = "black")

# First pair: ICE car and truck consumption
par(mfrow = c(1, 2))  # Arrange plots side by side
hist(ice_car_consumption,
     main = "ICE Car Consumption",
     xlab = "Consumption (kWh)",
     col = "skyblue",
     breaks = 30,
     border = "white")
hist(ice_truck_consumption,
     main = "ICE Truck Consumption",
     xlab = "Consumption (kWh)",
     col = "skyblue",
     breaks = 30,
     border = "white")

# Second pair: BEV car and truck consumption
par(mfrow = c(1, 2))  # Arrange plots side by side
hist(bev_car_consumption,
     main = "BEV Car Consumption",
     xlab = "Consumption (kWh)",
     col = "lightgreen",
     breaks = 30,
     border = "white")
hist(bev_truck_consumption,
     main = "BEV Truck Consumption",
     xlab = "Consumption (kWh)",
     col = "lightgreen",
     breaks = 30,
     border = "white")

# Third pair: BEV grid demand for cars and trucks
par(mfrow = c(1, 2))  # Arrange plots side by side
hist(bev_car_grid_demand,
     main = "BEV Car Grid Demand",
     xlab = "Grid Demand (Inverse Efficiency)",
     col = "orange",
     breaks = 30,
     border = "white")
hist(bev_truck_grid_demand,
     main = "BEV Truck Grid Demand",
     xlab = "Grid Demand (Inverse Efficiency)",
     col = "orange",
     breaks = 30,
     border = "white")

# Reset plotting layout to default
par(mfrow = c(1, 1))

print(dim(bev_car_battery_production_matrix))  # Should be months x iterations

## [1]   60 1200

print(dim(bev_truck_battery_production_matrix))  # Should be months x iterations

## [1]   60 1200

print(dim(bev_car_grid_co2))  # Should be months x iterations

## [1]   60 1200

print(dim(bev_truck_grid_co2))  # Should be months x iterations

## [1]   60 1200

print(dim(bev_car_co2_matrix))  # Should be months x iterations

## [1]   60 1200

print(dim(bev_truck_co2_matrix))  # Should be months x iterations

## [1]   60 1200

******************************************************************************

***** Single Vehicle Monte Carlo Method Code *****

******************************************************************************

# Retain all previous libraries and settings
# Load necessary libraries
library(dplyr)
library(triangle)
library(ggplot2)
library(stats)

# Set global seed for reproducibility
set.seed(123)

Initialize the number of months and iterations to run

# Set number of iterations
iterations <- 1000
months <- 60

Initialize matrices to stor miles and emissions for each iteration

# Initialize arrays to store cumulative miles and CO₂ emissions per iteration
cumulative_miles_iterations <- matrix(0, nrow = months, ncol = iterations)
bev_car_co2_iterations <- matrix(0, nrow = months, ncol = iterations)
ice_car_co2_iterations <- matrix(0, nrow = months, ncol = iterations)

Set global seed

# Run simulation across iterations
set.seed(123)  # Ensure reproducibility

Run the actual simulaiton using “for” loops instead of native matrices

for (i in 1:iterations) {
  # Sample annual and monthly mileage for this iteration
  set.seed(1)
  annual_miles_singular <- rtriangle(1, 4785, 18858, 13476)
  monthly_miles_singular <- annual_miles_singular / 12
  
  # Initialize cumulative mileage and emissions for this iteration
  cumulative_miles <- numeric(months)
  bev_car_co2 <- numeric(months)
  ice_car_co2 <- numeric(months)
  
  # Apply initial battery production CO₂ for BEV
  set.seed(2)
  bev_battery_size_kwh <- rtriangle(1, 45, 135, 82)
  bev_battery_co2_per_kwh <- rtriangle(1, 0.05, 0.2, 0.14)
  bev_battery_co2 <- bev_battery_size_kwh * bev_battery_co2_per_kwh
  bev_car_co2[1] <- bev_battery_co2
  
  # Initial ICE emissions
  set.seed(3)
  car_co2_per_10k_miles <- rnorm(1, mean = (2.6 + 5.2) / 2, sd = (5.2 - 2.6) / 4)
  cumulative_miles[1] <- monthly_miles_singular
  ice_car_co2[1] <- (monthly_miles_singular / 10000) * car_co2_per_10k_miles
  
  # Loop over months
  for (month in 2:months) {
    cumulative_miles[month] <- cumulative_miles[month - 1] + monthly_miles_singular
    
    # BEV energy consumption and grid emissions
    set.seed(4)
    bev_energy_consumption_per_mile <- runif(1, 0.24, 0.35)
    bev_energy_consumption <- bev_energy_consumption_per_mile * monthly_miles_singular
    bev_car_co2[month] <- bev_car_co2[month - 1] + (bev_energy_consumption * 0.000233)
    
    # ICE emissions
    ice_car_co2[month] <- ice_car_co2[month - 1] + (monthly_miles_singular / 10000) * car_co2_per_10k_miles
  }
  
  # Store results in iteration matrices
  cumulative_miles_iterations[, i] <- cumulative_miles
  bev_car_co2_iterations[, i] <- bev_car_co2
  ice_car_co2_iterations[, i] <- ice_car_co2
}

Calcualte and store miles and emission results across all iterations

# Calculate mean cumulative miles and emissions across all iterations
mean_cumulative_miles <- rowMeans(cumulative_miles_iterations)
mean_bev_car_co2 <- rowMeans(bev_car_co2_iterations)
mean_ice_car_co2 <- rowMeans(ice_car_co2_iterations)

Manually calculate confidencie intervals instead of creating a funciton

# Optionally calculate confidence intervals (95%)
bev_car_co2_se <- apply(bev_car_co2_iterations, 1, sd) / sqrt(iterations)
ice_car_co2_se <- apply(ice_car_co2_iterations, 1, sd) / sqrt(iterations)

bev_car_co2_lower <- mean_bev_car_co2 - qt(0.975, iterations - 1) * bev_car_co2_se
bev_car_co2_upper <- mean_bev_car_co2 + qt(0.975, iterations - 1) * bev_car_co2_se

ice_car_co2_lower <- mean_ice_car_co2 - qt(0.975, iterations - 1) * ice_car_co2_se
ice_car_co2_upper <- mean_ice_car_co2 + qt(0.975, iterations - 1) * ice_car_co2_se

Create data frames and plot the data

# Create a data frame for plotting
df_bev_vs_ice_avg <- data.frame(
  Miles = mean_cumulative_miles,
  BEV_CO2_Mean = mean_bev_car_co2,
  BEV_CO2_Lower = bev_car_co2_lower,
  BEV_CO2_Upper = bev_car_co2_upper,
  ICE_CO2_Mean = mean_ice_car_co2,
  ICE_CO2_Lower = ice_car_co2_lower,
  ICE_CO2_Upper = ice_car_co2_upper
)


library(ggplot2)

ggplot(df_bev_vs_ice_avg, aes(x = Miles)) +
  # BEV Emissions
  geom_line(aes(y = BEV_CO2_Mean), color = "green", linetype = "solid") +
  geom_ribbon(aes(ymin = BEV_CO2_Lower, ymax = BEV_CO2_Upper), fill = "green", alpha = 0.2) +
  
  # ICE Emissions
  geom_line(aes(y = ICE_CO2_Mean), color = "blue", linetype = "dashed") +
  geom_ribbon(aes(ymin = ICE_CO2_Lower, ymax = ICE_CO2_Upper), fill = "blue", alpha = 0.2) +
  
  # Plot Labels
  labs(
    title = "Average CO₂ Emissions: BEV vs ICE (Multiple Iterations)",
    x = "Cumulative Miles Driven",
    y = "CO₂ Emissions (Tons)",
    caption = "Green Solid: BEV (including battery production CO₂ at zero miles)\nBlue Dashed: ICE"
  ) +
  
  theme_minimal() +
  xlim(0, max(mean_cumulative_miles)) +
  ylim(0, max(c(mean_bev_car_co2, mean_ice_car_co2)))

******************************************************************************

***** RSM Three Variable Code *****

******************************************************************************

Install library and build three scales of input data

# Install necessary libraries if not already installed
# install.packages("plotly")
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# Generating data (BEV battery size, ICE CO2 per 10,000 miles, and Miles to Offset)
bev_battery_size <- seq(40, 212, length.out = 30)  # BEV battery size from 40 to 212 kWh
ice_co2_output <- seq(2.1, 8, length.out = 30)    # ICE CO2 output from 2.1 to 8 tons
co2_per_kwh <- seq(0.05, 0.2, length.out = 10)    # CO2 per kWh of battery production

Create a grid for BEV battery size and ICE CO2 per mile

# Create grid of values for battery size and ICE CO2
battery_ice_grid <- expand.grid(Battery = bev_battery_size, ICE_CO2 = ice_co2_output)

# Function to calculate Miles to Offset
miles_to_offset <- function(battery_size, ice_co2, co2_kwh) {
  # Constants
  co2_per_mile_bev <- 0.25 * 0.000233  # BEV CO2 per mile (fixed)
  
  # BEV battery CO2 emissions
  bev_battery_co2 <- battery_size * co2_kwh
  
  # Calculate miles needed to offset
  miles <- bev_battery_co2 / (ice_co2 / 10000 - co2_per_mile_bev)
  return(miles)
}

# Prepare an empty list to store plot data
plot_data <- list()

Create Surface Plots for a third variable

# Generate surface plots for different values of CO2 per kWh (third variable)
for (co2_kwh in co2_per_kwh) {
  # Calculate Miles to Offset for each combination of BEV battery size and ICE CO2
  battery_ice_grid$Miles <- with(battery_ice_grid, miles_to_offset(Battery, ICE_CO2, co2_kwh))
  
  # Add plot data for each CO2 per kWh scenario
  plot_data[[as.character(co2_kwh)]] <- plot_ly(
    x = ~battery_ice_grid$Battery,
    y = ~battery_ice_grid$ICE_CO2,
    z = ~battery_ice_grid$Miles,
    type = "scatter3d",
    mode = "markers",
    marker = list(size = 3, color = ~battery_ice_grid$Miles, colorscale = "Viridis"),
    text = paste("BEV Battery:", battery_ice_grid$Battery, "kWh", "<br>ICE CO2 Output:", battery_ice_grid$ICE_CO2,
                 "<br>Miles to Offset:", battery_ice_grid$Miles, "miles")
  )
}

# Render the first plot for one CO2 per kWh value
plot_data[[1]]

******************************************************************************

***** Chi-Square Testing *****

******************************************************************************

Triangular: BEV Car Battery Size

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 40  # Minimum value
max_val <- 118  # Maximum value
mode_val <- 82  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 4.0981, df = 9, p-value = 0.9048

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("skyblue", "orange"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Triangular: BEV Truck Battery Size

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 92  # Minimum value
max_val <- 212  # Maximum value
mode_val <- 135  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 0.70284, df = 9, p-value = 0.9999

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("blue", "yellow"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Triangular: Annual Operating Mileage (DOT, 2022)

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 4785  # Minimum value
max_val <- 18858  # Maximum value
mode_val <- 13476  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 2.4652, df = 9, p-value = 0.9818

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("skyblue", "orange"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Triangular: Upstream BEV Battery Emissions (Tons CO2 per kWh)

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 0.05  # Minimum value
max_val <- 0.2  # Maximum value
mode_val <- 0.14  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 3.0583, df = 9, p-value = 0.9619

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("blue", "yellow"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Triangular: ICE Car Energy Consumption

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 3000000  # Minimum value
max_val <- 5000000  # Maximum value
mode_val <- 4000000  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 2.9368, df = 9, p-value = 0.9667

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("lightblue", "pink"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Triangular: ICE Truck Energy Consumption

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 6000000  # Minimum value
max_val <- 8000000  # Maximum value
mode_val <- 7000000  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 2.9368, df = 9, p-value = 0.9667

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("blue", "yellow"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Triangular: Grid & Charing Efficiency Factor

# Load the triangle package
if (!require(triangle)) {
  install.packages("triangle")
}
library(triangle)

# Set seed for reproducibility
set.seed(123)

# Generate random samples from a triangular distribution
n <- 1000  # Number of samples
min_val <- 0.6  # Minimum value
max_val <- 0.89  # Maximum value
mode_val <- 0.8  # Mode (peak) value

# Generate the data
observed_samples <- rtriangle(n, a = min_val, b = max_val, c = mode_val)

# Create bins for the range of the distribution
num_bins <- 10  # Number of bins
breaks <- seq(min_val, max_val, length.out = num_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(observed_samples, breaks = breaks, plot = FALSE)$counts

# Function for Triangular PDF
triangular_pdf <- function(x, a, b, c) {
  if (x < a || x > b) return(0)
  if (x <= c) {
    return(2 * (x - a) / ((b - a) * (c - a)))
  } else {
    return(2 * (b - x) / ((b - a) * (b - c)))
  }
}

# Calculate expected frequencies
expected_freq <- sapply(1:num_bins, function(i) {
  lower <- breaks[i]
  upper <- breaks[i + 1]
  mid <- (lower + upper) / 2  # Approximate PDF using the bin midpoint
  prob <- triangular_pdf(mid, min_val, max_val, mode_val) * (upper - lower)
  return(prob)
})

# Scale probabilities to match the total number of observations
expected_freq <- expected_freq * n

# Perform the Chi-Square Goodness-of-Fit Test
chisq_result <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Display the result
print(chisq_result)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 0.47266, df = 9, p-value = 1

# Visualization
barplot(rbind(observed_freq, expected_freq),
        beside = TRUE,
        col = c("lightblue", "lightgreen"),
        names.arg = round(breaks[-1], 1),
        legend.text = c("Observed", "Expected"),
        args.legend = list(x = "topright"),
        main = "Observed vs Expected Frequencies (Triangular Distribution)",
        xlab = "Bins",
        ylab = "Frequency")

Normal Distrubiton: ICE Truck Operational Emissions

# Load required libraries
library(ggplot2)

# Set parameters for the normal distribution
set.seed(123)
n_samples <- 5000  # Number of samples
mean_value <- 5  # Mean of the normal distribution
sd_value <- 0.4   # Standard deviation of the normal distribution

# Generate random samples from a normal distribution
sample_data <- rnorm(n_samples, mean = mean_value, sd = sd_value)

# Define the number of bins
n_bins <- 10

# Create bins for the histogram
breaks <- seq(min(sample_data), max(sample_data), length.out = n_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(sample_data, breaks = breaks, plot = FALSE)$counts

# Calculate expected frequencies using the normal PDF
expected_freq <- diff(pnorm(breaks, mean = mean_value, sd = sd_value)) * n_samples

# Ensure no expected frequency is zero (adjust if necessary)
if (any(expected_freq == 0)) {
  stop("Some expected frequencies are zero. Adjust the number of bins or sample size.")
}

# Perform the Chi-Square Goodness-of-Fit Test
chi_sq_test <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Print the test results
print(chi_sq_test)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 5.8746, df = 9, p-value = 0.7524

# Plot the observed vs expected frequencies
barplot(
  rbind(observed_freq, expected_freq),
  beside = TRUE,
  col = c("skyblue", "orange"),
  legend.text = c("Observed", "Expected"),
  args.legend = list(x = "topright"),
  main = "Observed vs Expected Frequencies (Normal Distribution)",
  xlab = "Bins",
  ylab = "Frequency"
)

Normal Distrubiton: ICE Car Operational Emissions

# Load required libraries
library(ggplot2)

# Set parameters for the normal distribution
set.seed(123)
n_samples <- 5000  # Number of samples
mean_value <- 4.15  # Mean of the normal distribution
sd_value <- 0.675   # Standard deviation of the normal distribution

# Generate random samples from a normal distribution
sample_data <- rnorm(n_samples, mean = mean_value, sd = sd_value)

# Define the number of bins
n_bins <- 10

# Create bins for the histogram
breaks <- seq(min(sample_data), max(sample_data), length.out = n_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(sample_data, breaks = breaks, plot = FALSE)$counts

# Calculate expected frequencies using the normal PDF
expected_freq <- diff(pnorm(breaks, mean = mean_value, sd = sd_value)) * n_samples

# Ensure no expected frequency is zero (adjust if necessary)
if (any(expected_freq == 0)) {
  stop("Some expected frequencies are zero. Adjust the number of bins or sample size.")
}

# Perform the Chi-Square Goodness-of-Fit Test
chi_sq_test <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Print the test results
print(chi_sq_test)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 5.8746, df = 9, p-value = 0.7524

# Plot the observed vs expected frequencies
barplot(
  rbind(observed_freq, expected_freq),
  beside = TRUE,
  col = c("skyblue", "orange"),
  legend.text = c("Observed", "Expected"),
  args.legend = list(x = "topright"),
  main = "Observed vs Expected Frequencies (Normal Distribution)",
  xlab = "Bins",
  ylab = "Frequency"
)

Normal Distrubiton: BEV Car Energy Consumption

# Load required libraries
library(ggplot2)

# Set parameters for the normal distribution
set.seed(123)
n_samples <- 5000  # Number of samples
mean_value <- 990000  # Mean of the normal distribution
sd_value <- 156174   # Standard deviation of the normal distribution

# Generate random samples from a normal distribution
sample_data <- rnorm(n_samples, mean = mean_value, sd = sd_value)

# Define the number of bins
n_bins <- 10

# Create bins for the histogram
breaks <- seq(min(sample_data), max(sample_data), length.out = n_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(sample_data, breaks = breaks, plot = FALSE)$counts

# Calculate expected frequencies using the normal PDF
expected_freq <- diff(pnorm(breaks, mean = mean_value, sd = sd_value)) * n_samples

# Ensure no expected frequency is zero (adjust if necessary)
if (any(expected_freq == 0)) {
  stop("Some expected frequencies are zero. Adjust the number of bins or sample size.")
}

# Perform the Chi-Square Goodness-of-Fit Test
chi_sq_test <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Print the test results
print(chi_sq_test)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 5.8746, df = 9, p-value = 0.7524

# Plot the observed vs expected frequencies
barplot(
  rbind(observed_freq, expected_freq),
  beside = TRUE,
  col = c("skyblue", "orange"),
  legend.text = c("Observed", "Expected"),
  args.legend = list(x = "topright"),
  main = "Observed vs Expected Frequencies (Normal Distribution)",
  xlab = "Bins",
  ylab = "Frequency"
)

Normal Distrubiton: BEV Truck Energy Consumption

# Load required libraries
library(ggplot2)

# Set parameters for the normal distribution
set.seed(123)
n_samples <- 5000  # Number of samples
mean_value <- 1800000  # Mean of the normal distribution
sd_value <- 120000   # Standard deviation of the normal distribution

# Generate random samples from a normal distribution
sample_data <- rnorm(n_samples, mean = mean_value, sd = sd_value)

# Define the number of bins
n_bins <- 10

# Create bins for the histogram
breaks <- seq(min(sample_data), max(sample_data), length.out = n_bins + 1)

# Calculate observed frequencies
observed_freq <- hist(sample_data, breaks = breaks, plot = FALSE)$counts

# Calculate expected frequencies using the normal PDF
expected_freq <- diff(pnorm(breaks, mean = mean_value, sd = sd_value)) * n_samples

# Ensure no expected frequency is zero (adjust if necessary)
if (any(expected_freq == 0)) {
  stop("Some expected frequencies are zero. Adjust the number of bins or sample size.")
}

# Perform the Chi-Square Goodness-of-Fit Test
chi_sq_test <- chisq.test(x = observed_freq, p = expected_freq / sum(expected_freq), rescale.p = TRUE)

# Print the test results
print(chi_sq_test)

## 
##  Chi-squared test for given probabilities
## 
## data:  observed_freq
## X-squared = 5.8746, df = 9, p-value = 0.7524

# Plot the observed vs expected frequencies
barplot(
  rbind(observed_freq, expected_freq),
  beside = TRUE,
  col = c("skyblue", "orange"),
  legend.text = c("Observed", "Expected"),
  args.legend = list(x = "topright"),
  main = "Observed vs Expected Frequencies (Normal Distribution)",
  xlab = "Bins",
  ylab = "Frequency"
)

******************************************************************************

***** Sobol Indext Sensativity Testing *******************************

******************************************************************************

# ------------------------------
# Fleet-Wide Monte Carlo Sobol Sensitivity Analysis
# ------------------------------

# Install and load required libraries
if (!requireNamespace("sensitivity", quietly = TRUE)) install.packages("sensitivity")

## Registered S3 method overwritten by 'sensitivity':
##   method    from 
##   print.src dplyr

if (!requireNamespace("lhs", quietly = TRUE)) install.packages("lhs")

library(sensitivity)

## Warning: package 'sensitivity' was built under R version 4.4.2

## 
## Attaching package: 'sensitivity'

## The following object is masked from 'package:dplyr':
## 
##     src

library(lhs)

## Warning: package 'lhs' was built under R version 4.4.2

# 1. Define Parameter Ranges
param_ranges_fleet <- data.frame(
  ice_car_consumption = c(3e6, 5e6),       # Range for ICE car consumption (kWh)
  ice_truck_consumption = c(6e6, 8e6),     # Range for ICE truck consumption (kWh)
  bev_car_consumption = c(720000, 1260000),# Range for BEV car consumption (kWh)
  bev_truck_consumption = c(1.6e6, 2e6)    # Range for BEV truck consumption (kWh)
)

# 2. Generate Latin Hypercube Samples
n_samples <- 50000

lhs_samples_fleet <- randomLHS(n_samples, ncol(param_ranges_fleet))
colnames(lhs_samples_fleet) <- colnames(param_ranges_fleet)

# Scale samples to parameter ranges
lhs_scaled_fleet <- as.data.frame(
  sweep(lhs_samples_fleet, 2, param_ranges_fleet[,1], FUN = "*") +
  sweep(lhs_samples_fleet, 2, diff(t(param_ranges_fleet)), FUN = "*")
)

## Warning in sweep(lhs_samples_fleet, 2, diff(t(param_ranges_fleet)), FUN = "*"):
## STATS is longer than the extent of 'dim(x)[MARGIN]'

# 3. Define the Fleet-Wide Simulation Model
fleet_model <- function(params) {
  ice_car_energy <- params$ice_car_consumption * 0.8  # Arbitrary efficiency factor
  ice_truck_energy <- params$ice_truck_consumption * 0.85
  bev_car_energy <- params$bev_car_consumption * 0.9
  bev_truck_energy <- params$bev_truck_consumption * 0.95
  
  total_energy <- ice_car_energy + ice_truck_energy + bev_car_energy + bev_truck_energy
  return(total_energy)
}

# 4. Perform Sobol Sensitivity Analysis
sobol_result_fleet <- sobol(
  model = fleet_model,
  X1 = lhs_scaled_fleet,
  X2 = lhs_scaled_fleet,
  nboot = 100
)

# 5. Display and Visualize Results
print(sobol_result_fleet)

## 
## Call:
## sobol(model = fleet_model, X1 = lhs_scaled_fleet, X2 = lhs_scaled_fleet,     nboot = 100)
## 
## Model runs: 250000 
## 
## Sobol indices
##                       original         bias   std. error min. c.i. max. c.i.
## ice_car_consumption   1.000066 5.146731e-08 5.276134e-07  1.000065  1.000067
## ice_truck_consumption 1.000066 5.146731e-08 5.276134e-07  1.000065  1.000067
## bev_car_consumption   1.000066 5.146731e-08 5.276134e-07  1.000065  1.000067
## bev_truck_consumption 1.000066 5.146731e-08 5.276134e-07  1.000065  1.000067

plot(sobol_result_fleet, main = "Sobol Sensitivity Analysis: Fleet-Wide Monte Carlo Simulation")

# ------------------------------
# Single-Vehicle Monte Carlo Sobol Sensitivity Analysis
# ------------------------------

# Install and load required libraries
if (!requireNamespace("sensitivity", quietly = TRUE)) install.packages("sensitivity")
if (!requireNamespace("lhs", quietly = TRUE)) install.packages("lhs")

library(sensitivity)
library(lhs)

# 1. Define Parameter Ranges
param_ranges_single <- data.frame(
  annual_miles_singular = c(4785, 18858),  # Range for annual miles driven
  bev_battery_size = c(40, 135),          # Range for BEV battery size (kWh)
  car_co2_per_10k_miles = c(2.6, 5.2)     # Range for ICE CO2 emissions per 10,000 miles (tons)
)

# 2. Generate Latin Hypercube Samples
n_samples <- 1000

lhs_samples_single <- randomLHS(n_samples, ncol(param_ranges_single))
colnames(lhs_samples_single) <- colnames(param_ranges_single)

# Scale samples to parameter ranges
lhs_scaled_single <- as.data.frame(
  sweep(lhs_samples_single, 2, param_ranges_single[,1], FUN = "*") +
  sweep(lhs_samples_single, 2, diff(t(param_ranges_single)), FUN = "*")
)

## Warning in sweep(lhs_samples_single, 2, param_ranges_single[, 1], FUN = "*"):
## STATS does not recycle exactly across MARGIN

## Warning in sweep(lhs_samples_single, 2, diff(t(param_ranges_single)), FUN =
## "*"): STATS is longer than the extent of 'dim(x)[MARGIN]'

# 3. Define the Single-Vehicle Simulation Model
single_vehicle_model <- function(params) {
  bev_battery_co2 <- params$bev_battery_size * 0.14  # Upstream CO2 emissions per kWh
  ice_emissions <- params$car_co2_per_10k_miles * (params$annual_miles_singular / 10000)
  
  total_emissions <- bev_battery_co2 + ice_emissions
  return(total_emissions)
}

# 4. Perform Sobol Sensitivity Analysis
sobol_result_single <- sobol(
  model = single_vehicle_model,
  X1 = lhs_scaled_single,
  X2 = lhs_scaled_single,
  nboot = 100
)

# 5. Display and Visualize Results
print(sobol_result_single)

## 
## Call:
## sobol(model = single_vehicle_model, X1 = lhs_scaled_single, X2 = lhs_scaled_single,     nboot = 100)
## 
## Model runs: 4000 
## 
## Sobol indices
##                       original         bias   std. error min. c.i. max. c.i.
## annual_miles_singular 1.000545 2.913504e-06 2.300299e-05  1.000497  1.000585
## bev_battery_size      1.000545 2.913504e-06 2.300299e-05  1.000497  1.000585
## car_co2_per_10k_miles 1.000545 2.913504e-06 2.300299e-05  1.000497  1.000585

plot(sobol_result_single, main = "Sobol Sensitivity Analysis: Single-Vehicle Monte Carlo Simulation")

******************************************************************************

***** Wilcoxon Rank-Sum Test Sensativity Testing **********************

******************************************************************************

# Fleet-Wide Monte Carlo Wilcoxon Rank-Sum Test

# Ensure necessary libraries are loaded
library(ggplot2)

# Variables from Fleet-Wide Simulation
fleet_emissions_run1 <- rowSums(bev_car_co2_matrix)  # Total BEV car emissions per iteration
fleet_emissions_run2 <- rowSums(ice_car_co2_matrix)  # Total ICE car emissions per iteration

# Remove any NA values
fleet_emissions_run1 <- na.omit(fleet_emissions_run1)
fleet_emissions_run2 <- na.omit(fleet_emissions_run2)

# Perform Wilcoxon Rank-Sum Test
wilcox_fleet <- wilcox.test(fleet_emissions_run1, fleet_emissions_run2, 
                            alternative = "two.sided", 
                            exact = FALSE, 
                            correct = TRUE)

# Print test results
print(wilcox_fleet)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  fleet_emissions_run1 and fleet_emissions_run2
## W = 3600, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

# Visualization for Comparison
boxplot(fleet_emissions_run1, fleet_emissions_run2,
        names = c("BEV Fleet Emissions", "ICE Fleet Emissions"),
        col = c("skyblue", "lightgreen"),
        main = "Wilcoxon Test: Fleet-Wide Simulation Comparison",
        ylab = "CO₂ Emissions (Tons)")

# Interpretation:
if (wilcox_fleet$p.value < 0.05) {
  print("Significant difference detected between BEV and ICE Fleet Simulation Outputs.")
} else {
  print("No significant difference detected between BEV and ICE Fleet Simulation Outputs.")
}

## [1] "Significant difference detected between BEV and ICE Fleet Simulation Outputs."

Wilcoxon Test for the Single Vehicle Monte Carlo Method

# Single-Vehicle Monte Carlo Wilcoxon Rank-Sum Test


# Variables from Single-Vehicle Simulation
single_vehicle_emissions_run1 <- mean_bev_car_co2  # BEV single vehicle emissions
single_vehicle_emissions_run2 <- mean_ice_car_co2  # ICE single vehicle emissions

# Remove any NA values
single_vehicle_emissions_run1 <- na.omit(single_vehicle_emissions_run1)
single_vehicle_emissions_run2 <- na.omit(single_vehicle_emissions_run2)

# Perform Wilcoxon Rank-Sum Test
wilcox_single <- wilcox.test(single_vehicle_emissions_run1, single_vehicle_emissions_run2, 
                             alternative = "two.sided", 
                             exact = FALSE, 
                             correct = TRUE)

# Print test results
print(wilcox_single)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  single_vehicle_emissions_run1 and single_vehicle_emissions_run2
## W = 2524, p-value = 0.0001462
## alternative hypothesis: true location shift is not equal to 0

# Visualization for Comparison
boxplot(single_vehicle_emissions_run1, single_vehicle_emissions_run2,
        names = c("BEV Vehicle Emissions", "ICE Vehicle Emissions"),
        col = c("lightcoral", "lightblue"),
        main = "Wilcoxon Test: Single-Vehicle Simulation Comparison",
        ylab = "CO₂ Emissions (Tons)")

# Interpretation:
if (wilcox_single$p.value < 0.05) {
  print("Significant difference detected between BEV and ICE Single Vehicle Simulation Outputs.")
} else {
  print("No significant difference detected between BEV and ICE Single Vehicle Simulation Outputs.")
}

## [1] "Significant difference detected between BEV and ICE Single Vehicle Simulation Outputs."

******************************************************************************

***** Regressional Analysis for RSM in R *****************************

******************************************************************************

# Load necessary libraries
library(ggplot2)
library(dplyr)
library(car)  # For ANOVA analysis

## Warning: package 'car' was built under R version 4.4.2

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.4.2

## 
## Attaching package: 'car'

## The following object is masked from 'package:sensitivity':
## 
##     scatterplot

## The following object is masked from 'package:dplyr':
## 
##     recode

# Expand grid for all three variables
battery_ice_grid <- expand.grid(
  Battery = bev_battery_size,
  ICE_CO2 = ice_co2_output,
  co2_kwh = co2_per_kwh
)

# Calculate Miles for each combination
battery_ice_grid$Miles <- with(battery_ice_grid, miles_to_offset(Battery, ICE_CO2, co2_kwh))


# Fit the Response Surface Regression Model
rsm_model <- lm(Miles ~ Battery + ICE_CO2 + co2_kwh, data = battery_ice_grid)

# Display summary of the regression model
summary_rsm <- summary(rsm_model)
print(summary_rsm)

## 
## Call:
## lm(formula = Miles ~ Battery + ICE_CO2 + co2_kwh, data = battery_ice_grid)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -29762  -9706  -4932   5382 147897 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13776.713    830.727   16.58   <2e-16 ***
## Battery        341.806      3.445   99.22   <2e-16 ***
## ICE_CO2     -11256.302    100.426 -112.08   <2e-16 ***
## co2_kwh     344540.909   3694.147   93.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16780 on 8996 degrees of freedom
## Multiple R-squared:  0.7757, Adjusted R-squared:  0.7756 
## F-statistic: 1.037e+04 on 3 and 8996 DF,  p-value: < 2.2e-16