From 487abba0fead1f16033d881198ed52afd8e2b266 Mon Sep 17 00:00:00 2001
From: Giuseppe Brentino <brentinogiuseppe@gmail.com>
Date: Tue, 19 Sep 2023 10:14:52 +0200
Subject: [PATCH] Added merit parameter on air brakes, update rail
 configuration

---
 data/msa-toolkit                              |   2 +-
 .../Preflight_results/Results_preflight.txt   | Bin 1199 -> 1200 bytes
 simulator/configs/configLaunchRail.m          |   4 +-
 simulator/configs/configSimulator.m           |  22 +++++------
 simulator/configs/configWind.m                |  10 ++---
 simulator/mainMontecarlo_FAST.m               |  37 ++++++++++++++++--
 simulator/montecarlo/plotsMontecarlo.m        |  23 +++++++----
 simulator/src_FAST/std_run_FAST.m             |   3 ++
 8 files changed, 70 insertions(+), 31 deletions(-)

diff --git a/data/msa-toolkit b/data/msa-toolkit
index 6f73c9db..0742afba 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
GIT binary patch
delta 155
zcmZ3_xq)-TStdsF$>*5-j15f9^vul_a(TIYGxHP@3-Z%bQ>_#%4Gr{6jlrU+i3njs
z0|QGvqsjiv>Kuj!=6V+9X2z3inT@zi%?$KROcbJXV<*=#YjYQuq$pS^m|N%>8BKoA
utjTI@u4ia6S(ZhE#Z=GCV6rER60^Cv(d0Z9V?INL=v@8cSOs1#5Cj0jP$sAV

delta 154
zcmdnMxt?>wStdra$>*5-j13LU^bE}va(TIYGxHP@3-Z%bQ>_#%Elu?-48WqPi3njs
z15*n<gUSBP>g*;)dPc_PmXqt4jkrt<%=OF-6{2%vC)Y7+a~GGSC|D^Nnd+HYO#Z;E
t$!cV(XKprGjzxpTP|wV4vKNaIv!Sus<a`!mK4XRGT>au$1zs)?1OT9JCV>C|

diff --git a/simulator/configs/configLaunchRail.m b/simulator/configs/configLaunchRail.m
index 0e823cff..703994d0 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 bb87d7c9..49754083 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 c9abdb2f..c217bdcf 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 94ca95a3..06021641 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 8b5f1363..d0c95bd7 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 0cc3f6d7..9d9b9167 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');
 
-- 
GitLab