From a6836b8779c2deb0736daa118efb592d21c966dc Mon Sep 17 00:00:00 2001
From: "Pier Francesco A. Bachini" <pierfrancesco.bachini@skywarder.eu>
Date: Thu, 20 Feb 2025 00:27:55 +0100
Subject: [PATCH] [ADA] Update Montecarlo with run_ada comparison
---
simulator/mainMontecarlo.m | 24 ++++++-
simulator/montecarlo/plotsMontecarlo.m | 71 +++++++++++++-------
simulator/src/std_run_parts/std_subsystems.m | 4 +-
3 files changed, 73 insertions(+), 26 deletions(-)
diff --git a/simulator/mainMontecarlo.m b/simulator/mainMontecarlo.m
index e75023d1..9aeb135d 100644
--- a/simulator/mainMontecarlo.m
+++ b/simulator/mainMontecarlo.m
@@ -115,6 +115,7 @@ for alg_index = 4
motor_K = settings.motor.K;
parfor i = 1:N_sim
+ % for i = 1:N_sim
rocket_vec{i} = copy(rocket);
settings_mont = settings_mont_init;
settings_mont.motor.expThrust = stoch.thrust(i,:); % initialize the thrust vector of the current simulation (parfor purposes)
@@ -314,7 +315,25 @@ for alg_index = 4
p_50 = normcdf([settings.z_final-50, settings.z_final+50],apogee.altitude_mean,apogee.altitude_std);
apogee.accuracy_gaussian_50 =( p_50(2) - p_50(1) )*100;
-
+ %%
+ clc
+ for ii = 1:N_sim
+ if sum(save_thrust{ii}.sensors.ada.apo_times == -1) > 1
+ disp("Problems with ADA " + num2str(ii));
+ end
+ end
+
+ %%
+
+ ada_diff = zeros(N_sim, 1);
+ for ii = 1:N_sim
+ if save_thrust{ii}.sensors.old_ada.t_apogee ~= -1
+ ada_diff(ii) = save_thrust{ii}.sensors.ada.t_apogee - save_thrust{ii}.sensors.old_ada.t_apogee;
+ else
+ ada_diff(ii) = -1;
+ end
+ end
+
%% PLOTS
plotsMontecarlo;
@@ -343,6 +362,9 @@ for alg_index = 4
case {'2024_Lyra_Portugal_October', '2024_Lyra_Roccaraso_September'}
folder = [folder ; "MontecarloResults\"+mission.name+"\"+contSettings.algorithm+"\"+num2str(N_sim)+"sim_Mach"+num2str(100*rocket.airbrakes.maxMach)+"_"+simulationType_thrust+"_"+saveDate]; % offline
+
+ case {'2025_Orion_Portugal_October', '2025_Orion_Roccaraso_September'}
+ folder = [folder ; "MontecarloResults\"+mission.name+"\"+contSettings.algorithm+"\"+num2str(N_sim)+"sim_Mach"+num2str(100*rocket.airbrakes.maxMach)+"_"+simulationType_thrust+"_"+saveDate]; % offline
end
end
diff --git a/simulator/montecarlo/plotsMontecarlo.m b/simulator/montecarlo/plotsMontecarlo.m
index c20efe46..2e16b9dc 100644
--- a/simulator/montecarlo/plotsMontecarlo.m
+++ b/simulator/montecarlo/plotsMontecarlo.m
@@ -3,7 +3,7 @@ N_histCol = min(N_sim,25); % best looking if we don't go higher than 200, but if
sim_per_interval = N_sim/N_histCol;
%% PLOT HISTOGRAM
-montFigures.apogee_histogram_pdf = figure('Color','white','Name','Apogee histogram','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_histogram_pdf = figure('Color','white','Name','Apogee histogram','Units','Normalized','WindowState','maximized');
hold on; grid on;
@@ -16,7 +16,7 @@ ylabel('Number of apogees in the same interval')
title('Reached apogee distribution')
legend
-montFigures.apogee_histogram_cumulative = figure('Color','white','Name','Apogee histogram Cumulative','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_histogram_cumulative = figure('Color','white','Name','Apogee histogram Cumulative','Units','Normalized','WindowState','maximized');
hold on; grid on;
xline(settings.z_final-10, 'r--', 'LineWidth', 1,'DisplayName','-10m from target')
xline(settings.z_final+10, 'r--', 'LineWidth', 1,'DisplayName','+10m from target')
@@ -39,7 +39,7 @@ if ~(strcmp(contSettings.algorithm,'engine')||strcmp(contSettings.algorithm,'NoC
arb_deploy_time_MEAN = mean(arb_deploy_time_vec);
arb_deploy_time_MODE = mode(arb_deploy_time_vec);
% figure
-montFigures.arb_deploy_histogram_pdf = figure('Color','white','Name','Air brakes deployment histogram','Units','Normalized','Position',[0,0,1,1]);
+montFigures.arb_deploy_histogram_pdf = figure('Color','white','Name','Air brakes deployment histogram','Units','Normalized','WindowState','maximized');
histogram(arb_deploy_time_vec,N_histCol)
hold on; grid on;
xline(arb_deploy_time_MEAN,'r--')
@@ -49,7 +49,7 @@ montFigures.arb_deploy_histogram_pdf = figure('Color','white','Name','Air brakes
title("Airbrakes deployment time's distribution")
legend('Airbrakes time deploy','Mean', 'Median')
-montFigures.arb_deploy_histogram_cumulative = figure('Color','white','Name','Air brakes deployment histogram cumulative','Units','Normalized','Position',[0,0,1,1]);
+montFigures.arb_deploy_histogram_cumulative = figure('Color','white','Name','Air brakes deployment histogram cumulative','Units','Normalized','WindowState','maximized');
histogram(arb_deploy_time_vec,N_histCol,'Normalization','cdf')
hold on; grid on;
xline(arb_deploy_time_MEAN,'r--')
@@ -69,17 +69,17 @@ apogee_time_MEAN = mean(apogee_time_vec);
apogee_time_MODE = mode(apogee_time_vec);
% figure
-montFigures.apogee_time_histogram_pdf = figure('Color','white','Name','Apogee time histogram','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_time_histogram_pdf = figure('Color','white','Name','Apogee time histogram','Units','Normalized','WindowState','maximized');
histogram(apogee_time_vec,N_histCol,'DisplayName','Time')
hold on; grid on;
xline(apogee_time_MEAN,'r--','DisplayName','Average')
xline(apogee_time_MODE,'g--','DisplayName','Mode')
-xlabel('Apogee value (m)')
+xlabel('Apogee value (s)')
ylabel('Number of apogees in the same interval')
title('Apogee time distribution')
legend
-montFigures.apogee_time_histogram_cumulative = figure('Color','white','Name','Apogee time histogram cumulative','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_time_histogram_cumulative = figure('Color','white','Name','Apogee time histogram cumulative','Units','Normalized','WindowState','maximized');
histogram(apogee_time_vec,N_histCol,'Normalization','cdf','DisplayName','Time')
hold on; grid on;
xline(apogee_time_MEAN,'r--','DisplayName','Average')
@@ -88,9 +88,32 @@ xlabel('Apogee value (m)')
ylabel('Number of apogees in the same interval')
title('Cumulative apogee time')
legend
+
+%% PLOT COMPARISON BETWEEN RUN_ADA AND MAJORITY VOTING ADA
+
+if contSettings.run_old_ada
+ old_ada_apogee_time_vec = zeros(N_sim, 1);
+ for ii = 1:N_sim
+ old_ada_apogee_time_vec(ii) = save_thrust{ii}.sensors.old_ada.t_apogee;
+ end
+ old_ada_apo_time_mean = mean(old_ada_apogee_time_vec);
+
+ montFigures.apogee_time_ada_comparison = figure('Color','white','Name','Apogee time histogram cumulative','Units','Normalized','WindowState','maximized');
+ histogram(apogee_time_vec, N_histCol, 'DisplayName', 'Majority ADA')
+ hold on; grid on;
+ histogram(old_ada_apogee_time_vec, N_histCol, 'DisplayName', 'run\_ADA');
+ xline(apogee_time_MEAN,'r--','DisplayName','Majority ADA Average')
+ xline(old_ada_apo_time_mean,'r--','DisplayName','run\_ADA Average')
+ xlabel('Apogee value (s)')
+ ylabel('Number of apogees in the same interval')
+ title('Apogee time distribution')
+ legend
+
+end
+
%% PLOT APOGEE MONTECARLO STATISTICS
if settings.scenario ~= "descent"
-montFigures.montecarlo_apogee_statistics = figure('Color','white','Name','Monte Carlo apogee statistics','Units','Normalized','Position',[0,0,1,1]);
+montFigures.montecarlo_apogee_statistics = figure('Color','white','Name','Monte Carlo apogee statistics','Units','Normalized','WindowState','maximized');
subplot(2,1,1)
plot(1:N_sim,apogee_mu,'DisplayName','Increasing mean')
ylabel("Mean value")
@@ -103,7 +126,7 @@ end
%% PLOT LANDING MONTECARLO STATISTICS
if (settings.scenario == "descent" || settings.scenario == "full flight") && settings.parafoil
-montFigures.montecarlo_landing_statistics = figure('Color','white','Name','Monte Carlo landing statistics','Units','Normalized','Position',[0,0,1,1]);
+montFigures.montecarlo_landing_statistics = figure('Color','white','Name','Monte Carlo landing statistics','Units','Normalized','WindowState','maximized');
subplot(2,1,1)
plot(1:N_sim,landing_mu,'DisplayName','Increasing mean')
ylabel("Mean value")
@@ -115,7 +138,7 @@ montFigures.montecarlo_landing_statistics = figure('Color','white','Name','Monte
end
%% PLOT CONTROL
-montFigures.control_action = figure('Color','white','Name','Control action','Units','Normalized','Position',[0,0,1,1]);
+montFigures.control_action = figure('Color','white','Name','Control action','Units','Normalized','WindowState','maximized');
for i = floor(linspace(1,N_sim,5))
plot(save_thrust{i}.t,save_thrust{i}.Y(:,14))
hold on; grid on;
@@ -126,7 +149,7 @@ ylabel('Servo angle \alpha (rad)')
legend(contSettings.algorithm);
%% PLOT APOGEE wrt THRUST
-montFigures.apogee_wrt_thrust = figure('Color','white','Name','Apogee wrt thrust','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_wrt_thrust = figure('Color','white','Name','Apogee wrt thrust','Units','Normalized','WindowState','maximized');
hold on; grid on;
scatter(thrust_percentage,apogee.altitude,'.')
yline(settings.z_final-10,'r--')
@@ -139,7 +162,7 @@ ylim([min(apogee.altitude)-20 max(apogee.altitude)+20])
legend(contSettings.algorithm);
%% PLOT APOGEE wrt MASS OFFSET
-montFigures.apogee_wrt_mass_offset = figure('Color','white','Name','Apogee wrt mass offset','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_wrt_mass_offset = figure('Color','white','Name','Apogee wrt mass offset','Units','Normalized','WindowState','maximized');
hold on; grid on;
scatter(stoch.mass_offset,apogee.altitude,'.')
yline(settings.z_final-10,'r--')
@@ -152,7 +175,7 @@ ylim([min(apogee.altitude)-20 max(apogee.altitude)+20])
legend(contSettings.algorithm);
%% PLOT APOGEE wrt RAIL OMEGA
-montFigures.apogee_wrt_rail_OMEGA= figure('Color','white','Name','Apogee wrt rail elevation','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_wrt_rail_OMEGA= figure('Color','white','Name','Apogee wrt rail elevation','Units','Normalized','WindowState','maximized');
hold on; grid on;
scatter(rad2deg(stoch.OMEGA_rail),apogee.altitude,'.')
yline(settings.z_final-10,'r--')
@@ -167,19 +190,19 @@ legend(contSettings.algorithm);
%% PLOT SHUTDOWN TIME 2D
%%% t_shutdown histogram
if settings.scenario~= "descent"
- montFigures.t_shutdown_histogram_pdf = figure('Color','white','Name','Shutdown time histogram','Units','Normalized','Position',[0,0,1,1]);
+ montFigures.t_shutdown_histogram_pdf = figure('Color','white','Name','Shutdown time histogram','Units','Normalized','WindowState','maximized');
histogram(t_shutdown.value,N_histCol)
xlabel('Shutdown time (s)')
ylabel('Number of shutdowns in the same time interval')
title('Engine shutdown time distribution')
- montFigures.t_shutdown_histogram_cumulative = figure('Color','white','Name','Shutdown time histogram cumulative','Units','Normalized','Position',[0,0,1,1]);
+ montFigures.t_shutdown_histogram_cumulative = figure('Color','white','Name','Shutdown time histogram cumulative','Units','Normalized','WindowState','maximized');
histogram(t_shutdown.value,N_histCol,'Normalization','cdf')
xlabel('Shutdown time (s)')
ylabel('Shutdown time cdf')
title('Engine shutdown time cumulative distribution')
%%% t_shutdown wrt wind
- montFigures.tShutdown_wind = figure('Color','white','Name','Shutdown time wrt wind','Units','Normalized','Position',[0,0,1,1]);
+ montFigures.tShutdown_wind = figure('Color','white','Name','Shutdown time wrt wind','Units','Normalized','WindowState','maximized');
subplot(1,3,1)
for i = 1:N_sim
plot(wind_el(i),save_thrust{i}.sensors.mea.t_shutdown,'.')
@@ -238,7 +261,7 @@ end
%% Predicted apogee
if (strcmp(contSettings.algorithm,'engine') || strcmp(contSettings.algorithm,'complete'))
- montFigures.apogee_prediction = figure('Color','white','Name','Apogee prediction','Units','Normalized','Position',[0,0,1,1]);
+ montFigures.apogee_prediction = figure('Color','white','Name','Apogee prediction','Units','Normalized','WindowState','maximized');
scatter(apogee.prediction_last_time,apogee.prediction,'k','DisplayName','Prediction');
hold on; grid on;
scatter(apogee.prediction_last_time,apogee.altitude,'r','DisplayName','Simulated');
@@ -249,7 +272,7 @@ if (strcmp(contSettings.algorithm,'engine') || strcmp(contSettings.algorithm,'co
end
%% PLOT APOGEE 3D
-montFigures.apogee_3D_wind_thrust = figure('Color','white','Name','Apogee 3D wind thrust','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_3D_wind_thrust = figure('Color','white','Name','Apogee 3D wind thrust','Units','Normalized','WindowState','maximized');
%%%%%%%%%% wind magnitude - thrust - apogee
subplot(2,2,1)
hold on; grid on;
@@ -327,7 +350,7 @@ for i =1:N_sim
max_force_kg(i) = max(abs(force))/9.81;
end
-montFigures.dynamic_pressure_and_forces = figure('Color','white','Name','Dynamic pressure and forces','Units','Normalized','Position',[0,0,1,1]);
+montFigures.dynamic_pressure_and_forces = figure('Color','white','Name','Dynamic pressure and forces','Units','Normalized','WindowState','maximized');
%%%%%%%%%%% max dynamic pressure (Pa)
subplot(2,1,1)
histogram(qdyn_max,N_histCol)
@@ -348,7 +371,7 @@ for i = 1:length(save_thrust)
est_mass(i) = save_thrust{i}.sensors.mea.mass(end);
true_mass(i) = save_thrust{i}.sensors.mea.true_mass_at_shutdown;
end
-montFigures.estimated_mass_histogram_pdf = figure('Color','white','Name','Estimated mass histogram','Units','Normalized','Position',[0,0,1,1]);
+montFigures.estimated_mass_histogram_pdf = figure('Color','white','Name','Estimated mass histogram','Units','Normalized','WindowState','maximized');
histogram(est_mass,N_histCol)
hold on;
histogram(true_mass,N_histCol,'DisplayName','Simulated')
@@ -357,7 +380,7 @@ xlabel('Mass [kg]')
ylabel('Number of simulations')
title('Estimated final mass')
-montFigures.estimated_mass_histogram_cumulative = figure('Color','white','Name','Estimated mass histogram cumulative','Units','Normalized','Position',[0,0,1,1]);
+montFigures.estimated_mass_histogram_cumulative = figure('Color','white','Name','Estimated mass histogram cumulative','Units','Normalized','WindowState','maximized');
histogram(est_mass,N_histCol,'DisplayName','Estimated','Normalization','cdf')
hold on;
histogram(true_mass,N_histCol,'DisplayName','Simulated','Normalization','cdf')
@@ -377,7 +400,7 @@ for ii = 1:N_sim
end
%%%%%%%%%%%%%%%% landing position w.r.t. target
-montFigures.landing_ellipses = figure('Color','white','Name','Landing ellipses','Units','Normalized','Position',[0,0,1,1]);
+montFigures.landing_ellipses = figure('Color','white','Name','Landing ellipses','Units','Normalized','WindowState','maximized');
scatter(landing.position(:,1),landing.position(:,2),'k.','DisplayName','Landings')
hold on;
scatter(settings.payload.target(1),settings.payload.target(2),'DisplayName','Target')
@@ -394,7 +417,7 @@ legend
%% Parafoil trajectories
-montFigures.prf_trajs = figure('Color','white','Name','Parafoil trajectories','Units','Normalized','Position',[0,0,1,1]);
+montFigures.prf_trajs = figure('Color','white','Name','Parafoil trajectories','Units','Normalized','WindowState','maximized');
for ii = 1:N_sim
plot3(save_thrust{ii}.Y(idx_dep:end, 2), save_thrust{ii}.Y(idx_dep:end, 1), -save_thrust{ii}.Y(idx_dep:end, 3));
hold on;
@@ -409,7 +432,7 @@ title('Parafoil trajectories')
% useful for realistic opening loads
% TODO:
% - modify into histogram
-montFigures.apogee_velocity = figure('Color','white','Name','Apogee velocity','Units','Normalized','Position',[0,0,1,1]);
+montFigures.apogee_velocity = figure('Color','white','Name','Apogee velocity','Units','Normalized','WindowState','maximized');
subplot(3,1,1)
plot(apogee.horizontalSpeedX)
title('vn')
diff --git a/simulator/src/std_run_parts/std_subsystems.m b/simulator/src/std_run_parts/std_subsystems.m
index 025f9ceb..1cd97027 100644
--- a/simulator/src/std_run_parts/std_subsystems.m
+++ b/simulator/src/std_run_parts/std_subsystems.m
@@ -8,7 +8,9 @@ This function runs all subsystems in a simulated environment
%% ADA
if settings.flagADA && settings.dataNoise
- % [sensorData, sensorTot, settings.ada, flagApogee, flagOpenPara] = run_ADA(sensorData, sensorTot, settings,t1);
+ if contSettings.run_old_ada
+ [sensorData, sensorTot] = run_ADA(sensorData, sensorTot, settings,t1);
+ end
[sensorData, sensorTot, settings.ada, flagApogee, flagOpenPara] = run_Majority_ADA(Tf, contSettings.ADA_N_instances, settings.ada, settings.frequencies.ADAFrequency, sensorData, sensorTot, currentState, engineT0);
end
--
GitLab