ROC v2.0: Improved computational efficiency under uncertainty

Author: Abolfazl-Simorgh

.DS_Store

@@ -12,8 +12,6 @@ The Python Library roc is a software package developed by UC3M. It is a tool for
 
 **Support:** Support of all general technical questions on roc, i.e., installation, application, and development, will be provided by Abolfazl Simorgh (abolfazl.simorgh@uc3m.es). 
 
-**Core developer team:** Daniel González Arribas, Abolfazl Simorgh, and Manuel Soler. 
-
 Copyright (C) 2022, Universidad Carlos III de Madrid
 
 ## Citation info
@@ -33,7 +31,7 @@ conda create -n env_roc python=3.9
 conda activate env_roc
 ```
 
-1. Clone or download the repository. The roost source code is available on a public GitHub repository: https://github.com/Abolfazl-Simorgh/roc. The easiest way to obtain it is to clone the repository using git: git clone https://github.com/Abolfazl-Simorgh/roc.git.
+1. Clone or download the repository. The roost source code is available on a public GitHub repository: https://github.com/Aircraft-Operations-Lab/roc. The easiest way to obtain it is to clone the repository using git: git clone https://github.com/Aircraft-Operations-Lab/roc.git.
 
 2. Locate yourself in the roost (library folder) path, and run the following line, using terminal (MacOS and Linux) or cmd (Windows), which will install all dependencies:
 ```python

README.md

@@ -29,10 +29,80 @@ def aCCF_dCont(olr):
 
 
 def aCCF_nCont(t):
-    accf_nCont = 13.9 * 0.0151 * (1e-10 * (0.0073 * 10 ** (0.0107 * t) - 1.03)) * 0.42/ 3
+    accf_nCont = 13.9 * 0.0151 * (1e-10 * (0.0073 * np.power(10, 0.0107 * t) - 1.03)) * 0.42/ 3
     return accf_nCont
 
 
+def G (P):
+    G_ = (1.23 * 1004 * P)/(0.6222 * 43130000 * 0.7)
+    return G_
+
+def T_crit (G):
+    T_c = -46.46 + 9.43 * np.log(G - 0.053) + 0.72 * np.power(np.log(G - 0.053),2) + 273.15
+    return T_c
+
+def rlw (q, P, el_t_amb):
+    r_h = (q * P * (461.51/287.05)) / el_t_amb
+    return r_h
+
+def rh(q, T, p):
+    R_d = 287.05
+    R_v = 461.51
+    rh = (q * p * (R_v / R_d)) / e_sat_liquid(T)
+    return rh
+
+def e_sat_liquid (t):
+    eL_T = 100 * np.exp(-6096.9385/t + 16.635794 - 0.02711193 * t + 1.673952 * 1e-5 * np.power(t, 2) + 2.433502 * np.log(t))
+    return eL_T
+
+
+
+def rhi(q, T, p):
+    R_d = 287.05
+    R_v = 461.51
+    return (q * p * (R_v /R_d)) / e_sat_ice(T)
+
+
+
+def e_sat_ice(t):
+    return 100 * np.exp(  # type: ignore[return-value]
+        (-6024.5282 / t)
+        + 24.7219
+        + (0.010613868 * t)
+        - (1.3198825e-5 * (t**2))
+        - 0.49382577 * np.log(t)
+    )
+
+def r_cont (G, t, t_c, el_t_c, el_t_amb):
+    r_sac = (G * (t - t_c) + el_t_c) / el_t_amb
+    return r_sac
+
+
+T0_MSL = 288.15
+P0_MSL = 101325.0
+lapse_rate_troposphere = -0.0065
+
+def Hp2P(Hp):
+    return P0_MSL*(1 + lapse_rate_troposphere*Hp/T0_MSL)**5.2561
+
+
+def FL2P (FL):
+    Hp = FL * 30.48
+    P = Hp2P(Hp)/100 #hpa
+    return P
+
+
+def sigmd (x, c, x_cut):
+    return 1 / (1 + np.exp( -c * (x-x_cut) ) )
+
+
+def sac (r_sac, r_hm, T, T_crit):
+    T_condition = sigmd (T_crit, 3, T)
+    R_condition = sigmd (r_hm, 20, r_sac)
+    sac_ = R_condition * T_condition
+    return sac_
+    
+
 def pcfa(r, t, r_thr, t_thr, r_pa, t_pa):
     t_condition = 1/(1+np.exp(-t_pa*(-t+t_thr)))
     r_condition = 1/(1+np.exp(-r_pa*( r-r_thr)))

roc3/accf.py

@@ -266,7 +266,13 @@ def patch_mass_and_altitude(fpp: FPP, trj, m0, n_nodes=8000):
         ig.p['min_atr'] = 0
         ig.p['max_atr'] = 0
     if fpp.pcfg['airspeed'] != 'constant' and fpp.pcfg['climate_impact'] and fpp.pcfg['climate_modeling'] == 'mayer' and not fpp.pcfg['climate_std']: 
