diff --git a/data/msa-toolkit b/data/msa-toolkit index 6f73c9db5ea4a2e38e6cff1c65931289b99382ac..0742afba9df49ed58374fe707eb787dbd0db47ad 160000 --- a/data/msa-toolkit +++ b/data/msa-toolkit @@ -1 +1 @@ -Subproject commit 6f73c9db5ea4a2e38e6cff1c65931289b99382ac +Subproject commit 0742afba9df49ed58374fe707eb787dbd0db47ad diff --git a/simulator/Preflight_results/Results_preflight.txt b/simulator/Preflight_results/Results_preflight.txt index 346ac378aadbfe9eb761253ee397bf564ece1fa8..27222b4e0d73f2987ad5bcad76cd1092360f704d 100644 Binary files a/simulator/Preflight_results/Results_preflight.txt and b/simulator/Preflight_results/Results_preflight.txt differ diff --git a/simulator/configs/configLaunchRail.m b/simulator/configs/configLaunchRail.m index 0e823cffe197f7e8417d19cd76cd4cad27e592ff..703994d03308764b478d455d6c338e426fbfbea6 100644 --- a/simulator/configs/configLaunchRail.m +++ b/simulator/configs/configLaunchRail.m @@ -7,5 +7,5 @@ launch rail config script % launchpad directions % for a single run the maximum and the minimum value of the following angles must be the same. -settings.OMEGA = 80*pi/180; % [rad] Minimum Elevation Angle, user input in degrees (ex. 80) -settings.PHI = 168*pi/180; % [rad] Minimum Azimuth Angle from North Direction, user input in degrees (ex. 90) +settings.OMEGA = 77*pi/180; % [rad] Minimum Elevation Angle, user input in degrees (ex. 80) +settings.PHI = 177*pi/180; % [rad] Minimum Azimuth Angle from North Direction, user input in degrees (ex. 90) diff --git a/simulator/configs/configSimulator.m b/simulator/configs/configSimulator.m index bb87d7c91a5a3d1edbfa9fa7aed7fabd7c751439..49754083f4b8a16b95c6393e712964a1c14d25f1 100644 --- a/simulator/configs/configSimulator.m +++ b/simulator/configs/configSimulator.m @@ -41,17 +41,17 @@ configWind; configPlots; %% set real time parameters -% I_TOT = 4718.86; -% BURNING_TIME = 2.90897; -% %%% ----------------------------- -% BURNING_TIME = BURNING_TIME + settings.tIGN + settings.tCO; -% timeNew = settings.motor.expTime * BURNING_TIME/settings.tb; -% Itemp = trapz(timeNew, settings.motor.expThrust); -% ThrustNew = settings.motor.expThrust * I_TOT/Itemp; -% settings.State.xcgTime = settings.State.xcgTime * timeNew(end)/settings.tb; -% settings.tb = timeNew(end); -% settings.motor.expTime = timeNew; -% settings.motor.expThrust = ThrustNew; +I_TOT = 4861.46; +BURNING_TIME = 3.13561; +%%% ----------------------------- +BURNING_TIME = BURNING_TIME + settings.tIGN + settings.tCO; +timeNew = settings.motor.expTime * BURNING_TIME/settings.tb; +Itemp = trapz(timeNew, settings.motor.expThrust); +ThrustNew = settings.motor.expThrust * I_TOT/Itemp; +settings.State.xcgTime = settings.State.xcgTime * timeNew(end)/settings.tb; +settings.tb = timeNew(end); +settings.motor.expTime = timeNew; +settings.motor.expThrust = ThrustNew; %% montecarlo settings configMontecarlo; diff --git a/simulator/configs/configWind.m b/simulator/configs/configWind.m index c9abdb2fbc1e7ac42418917c5c13c19eec2c0647..c217bdcf82d53b61548b35cc72775d81e2684a59 100644 --- a/simulator/configs/configWind.m +++ b/simulator/configs/configWind.m @@ -41,14 +41,14 @@ if conf.script == "simulator" % Wind is generated for every altitude interpolating with the coefficient defined below settings.wind.input = true; settings.wind.model = false; - settings.wind.inputGround = 6; % [m/s] Wind magnitude at the ground - settings.wind.inputAlt = 2*linspace(0, 2000, 100); % [m] Altitude vector + settings.wind.inputGround = 7; % [m/s] Wind magnitude at the ground + settings.wind.inputAlt = linspace(0, 2000, 100); % [m] Altitude vector % settings.wind.inputMult = [0 1 2 3 4 4.5 5 5.5 6]*50; % [-] Percentage of increasing magnitude at each altitude % settings.wind.inputAzimut = 200*pi/180*ones(1,9); % [deg] Wind azimut angle at each altitude (toward wind incoming direction) % settings.wind.input_uncertainty = [0,0]; - settings.wind.inputMagnitude = [2 * ones(1, 15), 3 * ones(1, 5), 4 * ones(1, 10), 5*ones(1, 15), 6*ones(1, 20), 7 * ones(1, 35)]; - settings.wind.inputAzimut = [150 * ones(1, 15), 310 * ones(1, 85)]; % wind azimut angle at each altitude (toward wind incoming direction) [deg] - settings.wind.input_uncertainty = [50, 30]; + settings.wind.inputMagnitude = [7 * ones(1, 15), 8 * ones(1, 5), 8 * ones(1, 10), 11*ones(1, 15), 13*ones(1, 20), 15 * ones(1, 35)]; + settings.wind.inputAzimut = [(177+180) * ones(1, 15), 270 * ones(1, 85)]; % wind azimut angle at each altitude (toward wind incoming direction) [deg] + settings.wind.input_uncertainty = [10, 20]; settings.wind.inputMult = (settings.wind.inputMagnitude./settings.wind.inputGround - 1) * 100; % settings.wind.input_uncertainty = [a,b]; wind uncertanties: % - a, wind magnitude percentage uncertanty: magn = magn *(1 +- a) diff --git a/simulator/mainMontecarlo_FAST.m b/simulator/mainMontecarlo_FAST.m index 94ca95a3c579179e2bc823b7fea61388b4b0365a..060216411bd29e6988f80c7a0d5b6c667a24e48b 100644 --- a/simulator/mainMontecarlo_FAST.m +++ b/simulator/mainMontecarlo_FAST.m @@ -188,8 +188,31 @@ apogee.horizontalSpeed_min = min(apogee.horizontalSpeed); apogee.accuracy_10 = N_ApogeeWithinTarget_10/N_sim*100; % percentage, so*100 apogee.accuracy_50 = N_ApogeeWithinTarget_50/N_sim*100; % percentage, so*100 +%% compute how many simulations open air brakes +on = 0; +on_max = 0; +off = 0; +for i = 1:N_sim + if max(save_thrust{i}.ARB.cmdPosition) > 0 + on = on+1; + if max(save_thrust{i}.ARB.cmdPosition) >= settings.servo.maxAngle-1e-3 + on_max = on_max + 1; + end + else + off = off+1; + end +end - +meritParam = zeros(N_sim,1); +for i = 1:N_sim + totalABKTime = save_thrust{i}.ARB.cmdTime(find(save_thrust{i}.ARB.cmdPosition>0,1,'last'))-save_thrust{i}.ARB.cmdTime(find(save_thrust{i}.ARB.cmdPosition>0,1,'first')); + totalABKPosition = mean(save_thrust{i}.ARB.cmdPosition)/settings.servo.maxAngle; + if ~isnan(totalABKPosition) && ~isempty(totalABKTime) + meritParam(i) = totalABKPosition*totalABKTime/(max(save_thrust{i}.t)-contSettings.ABK_shadowmode); + else + meritParam(i) = 0; + end +end %% PLOTS @@ -197,9 +220,11 @@ plotsMontecarlo; %% print fprintf('\n') -fprintf('MIN shutdown time = %.3f\n',max(t_shutdown.value)) -fprintf('MAX shutdown time = %.3f\n',min(t_shutdown.value)) -fprintf('MEAN shutdown time = %.3f\n',t_shutdown.mean) +fprintf('MIN shutdown time = %.3f\n',min(t_shutdown.value)) +fprintf('MAX shutdown time = %.3f\n',max(t_shutdown.value)) +fprintf('MEAN shutdown time = %.3f\n\n',t_shutdown.mean) +fprintf('%% simulations abk open = %.2f%%\n',100*on/N_sim) +fprintf('%% simulations merit param above 0.1= %.2f%%\n\n',100*sum(meritParam>0.1)/N_sim) fprintf('Computation time %.3f\n', toc) @@ -282,3 +307,7 @@ end fclose(fid); +%% save file +% saveas(gca,azwind+config"om"+num2str(settings.OMEGA)) + + diff --git a/simulator/montecarlo/plotsMontecarlo.m b/simulator/montecarlo/plotsMontecarlo.m index 8b5f1363a764d8beb75a810f3bbd3370f84e17e4..d0c95bd7163b87ec0df3b6387ed9dbe50f2a6e2e 100644 --- a/simulator/montecarlo/plotsMontecarlo.m +++ b/simulator/montecarlo/plotsMontecarlo.m @@ -19,14 +19,21 @@ legend % PLOT CONTROL subplot(2,2,2) -for i = floor(linspace(1,N_sim,10)) - plot(save_thrust{i}.t,save_thrust{i}.Y(:,14)) - hold on; grid on; -end -title('Control action') -xlabel('Time [s]') -ylabel('Servo angle \alpha [rad]') -legend(contSettings.algorithm); +scatter(1:N_sim,meritParam) +hold on; +yline(mean(meritParam),'r--',['Mean = ',num2str(mean(meritParam))]) + + + +% for i = floor(linspace(1,N_sim,10)) +% plot(save_thrust{i}.t,save_thrust{i}.Y(:,14)) +% hold on; grid on; +% end +% title('Control action') +% xlabel('Time [s]') +% ylabel('Servo angle \alpha [rad]') +% legend(contSettings.algorithm); + %% t_shutdown % I am adding 0.3 seconds because we do not consider it in simulation, while in reality there is a delay between the command and the actuation subplot(2,2,3) diff --git a/simulator/src_FAST/std_run_FAST.m b/simulator/src_FAST/std_run_FAST.m index 0cc3f6d7665760705ddee2addf1eaac252891b5a..9d9b91671f702cb74ed59c6841693092974bb62f 100644 --- a/simulator/src_FAST/std_run_FAST.m +++ b/simulator/src_FAST/std_run_FAST.m @@ -462,6 +462,9 @@ if settings.montecarlo % air brakes (ARB) struct_out.ARB.cmdPosition = interp1(struct_out.ARB.cmdTime,struct_out.ARB.cmdPosition,t_vec); + if any(isnan(struct_out.ARB.cmdPosition)) + struct_out.ARB.cmdPosition(isnan(struct_out.ARB.cmdPosition)) = 0; + end struct_out.ARB.cmdTime = t_vec; struct_out.ARB = rmfield(struct_out.ARB, 'allowanceIdx');