I am computing the probablility of collision with the Monte Carlo method in MATLAB. The following is my code for it. The results don't seem to be accurate. Is there something wrong with my code? Also, it takes more than 5 min to run this code. Can you tell me why? %Initial Mean and Covariance mu_01 = [153446.180; 41874155.872;0;3066.875;-11.374;0]/1000; P_01 = [6494.080 -376.139 0 0.0160 -0.494 0; -376.139 22.560 0 -9.883e-4 0.0286 0; 0 0 1.205 0 0 -6.071e-5; 0.0160 -9.883e-4 0 4.437e-8 -1.212e-6 0; -0.494 0.0286 0 -1.212e-6 3.762e-5 0; 0 0 -6.071e-5 0 0 3.390e-9]*1e-6; mu_02 = [153446.679; 41874156.372;5.000;3066.865;-11.374;-1.358e-6]/1000; P_02 = [6494.224 -376.156 -4.492e-5 0.0160 -0.494 -5.902e-8; -376.156 22.561 2.550e-6 -9.885e-3 0.0286 3.419e-9; -4.492e-5 2.550e-6 1.205 -1.180e-10 3.419e-9 -6.072e-5; 0.0160 -9.883e-3 -1.180e-10 4.438e-8 -1.212e-6 -1.448e-13; -0.494 0.0286 3.419e-9 -1.212e-6 3.762e-5 4.492e-12; -5.902e-8 3.419e-9 -6.072e-5 -1.448e-13 4.492e-12 3.392e-9]*1e-6; %Collision parameters R = 6378.0; %km mu = 398600.4415; %km^3/s^2 a1 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period1 = 2*pi*sqrt((a1^3)/mu); t_total1 = 0:10:2*period1; a2 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period2 = 2*pi*sqrt((a2^3)/mu); t_total2 = 0:10:2*period2; rho_AB = 15e-3; function [PC, n_collide] = MonteCarlo_PC(tm1,P_01,mu_01,tm2,P_02,mu_02,rho_AB) %Creating Eigen Matrix and Eigen Vector [R1,Lam1] = eig(P_01); [R2,Lam2] = eig(P_02); %Creating state vector x_01 = inv(R1)*(R1*mu_01 + sqrt(abs(diag(Lam1)))*randn(1,1000)); x_02 = inv(R2)*(R2*mu_02 + sqrt(abs(diag(Lam2)))*randn(1,1000)); x_final1 = zeros(6,1000); x_final2 = zeros(6,1000); % Propogate Orbit for i=1:1000 % Integrating the state using ode45 integrator init_x = x_01(:,i); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm1, init_x, options); x_final1(:,i) = x_out(length(tm1),:); end for j=1:1000 % Integrating the state using ode45 integrator init_x = x_02(:,j); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm2, init_x, options); x_final2(:,j) = x_out(length(tm2),:); end n_sampleA = 1000; n_sampleB = 1000; n_collide = 0; for idx1=1:n_sampleA for idx2=1:n_sampleB x1 = x_final1(1:3,idx1); x2 = x_final2(1:3,idx2); dist = sqrt((x1(1)-x2(1))^2 + (x1(2)-x2(2))^2 + (x1(3)-x2(3))^2); if dist <= rho_AB n_collide = n_collide + 1; end end end PC = n_collide/(n_sampleA*n_sampleB); % Plot fig1 = figure(); set(fig1, 'color', 'white'); scatter3(x_final1(1,:),x_final1(2,:),x_final1(3,:)); fig2 = figure(); set(fig2, 'color', 'white'); scatter3(x_final2(1,:),x_final2(2,:),x_final2(3,:)); end
I am computing the probablility of collision with the Monte Carlo method in MATLAB. The following is my code for it. The results don't seem to be accurate. Is there something wrong with my code? Also, it takes more than 5 min to run this code. Can you tell me why? %Initial Mean and Covariance mu_01 = [153446.180; 41874155.872;0;3066.875;-11.374;0]/1000; P_01 = [6494.080 -376.139 0 0.0160 -0.494 0; -376.139 22.560 0 -9.883e-4 0.0286 0; 0 0 1.205 0 0 -6.071e-5; 0.0160 -9.883e-4 0 4.437e-8 -1.212e-6 0; -0.494 0.0286 0 -1.212e-6 3.762e-5 0; 0 0 -6.071e-5 0 0 3.390e-9]*1e-6; mu_02 = [153446.679; 41874156.372;5.000;3066.865;-11.374;-1.358e-6]/1000; P_02 = [6494.224 -376.156 -4.492e-5 0.0160 -0.494 -5.902e-8; -376.156 22.561 2.550e-6 -9.885e-3 0.0286 3.419e-9; -4.492e-5 2.550e-6 1.205 -1.180e-10 3.419e-9 -6.072e-5; 0.0160 -9.883e-3 -1.180e-10 4.438e-8 -1.212e-6 -1.448e-13; -0.494 0.0286 3.419e-9 -1.212e-6 3.762e-5 4.492e-12; -5.902e-8 3.419e-9 -6.072e-5 -1.448e-13 4.492e-12 3.392e-9]*1e-6; %Collision parameters R = 6378.0; %km mu = 398600.4415; %km^3/s^2 a1 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period1 = 2*pi*sqrt((a1^3)/mu); t_total1 = 0:10:2*period1; a2 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period2 = 2*pi*sqrt((a2^3)/mu); t_total2 = 0:10:2*period2; rho_AB = 15e-3; function [PC, n_collide] = MonteCarlo_PC(tm1,P_01,mu_01,tm2,P_02,mu_02,rho_AB) %Creating Eigen Matrix and Eigen Vector [R1,Lam1] = eig(P_01); [R2,Lam2] = eig(P_02); %Creating state vector x_01 = inv(R1)*(R1*mu_01 + sqrt(abs(diag(Lam1)))*randn(1,1000)); x_02 = inv(R2)*(R2*mu_02 + sqrt(abs(diag(Lam2)))*randn(1,1000)); x_final1 = zeros(6,1000); x_final2 = zeros(6,1000); % Propogate Orbit for i=1:1000 % Integrating the state using ode45 integrator init_x = x_01(:,i); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm1, init_x, options); x_final1(:,i) = x_out(length(tm1),:); end for j=1:1000 % Integrating the state using ode45 integrator init_x = x_02(:,j); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm2, init_x, options); x_final2(:,j) = x_out(length(tm2),:); end n_sampleA = 1000; n_sampleB = 1000; n_collide = 0; for idx1=1:n_sampleA for idx2=1:n_sampleB x1 = x_final1(1:3,idx1); x2 = x_final2(1:3,idx2); dist = sqrt((x1(1)-x2(1))^2 + (x1(2)-x2(2))^2 + (x1(3)-x2(3))^2); if dist <= rho_AB n_collide = n_collide + 1; end end end PC = n_collide/(n_sampleA*n_sampleB); % Plot fig1 = figure(); set(fig1, 'color', 'white'); scatter3(x_final1(1,:),x_final1(2,:),x_final1(3,:)); fig2 = figure(); set(fig2, 'color', 'white'); scatter3(x_final2(1,:),x_final2(2,:),x_final2(3,:)); end
Elements Of Electromagnetics
7th Edition
ISBN:9780190698614
Author:Sadiku, Matthew N. O.
Publisher:Sadiku, Matthew N. O.
ChapterMA: Math Assessment
Section: Chapter Questions
Problem 1.1MA
Related questions
Question
I am computing the probablility of collision with the Monte Carlo method in MATLAB. The following is my code for it. The results don't seem to be accurate. Is there something wrong with my code? Also, it takes more than 5 min to run this code. Can you tell me why?
%Initial Mean and Covariance
mu_01 = [153446.180; 41874155.872;0;3066.875;-11.374;0]/1000;
P_01 = [6494.080 -376.139 0 0.0160 -0.494 0;
-376.139 22.560 0 -9.883e-4 0.0286 0;
0 0 1.205 0 0 -6.071e-5;
0.0160 -9.883e-4 0 4.437e-8 -1.212e-6 0;
-0.494 0.0286 0 -1.212e-6 3.762e-5 0;
0 0 -6.071e-5 0 0 3.390e-9]*1e-6;
mu_02 = [153446.679; 41874156.372;5.000;3066.865;-11.374;-1.358e-6]/1000;
P_02 = [6494.224 -376.156 -4.492e-5 0.0160 -0.494 -5.902e-8;
-376.156 22.561 2.550e-6 -9.885e-3 0.0286 3.419e-9;
-4.492e-5 2.550e-6 1.205 -1.180e-10 3.419e-9 -6.072e-5;
0.0160 -9.883e-3 -1.180e-10 4.438e-8 -1.212e-6 -1.448e-13;
-0.494 0.0286 3.419e-9 -1.212e-6 3.762e-5 4.492e-12;
-5.902e-8 3.419e-9 -6.072e-5 -1.448e-13 4.492e-12 3.392e-9]*1e-6;
%Collision parameters
R = 6378.0; %km
mu = 398600.4415; %km^3/s^2
a1 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu));
period1 = 2*pi*sqrt((a1^3)/mu);
t_total1 = 0:10:2*period1;
a2 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu));
period2 = 2*pi*sqrt((a2^3)/mu);
t_total2 = 0:10:2*period2;
rho_AB = 15e-3;
function [PC, n_collide] = MonteCarlo_PC(tm1,P_01,mu_01,tm2,P_02,mu_02,rho_AB)
%Creating Eigen Matrix and Eigen Vector
[R1,Lam1] = eig(P_01);
[R2,Lam2] = eig(P_02);
%Creating state vector
x_01 = inv(R1)*(R1*mu_01 + sqrt(abs(diag(Lam1)))*randn(1,1000));
x_02 = inv(R2)*(R2*mu_02 + sqrt(abs(diag(Lam2)))*randn(1,1000));
x_final1 = zeros(6,1000);
x_final2 = zeros(6,1000);
% Propogate Orbit
for i=1:1000
% Integrating the state using ode45 integrator
init_x = x_01(:,i);
options = odeset('RelTol',1e-12, 'AbsTol',1e-12);
[~, x_out] = ode45(@TwoBody, tm1, init_x, options);
x_final1(:,i) = x_out(length(tm1),:);
end
for j=1:1000
% Integrating the state using ode45 integrator
init_x = x_02(:,j);
options = odeset('RelTol',1e-12, 'AbsTol',1e-12);
[~, x_out] = ode45(@TwoBody, tm2, init_x, options);
x_final2(:,j) = x_out(length(tm2),:);
end
n_sampleA = 1000;
n_sampleB = 1000;
n_collide = 0;
for idx1=1:n_sampleA
for idx2=1:n_sampleB
x1 = x_final1(1:3,idx1);
x2 = x_final2(1:3,idx2);
dist = sqrt((x1(1)-x2(1))^2 + (x1(2)-x2(2))^2 + (x1(3)-x2(3))^2);
if dist <= rho_AB
n_collide = n_collide + 1;
end
end
end
PC = n_collide/(n_sampleA*n_sampleB);
% Plot
fig1 = figure();
set(fig1, 'color', 'white');
scatter3(x_final1(1,:),x_final1(2,:),x_final1(3,:));
fig2 = figure();
set(fig2, 'color', 'white');
scatter3(x_final2(1,:),x_final2(2,:),x_final2(3,:));
end
Expert Solution
This question has been solved!
Explore an expertly crafted, step-by-step solution for a thorough understanding of key concepts.
Step by step
Solved in 2 steps
Recommended textbooks for you
Elements Of Electromagnetics
Mechanical Engineering
ISBN:
9780190698614
Author:
Sadiku, Matthew N. O.
Publisher:
Oxford University Press
Mechanics of Materials (10th Edition)
Mechanical Engineering
ISBN:
9780134319650
Author:
Russell C. Hibbeler
Publisher:
PEARSON
Thermodynamics: An Engineering Approach
Mechanical Engineering
ISBN:
9781259822674
Author:
Yunus A. Cengel Dr., Michael A. Boles
Publisher:
McGraw-Hill Education
Elements Of Electromagnetics
Mechanical Engineering
ISBN:
9780190698614
Author:
Sadiku, Matthew N. O.
Publisher:
Oxford University Press
Mechanics of Materials (10th Edition)
Mechanical Engineering
ISBN:
9780134319650
Author:
Russell C. Hibbeler
Publisher:
PEARSON
Thermodynamics: An Engineering Approach
Mechanical Engineering
ISBN:
9781259822674
Author:
Yunus A. Cengel Dr., Michael A. Boles
Publisher:
McGraw-Hill Education
Control Systems Engineering
Mechanical Engineering
ISBN:
9781118170519
Author:
Norman S. Nise
Publisher:
WILEY
Mechanics of Materials (MindTap Course List)
Mechanical Engineering
ISBN:
9781337093347
Author:
Barry J. Goodno, James M. Gere
Publisher:
Cengage Learning
Engineering Mechanics: Statics
Mechanical Engineering
ISBN:
9781118807330
Author:
James L. Meriam, L. G. Kraige, J. N. Bolton
Publisher:
WILEY