-        ig.x['atr'] = I(s, atr)  
+        ig.x['atr'] = I(s, atr)
+    ig.p['wp1'] = 593.6373951051227
+    ig.p['wp2'] = 883.8547170355159
+    # ig.p['wp3'] = 1100
+    
+    
+    
     ig.p['lft'] = t[-1]
     ig.p['eft'] = t[-1]
     ig.u['tas_slope'] = I(s, tas_slopes)
@@ -498,6 +504,11 @@ def avg_crs(course_1, course_2):
     ig.x['tas'] = I(distances, VTas)
     ig.x['t_0'] = I(distances, times)
     ig.x['m_0'] = I(distances, masses)
+    ig.p['wp1'] = 593.6373951051227
+    ig.p['wp2'] = 883.8547170355159
+    # ig.p['wp3'] = 1100
+
+
     ig.p['lft'] = times[-1]
     ig.p['eft'] = times[-1]
     ig.u['tas_slope'] = I(distances, accelerations)

roc3/initial_guess.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import casadi
+
+from roc3.bada4 import *
+from roc3.apm import *
+
+from casadi import interpolant, vertcat
+
+class Label(object):
+    def __init__(self, path):
+        self.arrs = np.load(path)
+        self.CL_coords = self.arrs['CL'][0,:]
+        self.M_coords = self.arrs['MACH'][0,:]
+        self.D_values = self.arrs['jetCR.D'][0,:]
+    def get_CD_interpolant(self):
+        xx = [self.CL_coords, self.M_coords]
+        y = self.D_values.T.flatten()
+        I = interpolant('CD', 'bspline', xx, y, {})
+        return lambda CL, M: I(vertcat(CL, M))
+
+class ModifiedBADA(BADA4_jet_CR):
+    def __init__(self, label, array_path):
+        self.isa = ISA()
+        self.label_tabular = Label(array_path)
+        super().__init__(label)
+        self.I_CD = self.label_tabular.get_CD_interpolant()
+   # def fc_from_thrust_P_M_T(self, thrust, P, M, T):
+   #     Hp = P2Hp(P)
+   #     T_ISA = self.isa.IT(Hp)
+   #     uncorrected_fc = self.I_fc(vertcat(thrust, Hp, M))
+   #     return uncorrected_fc*(T/T_ISA)**.5
+    def CD_from_CL_M(self, CL, M):
+        return self.I_CD(CL, M)

roc3/modified_bada.py

@@ -2,9 +2,9 @@
 
 import numpy as np
 
-from roc3.weather import WeatherModel
-from roc3.rfp import *
-from roc3.initial_guess import *
+from .weather import WeatherModel
+from .rfp import *
+from .initial_guess import *
 
 
 class Routing2DStandard(object):
@@ -23,12 +23,14 @@ def __init__(self, apm, weather_model, problem_config, n_nodes, hazards={}, conf
         self.wm_det = weather_model.get_slice_with_first_member()
         self.pcfg = problem_config
         self.n_nodes = n_nodes
+        
         self.config = {
             'n_ig': 200,
-            'n_det1': 20,
-            'n_det2': 40,
-            'n_robust1': 40,
-            }
+            'n_det1': 30,
+            'n_det2': 60,
+            'n_robust1': 60,
+            }        
+        
         self.config.update(config)
         # solver_options ['mu_init'] = 10**-9
         self.solver_options = solver_options
@@ -133,6 +135,7 @@ def solve(self, stop=False, estimations={}):
         self.check_trj_in_wm_bounds(new_ig, self.wm_det)
         fpp_det.build_ocp()
         trj = fpp_det.solve(new_ig, n_nodes=self.config['n_det2'], solver_options=self.solver_options)
+        # trj = fpp_det.solve(new_ig, n_nodes=self.config['n_det2'], solver_options={'max_iter':0})
         # fpp_det.autoplot()
         self.trjs_steps.append(trj)
 
@@ -142,14 +145,11 @@ def solve(self, stop=False, estimations={}):
 
         new_ig = trj.get_interpolator(patch_tu=True)
         so = self.solver_options.copy()
-        so['mu_init'] = 10**-3.8
-        so['max_iter'] = 1500
-        trj = fpp.solve(new_ig, n_nodes=self.config['n_robust1'], solver_options=so)
+        so['max_iter'] = 2000
+        trj = fpp.solve(new_ig, n_nodes=self.config['n_det2'], solver_options=so)
         self.trjs_steps.append(trj)
 
-
         new_ig = trj.get_interpolator(patch_tu=True)
-        so['mu_init'] = 10**-3.8
         trj = fpp.solve(new_ig, n_nodes=self.n_nodes, solver_options=so)
         self.trjs_steps.append(trj)