Coverage for src/amisc/examples/models.py: 94%
175 statements
« prev ^ index » next coverage.py v7.6.1, created at 2025-01-24 04:51 +0000
1"""Provides some example test functions for demonstration purposes."""
2import pickle
3import uuid
4from pathlib import Path
6import numpy as np
8from amisc.system import System
11def f1(x):
12 """Figure 5 in Jakeman 2022."""
13 y1 = x * np.sin(np.pi * x)
14 return y1
17def f2(y1):
18 """Figure 5 in Jakeman 2022."""
19 y2 = 1 / (1 + 25*y1**2)
20 return y2
23def f3(x, y2):
24 """Hypothetical third function for Figure 5 in Jakeman 2022."""
25 y3 = x * np.cos(np.pi * y2)
26 return y3
29def tanh_func(inputs, A=2, L=1, frac=4):
30 """Simple tunable tanh function."""
31 return {'y': A * np.tanh(2 / (L/frac) * (inputs['x'] - L/2)) + A + 1}
34def borehole_func(inputs):
35 """Model found at https://www.sfu.ca/~ssurjano/borehole.html
37 :param inputs: `dict` of input variables `rw`, `r`, `Tu`, `Hu`, `Tl`, `Hl`, `L`, and `Kw`
38 :returns vdot: Water flow rate in m^3/yr
39 """
40 rw = inputs['rw'] # Radius of borehole (m)
41 r = inputs['r'] # Radius of influence (m)
42 Tu = inputs['Tu'] # Transmissivity (m^2/yr)
43 Hu = inputs['Hu'] # Potentiometric head (m)
44 Tl = inputs['Tl'] # Transmissivity (m^2/yr)
45 Hl = inputs['Hl'] # Potentiometric head (m)
46 L = inputs['L'] # Length of borehole (m)
47 Kw = inputs['Kw'] # Hydraulic conductivity (m/yr)
48 vdot = 2*np.pi*Tu*(Hu-Hl) / (np.log(r/rw) * (1 + (2*L*Tu/(np.log(r/rw)*Kw*rw**2)) + (Tu/Tl)))
50 return {'vdot': vdot}
53def wing_weight_func(inputs):
54 """Model found at https://www.sfu.ca/~ssurjano/wingweight.html
56 :param inputs: `dict` of input variables `Sw`, `Wfw`, `A`, `Lambda`, `q`, `lamb`, `tc`, `Nz`, `Wdg`, and `Wp`
57 :returns Wwing: the weight of the airplane wing (lb)
58 """
59 Sw = inputs['Sw'] # Wing area (ft^2)
60 Wfw = inputs['Wfw'] # Weight of fuel (lb)
61 A = inputs['A'] # Aspect ratio
62 Lambda = inputs['Lambda'] # Quarter-chord sweep (deg)
63 q = inputs['q'] # Dynamic pressure (lb/ft^2)
64 lamb = inputs['lambda'] # taper ratio
65 tc = inputs['tc'] # Aerofoil thickness to chord ratio
66 Nz = inputs['Nz'] # Ultimate load factor
67 Wdg = inputs['Wdg'] # Design gross weight (lb)
68 Wp = inputs['Wp'] # Paint weight (lb/ft^2)
69 Lambda = Lambda*(np.pi/180)
70 Wwing = 0.036*(Sw**0.758)*(Wfw**0.0035)*((A/(np.cos(Lambda))**2)**0.6)*(q**0.006)*(lamb**0.04)*\
71 (100*tc/np.cos(Lambda))**(-0.3)*((Nz*Wdg)**0.49) + Sw*Wp
73 return {'Wwing': Wwing}
76def nonlinear_wave(inputs, env_var=0.1**2, wavelength=0.5, wave_amp=0.1, tanh_amp=0.5, L=1, t=0.25):
77 """Custom nonlinear model of a traveling Gaussian wave for testing.
79 :param inputs: `dict` of input variables `d` and `theta`
80 :param env_var: variance of Gaussian envelope
81 :param wavelength: sinusoidal perturbation wavelength
82 :param wave_amp: amplitude of perturbation
83 :param tanh_amp: amplitude of tanh(x)
84 :param L: domain length of underlying tanh function
85 :param t: transition length of tanh function (as fraction of L)
86 :returns: `dict` with model output `y`
87 """
88 d = inputs['d']
89 theta = inputs['theta']
91 # Traveling sinusoid with moving Gaussian envelope
92 env_range = [0.2, 0.6]
93 mu = env_range[0] + theta * (env_range[1] - env_range[0])
94 theta_env = 1 / (np.sqrt(2 * np.pi * env_var)) * np.exp(-0.5 * (d - mu) ** 2 / env_var)
95 ftheta = wave_amp * np.sin((2*np.pi/wavelength) * theta) * theta_env
97 # Underlying tanh dependence on d
98 fd = tanh_amp * np.tanh(2/(L*t)*(d - L/2)) + tanh_amp
100 # Compute model = f(theta, d) + f(d)
101 return {'y': ftheta + fd}
104def fire_sat_globals():
105 """Global variables for the fire satellite system."""
106 Re = 6378140. # Radius of Earth (m)
107 mu = 3.986e14 # Gravitational parameter (m^3 s^-2)
108 eta = 0.22 # Power efficiency
109 Id = 0.77 # Inherent degradation of the array
110 thetai = 0. # Sun incidence angle
111 LT = 15. # Spacecraft lifetime (years)
112 eps = 0.0375 # Power production degradation (%/year)
113 rlw = 3. # Length to width ratio
114 nsa = 3. # Number of solar arrays
115 rho = 700. # Mass density of arrays (kg/m^3)
116 t = 0.005 # Thickness (m)
117 D = 2. # Distance between panels (m)
118 IbodyX = 6200. # kg*m^2
119 IbodyY = 6200. # kg*m^2
120 IbodyZ = 4700. # kg*m^2
121 dt_slew = 760. # s
122 theta = 15. # Deviation of moment axis from vertical (deg)
123 As = 13.85 # Area reflecting radiation (m^2)
124 c = 2.9979e8 # Speed of light (m/s)
125 M = 7.96e15 # Magnetic moment of earth (A*m^2)
126 Rd = 5. # Residual dipole of spacecraft (A*m^2)
127 rhoa=5.148e-11 # Atmospheric density (kg/m^3) -- typo in Chaudhuri 2018 has this as 1e11 instead
128 A = 13.85 # Cross-section in flight (m^2)
129 Phold = 20. # Holding power (W)
130 omega = 6000. # Max vel of wheel (rpm)
131 nrw = 3. # Number of reaction wheels
133 return {k: v for k, v in locals().items() if not k.startswith('_') and isinstance(v, float)}
136def orbit_fun(inputs, output_path=None, pct_failure=0):
137 """Compute the orbit model for the fire satellite system.
139 :param inputs: `dict` of input variables `H` and `Φ`
140 :param output_path: where to save model outputs
141 :param pct_failure: probability of a failure
142 :returns: `dict` with model outputs `Vsat`, `To`, `Te`, and `Slew`
143 """
144 v = fire_sat_globals() # Global variables
146 H = inputs['H'] # Altitude (m)
147 phi = inputs[u'Φ'] # Target diameter (m)
148 vel = np.sqrt(v['mu'] / (v['Re'] + H)) # Satellite velocity (m/s)
149 dt_orbit = 2*np.pi*(v['Re'] + H) / vel # Orbit period (s)
150 dt_eclipse = (dt_orbit/np.pi)*np.arcsin(v['Re'] / (v['Re'] + H)) # Eclipse period (s)
151 theta_slew = np.arctan(np.sin(phi / v['Re']) / (1 - np.cos(phi / v['Re']) + H/v['Re'])) # Max slew angle (rad)
153 num_samples = np.atleast_1d(H).shape[0]
155 if np.random.rand() < pct_failure:
156 i = np.random.randint(0, num_samples)
157 i2 = np.random.randint(0, num_samples)
158 vel[i] = np.nan
159 theta_slew[i2] = np.nan
161 y = {'Vsat': vel, 'To': dt_orbit, 'Te': dt_eclipse, 'Slew': theta_slew}
163 if output_path is not None:
164 files = []
165 id = str(uuid.uuid4())
166 for index in range(num_samples):
167 fname = f'{id}_{index}.pkl'
168 with open(Path(output_path) / fname, 'wb') as fd:
169 pickle.dump({var: np.atleast_1d(y[var])[index] for var in y.keys()}, fd)
170 files.append(fname)
171 y['output_path'] = files
173 return y
176def power_fun(inputs, model_fidelity=(0,), *, output_path=None, pct_failure=0):
177 """Compute the power model for the fire satellite system.
179 :param inputs: `dict` of input variables `Po`, `Fs`, `To`, `Te`, and `Pat`
180 :param model_fidelity: fidelity index for the model
181 :param output_path: where to save model outputs
182 :param pct_failure: probability of a failure
183 :returns: `dict` with model outputs `Imin`, `Imax`, `Ptot`, and `Asa`
184 """
185 alpha = model_fidelity
186 pct = 1 - (2 - alpha[0]) * 0.04 if len(alpha) == 1 else 1 # extra pct error term
187 v = fire_sat_globals() # Global variables
189 Po = inputs['Po'] # Other power sources (W)
190 Fs = inputs['Fs'] # Solar flux (W/m^2)
191 dt_orbit = inputs['To'] # Orbit period (s)
192 dt_eclipse = inputs['Te'] # Eclipse period (s)
193 Pacs = inputs['Pat'] # Power from attitude control system (W)
195 Ptot = Po + Pacs
196 Pe = Ptot
197 Pd = Ptot
198 Xe = 0.6 # These are power efficiencies in eclipse and daylight
199 Xd = 0.8 # See Ch. 11 of Wertz 1999 SMAD
200 Te = dt_eclipse
201 Td = dt_orbit - Te
202 Psa = ((Pe*Te/Xe) + (Pd*Td/Xd)) / Td
203 Pbol = v['eta'] * Fs * v['Id'] * np.cos(v['thetai'])
204 Peol = Pbol * (1 - v['eps']) ** v['LT']
205 Asa = Psa / Peol # Total solar array size (m^2)
206 L = np.sqrt(Asa * v['rlw'] / v['nsa'])
207 W = np.sqrt(Asa / (v['rlw'] * v['nsa']))
208 msa = 2 * v['rho'] * L * W * v['t'] # Mass of solar array
209 Ix = msa*((1/12)*(L**2 + v['t']**2) + (v['D']+L/2)**2) + v['IbodyX']
210 Iy = (msa/12)*(L**2 + W**2) + v['IbodyY'] # typo in Zaman 2013 has this as H**2 instead of L**2
211 Iz = msa*((1/12)*(L**2 + W**2) + (v['D'] + L/2)**2) + v['IbodyZ']
212 Itot = np.concatenate((Ix[..., np.newaxis], Iy[..., np.newaxis], Iz[..., np.newaxis]), axis=-1)
213 Imin = np.min(Itot, axis=-1)
214 Imax = np.max(Itot, axis=-1)
216 num_samples = np.atleast_1d(Po).shape[0]
218 if np.random.rand() < pct_failure:
219 i = np.random.randint(0, num_samples)
220 i2 = np.random.randint(0, num_samples)
221 Imin[i2] = np.nan
222 Asa[i] = np.nan
224 y = {'Imin': Imin, 'Imax': Imax*pct, 'Ptot': Ptot*pct, 'Asa': Asa}
226 if output_path is not None:
227 files = []
228 id = str(uuid.uuid4())
229 for index in range(num_samples):
230 fname = f'{id}_{index}.pkl'
231 with open(Path(output_path) / fname, 'wb') as fd:
232 pickle.dump({var: np.atleast_1d(y[var])[index] for var in y.keys()}, fd)
233 files.append(fname)
234 y['output_path'] = files
236 return y
239def attitude_fun(inputs, model_fidelity=(0,)):
240 """Compute the attitude model for the fire satellite system.
242 :param inputs: `dict` of input variables `H`, `Fs`, `Lsp`, `q`, `La`, `Cd`, `Vsat`, and `Slew`
243 :param model_fidelity: fidelity index for the model
244 :returns: `dict` with model outputs `Pat` and `tau_tot`
245 """
246 alpha = model_fidelity
247 v = fire_sat_globals() # Global variables
248 pct = 1 - (2 - alpha[0])*0.04 if len(alpha) == 1 else 1 # extra model fidelity pct error term
250 H = inputs['H'] # Altitude (m)
251 Fs = inputs['Fs'] # Solar flux
252 Lsp = inputs['Lsp'] # Moment arm for solar radiation pressure
253 q = inputs['q'] # Reflectance factor
254 La = inputs['La'] # Moment arm for aerodynamic drag
255 Cd = inputs['Cd'] # Drag coefficient
256 vel = inputs['Vsat'] # Satellite velocity
257 theta_slew = inputs['Slew'] # Max slew angle
258 Imin = inputs['Imin'] # Minimum moment of inertia
259 Imax = inputs['Imax'] # Maximum moment of inertia
261 tau_slew = 4*theta_slew*Imax / v['dt_slew']**2
262 tau_g = 3 * v['mu'] * np.abs(Imax - Imin) * np.sin(2 * v['theta'] * (np.pi / 180)) / (2 * (v['Re'] + H) ** 3)
263 tau_sp = Lsp * Fs * v['As'] * (1 + q) * np.cos(v['thetai']) / v['c']
264 tau_m = 2 * v['M'] * v['Rd'] / (v['Re'] + H) ** 3
265 tau_a = (1 / 2) * La * v['rhoa'] * Cd * v['A'] * vel ** 2
266 tau_dist = np.sqrt(tau_g**2 + tau_sp**2 + tau_m**2 + tau_a**2)
267 tau_tot = np.max(np.concatenate((tau_slew[..., np.newaxis], tau_dist[..., np.newaxis]), axis=-1), axis=-1)
268 Pacs = tau_tot * (v['omega'] * (2 * np.pi / 60)) + v['nrw'] * v['Phold']
270 y = {'Pat': Pacs*pct, 'tau_tot': tau_tot*pct}
272 return y
275def fire_sat_system(save_dir=None):
276 """Fire satellite detection system model from Chaudhuri 2018.
278 !!! Note "Some modifications"
279 Orbit will save outputs all the time; Power won't because it is part of an FPI loop. Orbit and Power can
280 return `np.nan` some of the time (optional, to test robustness of surrogate). Power and Attitude have an
281 `alpha` fidelity index that controls accuracy of the returned values. `alpha=0,1,2` corresponds to `8,4,0`
282 percent error.
284 :param save_dir: where to save model outputs
285 :returns: a `System` object for the fire sat MD example system
286 """
287 return System.load_from_file(Path(__file__).parent / 'fire-sat.yml', root_dir=save_dir)