# Monte Carlo Simulation

## Problem 1 - Static Monte Carlo Project Cost Simulation

Problem Statement:

Consider a Project that is characterized by the following activity edge table. The table indicates: 1) the activity names; 2) the associated event node names (numbers); 3) the duration times (in days) for each activity; and 4) the hourly labor costs for each activity. Assume two people are required to complete each task and that each works 8 hours a day.

 Activity Nodes Duration Est (days) Labor Cost (\$/h) Min Max Likely Min Max Likely A 1, 2 4 8 6 \$100 \$150 \$120 B 1, 3 2 8 4 \$100 \$150 \$120 C 2, 4 1 7 3 \$100 \$150 \$120 D 3, 4 6 12 9 \$100 \$150 \$120 F 3, 5 5 15 10 \$100 \$150 \$120 I 3, 6 10 25 15 \$100 \$150 \$120 J 4, 7 5 12 9 \$90 \$110 \$100 K 5, 7 1 3 2 \$90 \$110 \$100 L 6, 8 2 6 3 \$90 \$110 \$100 M 7, 9 10 20 15 \$90 \$110 \$100 N 8, 9 6 11 9 \$90 \$110 \$100

(Assume triangular distribution)

Use Monte Carlo simulation to calculate:

• mean total cost of the project
• the associated standard deviation
• the associated standard error in the mean
• the cumulative running mean (CRM)

MATLAB Code: problem_1.m

Author: Yash Bansod
Date: 27th February, 2020
Problem 1 - Static Monte Carlo Project Cost Simulation

GitHub: https://github.com/YashBansod

### Clear the environment and the command line

``````clear;
clc;
close all;
``````

### Define the simulation parmeters

``````num_samples = 2000;
num_hist_bins = 51;

activity =  {'A' 'B' 'C' 'D' 'F' 'I' 'J' 'K' 'L' 'M' 'N'};
num_activities = size(activity, 2);

min_hours = [64 32 16 96 80 160 80 16 32 160 96];
max_hours = [128 128 112 192 240 400 192 48 96 320 176];
likely_hours = [96 64 48 144 160 240 144 32 48 240 144];

min_cost = [100 100 100 100 100 100 90 90 90 90 90];
max_cost = [150 150 150 150 150 150 110 110 110 110 110];
likely_cost = [120 120 120 120 120 120 100 100 100 100 100];

% Some dimension checks to make sure values were inputted correctly
assert(size(min_hours, 2) == num_activities);
assert(size(max_hours, 2) == num_activities);
assert(size(likely_hours, 2) == num_activities);

assert(size(min_cost, 2) == num_activities);
assert(size(max_cost, 2) == num_activities);
assert(size(likely_cost, 2) == num_activities);
``````

### Define the probability distributions to sample data from

``````% Create probability distribution samplers for labor hours and cost
prob_dist_hours = cell(size(likely_hours));
prob_dist_cost = cell(size(likely_cost));

for ind = 1:num_activities
prob_dist_hours{ind} = makedist('triangular', ...
'a', min_hours(ind), 'b', likely_hours(ind), 'c', max_hours(ind));
prob_dist_cost{ind} = makedist('triangular', ...
'a', min_cost(ind), 'b', likely_cost(ind), 'c', max_cost(ind));
end
``````

### Run the Monte-Carlo Simulation

``````% rand('seed', 5);                    % Seed to control randomization
tic;

rand_hours = zeros(num_activities, num_samples);
rand_cost = zeros(num_activities, num_samples);

for ind = 1:num_activities
rand_hours(ind, :) = random(prob_dist_hours{ind}, 1, num_samples);
rand_cost(ind, :) = random(prob_dist_cost{ind}, 1, num_samples);
end

sample_range = 1:num_samples;
cum_run_mean_proj_cost = cumsum(rand_proj_cost) ./ sample_range;
mean_proj_cost = cum_run_mean_proj_cost(num_samples);

toc;

std_dev_proj_cost = std(rand_proj_cost);
std_err_proj_cost = std_dev_proj_cost / sqrt(num_samples);

% Calculate 10%ile and 90%ile of the collected dataset
percent_10 = prctile(rand_proj_cost, 10);
percent_90 = prctile(rand_proj_cost, 90);

% 10%ile & 90%ile assuming data is distributed as normal distribution.
norm_percent_10 = norminv(0.1, mean_proj_cost, std_dev_proj_cost);
norm_percent_90 = norminv(0.9, mean_proj_cost, std_dev_proj_cost);
Elapsed time is 0.004723 seconds.
``````

### Plot the results

