I am trying to compute the probability of collision with the monte carlo method in MATLAB. But I get an error in teh following code. The error is the matrix P_01 must be semi definete and P_02 must be symmetric. How should I fix this problem? % 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 rho_AB = 15e-3; % Collision threshold (in km) % Time span setup a1 = 1/((2/norm(mu_01(1:3))) - ((norm(mu_01(4:6))^2) / mu)); period1 = 2 * pi * sqrt((a1^3) / mu); tm1 = 0:10:2 * period1; a2 = 1/((2/norm(mu_02(1:3))) - ((norm(mu_02(4:6))^2) / mu)); period2 = 2 * pi * sqrt((a2^3) / mu); tm2 = 0:10:2 * period2; % Monte Carlo Simulation [PC, n_collide] = MonteCarlo_PC(tm1, P_01, mu_01, tm2, P_02, mu_02, rho_AB); function [PC, n_collide] = MonteCarlo_PC(tm1, P_01, mu_01, tm2, P_02, mu_02, rho_AB) n_samples = 1000; % Generate samples using multivariate normal distribution x_01 = mvnrnd(mu_01, P_01, n_samples)'; x_02 = mvnrnd(mu_02, P_02, n_samples)'; x_final1 = zeros(6, n_samples); x_final2 = zeros(6, n_samples); % Parallelized propagation of orbits parfor i = 1:n_samples options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12); [~, x_out1] = ode45(@TwoBody, tm1, x_01(:, i), options); x_final1(:, i) = x_out1(end, :); [~, x_out2] = ode45(@TwoBody, tm2, x_02(:, i), options); x_final2(:, i) = x_out2(end, :); end % Collision detection using vectorized operations distances = sqrt(sum((x_final1(1:3, :) - x_final2(1:3, :)).^2, 1)); n_collide = sum(distances <= rho_AB); % Probability of collision PC = n_collide / n_samples; % Plot results figure(); scatter3(x_final1(1, :), x_final1(2, :), x_final1(3, :), 'b'); title('Sample 1 Propagated Positions'); xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)'); figure(); scatter3(x_final2(1, :), x_final2(2, :), x_final2(3, :), 'r'); title('Sample 2 Propagated Positions'); xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)'); end function dxdt = TwoBody(~, x) mu = 398600.4415; % km^3/s^2 (Gravitational parameter) r = sqrt(x(1)^2 + x(2)^2 + x(3)^2); dxdt = [x(4); x(5); x(6); -mu * x(1) / r^3; -mu * x(2) / r^3; -mu * x(3) / r^3]; end
I am trying to compute the probability of collision with the monte carlo method in MATLAB. But I get an error in teh following code. The error is the matrix P_01 must be semi definete and P_02 must be symmetric. How should I fix this problem? % 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 rho_AB = 15e-3; % Collision threshold (in km) % Time span setup a1 = 1/((2/norm(mu_01(1:3))) - ((norm(mu_01(4:6))^2) / mu)); period1 = 2 * pi * sqrt((a1^3) / mu); tm1 = 0:10:2 * period1; a2 = 1/((2/norm(mu_02(1:3))) - ((norm(mu_02(4:6))^2) / mu)); period2 = 2 * pi * sqrt((a2^3) / mu); tm2 = 0:10:2 * period2; % Monte Carlo Simulation [PC, n_collide] = MonteCarlo_PC(tm1, P_01, mu_01, tm2, P_02, mu_02, rho_AB); function [PC, n_collide] = MonteCarlo_PC(tm1, P_01, mu_01, tm2, P_02, mu_02, rho_AB) n_samples = 1000; % Generate samples using multivariate normal distribution x_01 = mvnrnd(mu_01, P_01, n_samples)'; x_02 = mvnrnd(mu_02, P_02, n_samples)'; x_final1 = zeros(6, n_samples); x_final2 = zeros(6, n_samples); % Parallelized propagation of orbits parfor i = 1:n_samples options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12); [~, x_out1] = ode45(@TwoBody, tm1, x_01(:, i), options); x_final1(:, i) = x_out1(end, :); [~, x_out2] = ode45(@TwoBody, tm2, x_02(:, i), options); x_final2(:, i) = x_out2(end, :); end % Collision detection using vectorized operations distances = sqrt(sum((x_final1(1:3, :) - x_final2(1:3, :)).^2, 1)); n_collide = sum(distances <= rho_AB); % Probability of collision PC = n_collide / n_samples; % Plot results figure(); scatter3(x_final1(1, :), x_final1(2, :), x_final1(3, :), 'b'); title('Sample 1 Propagated Positions'); xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)'); figure(); scatter3(x_final2(1, :), x_final2(2, :), x_final2(3, :), 'r'); title('Sample 2 Propagated Positions'); xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)'); end function dxdt = TwoBody(~, x) mu = 398600.4415; % km^3/s^2 (Gravitational parameter) r = sqrt(x(1)^2 + x(2)^2 + x(3)^2); dxdt = [x(4); x(5); x(6); -mu * x(1) / r^3; -mu * x(2) / r^3; -mu * x(3) / r^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 trying to compute the probability of collision with the monte carlo method in MATLAB. But I get an error in teh following code. The error is the matrix P_01 must be semi definete and P_02 must be symmetric. How should I fix this problem?
% 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
rho_AB = 15e-3; % Collision threshold (in km)
% Time span setup
a1 = 1/((2/norm(mu_01(1:3))) - ((norm(mu_01(4:6))^2) / mu));
period1 = 2 * pi * sqrt((a1^3) / mu);
tm1 = 0:10:2 * period1;
a2 = 1/((2/norm(mu_02(1:3))) - ((norm(mu_02(4:6))^2) / mu));
period2 = 2 * pi * sqrt((a2^3) / mu);
tm2 = 0:10:2 * period2;
% Monte Carlo Simulation
[PC, n_collide] = MonteCarlo_PC(tm1, P_01, mu_01, tm2, P_02, mu_02, rho_AB);
function [PC, n_collide] = MonteCarlo_PC(tm1, P_01, mu_01, tm2, P_02, mu_02, rho_AB)
n_samples = 1000;
% Generate samples using multivariate normal distribution
x_01 = mvnrnd(mu_01, P_01, n_samples)';
x_02 = mvnrnd(mu_02, P_02, n_samples)';
x_final1 = zeros(6, n_samples);
x_final2 = zeros(6, n_samples);
% Parallelized propagation of orbits
parfor i = 1:n_samples
options = odeset('RelTol', 1e-12, 'AbsTol', 1e-12);
[~, x_out1] = ode45(@TwoBody, tm1, x_01(:, i), options);
x_final1(:, i) = x_out1(end, :);
[~, x_out2] = ode45(@TwoBody, tm2, x_02(:, i), options);
x_final2(:, i) = x_out2(end, :);
end
% Collision detection using vectorized operations
distances = sqrt(sum((x_final1(1:3, :) - x_final2(1:3, :)).^2, 1));
n_collide = sum(distances <= rho_AB);
% Probability of collision
PC = n_collide / n_samples;
% Plot results
figure(); scatter3(x_final1(1, :), x_final1(2, :), x_final1(3, :), 'b');
title('Sample 1 Propagated Positions');
xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)');
figure(); scatter3(x_final2(1, :), x_final2(2, :), x_final2(3, :), 'r');
title('Sample 2 Propagated Positions');
xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)');
end
function dxdt = TwoBody(~, x)
mu = 398600.4415; % km^3/s^2 (Gravitational parameter)
r = sqrt(x(1)^2 + x(2)^2 + x(3)^2);
dxdt = [x(4); x(5); x(6);
-mu * x(1) / r^3;
-mu * x(2) / r^3;
-mu * x(3) / r^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