``````figure(1)
histogram(rand_proj_cost, num_hist_bins);
hold on
xline(mean_proj_cost, '--r');
xline(percent_10, '-.g');
xline(norm_percent_10, '-.c');
xline(percent_90, '-.g');
xline(norm_percent_90, '-.c');
hold off
legend('Histogram', 'Final Mean', ...
'10%->90% Actual Boundary', '10%->90% NormInv Boundary');
title('Histogram: Project Cost Distribution');
xlabel('Project Cost in \$');
ylabel('Count');

figure(2)
plot(sample_range, cum_run_mean_proj_cost);
hold on
yline(mean_proj_cost, '--r');
hold off
legend('CRM', 'Final Mean');
title('Cumulative Running Mean (CRM) of Project Cost');
xlabel('Num Samples');
ylabel('CRM of Project Cost in \$');
``````
``````fprintf('Project Cost - Mean: %.2f\n', mean_proj_cost);
fprintf('Project Cost - Standard Deviation: %.2f\n', std_dev_proj_cost)
fprintf('Project Cost - Standard Error: %.2f\n', std_err_proj_cost)
fprintf('Project Cost - 10%% -> 90%% Boundary: (%.2f, %.2f)\n',...
percent_10, percent_90);
fprintf('Project Cost - 10%% -> 90%% Norm Inv Boundary: (%.2f, %.2f)\n',...
norm_percent_10, norm_percent_90);
``````

Project Cost - Mean: 158972.94
Project Cost - Standard Deviation: 10503.04
Project Cost - Standard Error: 234.86
Project Cost - 10% -> 90% Boundary: (145860.56, 172652.69)
Project Cost - 10% -> 90% Norm Inv Boundary: (145512.75, 172433.14)

## Problem 2a - Time-Driven Dynamic MC Simulation of an AV

Consider an autonomous vehicle (AV) attempting to move in the “X direction” subject to random gusts of wind (blowing in the +/- Y direction). Assume the AV has a mass of 2 kg and is traveling at a speed of 60 km/h. Assume that the impulse (change in momentum) imparted to the AV by the wind gusts over a period of 5 seconds is characterized by a normal distribution with a mean of 0 kg km/h and a standard deviation of 10 kg km/hr.

a. Develop and submit a MATLAB program that:

1. Accepts the speed and mass of the AV, the time step size, the mean and standard deviation of the impulse delivered to the AV, the duration of the simulation (Ts), and number of bins. All units should be converted to MKS.
2. Provides the following output: i. The final location of the AV (at t = Ts). ii. A 2-D plot of the trajectory of AV (at dt time intervals).

MATLAB Code: problem_2a.m

Author: Yash Bansod
Date: 27th February, 2020
Problem 2a - Time-Driven Dynamic MC Simulation of an AV

GitHub: https://github.com/YashBansod

### Clear the environment and the command line

``````clear;
clc;
close all;
``````

### Define the simulation parmeters

``````num_samples = 1;

sim_time_len = 1000;
time_step = 5;
num_t_steps = floor(sim_time_len/time_step);

mass_av = 2;
init_v_x = 16.6666;             % Initial velocity of AV in x-axis
init_v_y = 0;                   % Initial velocity of AV in y-axis

init_s_x = 0;                   % Initial position of AV in x-axis
init_s_y = 0;                   % Initial position of AV in y-axis

delta_p_y_mean = 0;             % Mean of Change in momentum
delta_p_std = 2.7777;           % Standard deviation of change in momentum
``````

### Define the probability distributions to sample data from

``````delta_v_y_mean = delta_p_y_mean/mass_av;
delta_v_y_std = delta_p_std/mass_av;

delta_v_y_dist = makedist('Normal', 'mu', ...
delta_v_y_mean, 'sigma', delta_v_y_std);
``````

### Run the Monte-Carlo Simulation

``````% x-axis velocity trajectory of the AV
v_x = zeros(num_samples, num_t_steps + 1);
v_x(:, :) = init_v_x;

% x-axis position trajectory of the AV
delta_s_x = v_x(:, 1:end - 1) .* time_step;
s_x = zeros(num_samples, num_t_steps + 1);
s_x(:, 1) = init_s_x;
s_x(:, 2:end) = delta_s_x;
s_x = cumsum(s_x, 2);

tic;

% y-axis velocity trajectory of the AV
delta_v_y = random(delta_v_y_dist, num_samples, num_t_steps);
v_y = zeros(num_samples, num_t_steps + 1);
v_y(:, 1) = init_v_y;
v_y(:, 2:end) = delta_v_y;
v_y = cumsum(v_y, 2);

% y-axis positioin trajectory of the AB
delta_s_y = v_y(:, 1:end - 1) .* time_step + 0.5 * delta_v_y .* time_step;
s_y = zeros(num_samples, num_t_steps + 1);
s_y(:, 1) = init_s_y;
s_y(:, 2:end) = delta_s_y;
s_y = cumsum(s_y, 2);

end_y_cord = s_y(:, end);

toc;
Elapsed time is 0.000832 seconds.
``````

### Plot the results

``````figure(1)
plot(0:time_step:sim_time_len, s_y);
grid on;
title('Y-coordinate Trajectory of AV');
xlabel('Time (seconds)');
ylabel('Meters');

figure(2)
plot(s_x, s_y);
grid on;
title('Position Trajectory of AV');
xlabel('Meters');
ylabel('Meters');
``````
``````fprintf('Final Position Coordinate: (%.2f, %.2f)\n', s_x(end), s_y(end));
fprintf('Final Velocity: (%.2f, %.2f)\n', v_x(end), v_y(end));
``````

Final Position Coordinate: (16666.60, 2914.33) Final Velocity: (16.67, 10.74)

## Problem 2b - Time-Driven Dynamic MC Simulation of an AV

b. Develop and submit a MATLAB program that:

1. Accepts the speed and mass of the AV, the time step size, the mean and standard deviation of the impulse delivered to the AV, the duration of the simulation (Ts), the number of samples (Ns), and number of bins (Nb). All units should be converted to MKS.

2. Provides the following output:

i. The mean ending Y coordinate (in m), its associated standard deviation, and standard error in the mean. ii. A histogram graph showing the distribution of ending Y locations.

MATLAB Code: problem_2b.m

Author: Yash Bansod
Date: 27th February, 2020
Problem 2b - Time-Driven Dynamic MC Simulation of an AV

GitHub: https://github.com/YashBansod

### Clear the environment and the command line

``````clear;
clc;
close all;
``````

### Define the simulation parameters

``````num_samples = 5000;
num_hist_bins = 21;

sim_time_len = 400;
time_step = 5;
num_t_steps = floor(sim_time_len/time_step);

mass_av = 2;
init_v_x = 16.6666;             % Initial velocity of AV in x-axis
init_v_y = 0;                   % Initial velocity of AV in y-axis

init_s_x = 0;                   % Initial position of AV in x-axis
init_s_y = 0;                   % Initial position of AV in y-axis

delta_p_y_mean = 0;             % Mean of Change in momentum
delta_p_std = 2.7777;           % Standard deviation of change in momentum
``````

### Define the probability distributions to sample data from

``````delta_v_y_mean = delta_p_y_mean/mass_av;
delta_v_y_std = delta_p_std/mass_av;

delta_v_y_dist = makedist('Normal', 'mu', ...
delta_v_y_mean, 'sigma', delta_v_y_std);
``````

### Run the Monte-Carlo Simulation

``````% x-axis velocity trajectory of the AV
v_x = zeros(num_samples, num_t_steps + 1);
v_x(:, :) = init_v_x;

% x-axis position trajectory of the AV
delta_s_x = v_x(:, 1:end - 1) .* time_step;
s_x = zeros(num_samples, num_t_steps + 1);
s_x(:, 1) = init_s_x;
s_x(:, 2:end) = delta_s_x;
s_x = cumsum(s_x, 2);

tic;

% y-axis velocity trajectory of the AV
delta_v_y = random(delta_v_y_dist, num_samples, num_t_steps);
v_y = zeros(num_samples, num_t_steps + 1);
v_y(:, 1) = init_v_y;
v_y(:, 2:end) = delta_v_y;
v_y = cumsum(v_y, 2);

% y-axis positioin trajectory of the AB
delta_s_y = v_y(:, 1:end - 1) .* time_step + 0.5 * delta_v_y .* time_step;
s_y = zeros(num_samples, num_t_steps + 1);
s_y(:, 1) = init_s_y;
s_y(:, 2:end) = delta_s_y;
s_y = cumsum(s_y, 2);

end_y_cord = s_y(:, end);
mean_end_y_cord = mean(end_y_cord);
std_dev_end_y_cord = std(end_y_cord);
std_err_end_y_cord = std_dev_end_y_cord / sqrt(num_samples);

toc;
``````

Elapsed time is 0.017619 seconds.

### Plot the results

``````figure(1)
histogram(end_y_cord, num_hist_bins);
hold on
xline(mean_end_y_cord, '--r');
hold off
legend('Histogram', 'Mean');
title('Histogram: End Y Coordinate');
xlabel('Y Coordinate Position');
ylabel('Count');
``````
``````fprintf('Final Y-Coordinate - Mean: %.2f\n', mean_end_y_cord)
fprintf('Final Y-Coordinate - Standard Deviation: %.2f\n', ...
std_dev_end_y_cord)
fprintf('Final Y-Coordinate - Standard Error: %.2f\n', std_err_end_y_cord)
``````

Final Y-Coordinate - Mean: 12.63
Final Y-Coordinate - Standard Deviation: 2855.95
Final Y-Coordinate - Standard Error: 40.39