Skip to content

Commit

Permalink
Merge pull request OpenMDAO#578 from xjjiang/smooth_minimum
Browse files Browse the repository at this point in the history 
Smoothen minimum function in propeller performance
  • Loading branch information
@crecine
crecine authored Oct 23, 2024
2 parents f56d89a + 2f9519a commit 990527d
Show file tree
Hide file tree
Showing 3 changed files with 289 additions and 37 deletions.
10 changes: 5 additions & 5 deletions aviary/subsystems/propulsion/propeller/hamilton_standard.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -635,7 +635,7 @@ def compute(self, inputs, outputs):
TXCLI = np.zeros(6)
CTTT = np.zeros(4)
XXXFT = np.zeros(4)
act_factor = inputs[Aircraft.Engine.PROPELLER_ACTIVITY_FACTOR]
act_factor = inputs[Aircraft.Engine.PROPELLER_ACTIVITY_FACTOR][0]
for k in range(2):
AF_adj_CP[k], run_flag = _unint(Act_Factor_arr, AFCPC[k], act_factor)
AF_adj_CT[k], run_flag = _unint(Act_Factor_arr, AFCTC[k], act_factor)
Expand Down Expand Up @@ -667,7 +667,7 @@ def compute(self, inputs, outputs):
# flag that given lift coeff (cli) does not fall on a node point of CL_arr
CL_tab_idx_flg = 0 # NCL_flg
ifnd = 0
cli = inputs[Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICIENT]
cli = inputs[Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICIENT][0]
power_coefficient = inputs['power_coefficient'][i_node]
for ii in range(6):
cl_idx = ii
Expand Down Expand Up @@ -740,7 +740,7 @@ def compute(self, inputs, outputs):
CL_tab_idx = CL_tab_idx+1
if (CL_tab_idx_flg != 1):
PCLI, run_flag = _unint(
CL_arr[CL_tab_idx_begin:CL_tab_idx_begin+4], PXCLI[CL_tab_idx_begin:CL_tab_idx_begin+4], inputs[Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICIENT])
CL_arr[CL_tab_idx_begin:CL_tab_idx_begin+4], PXCLI[CL_tab_idx_begin:CL_tab_idx_begin+4], inputs[Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICIENT][0])
else:
PCLI = PXCLI[CL_tab_idx_begin]
# PCLI = CLI adjustment to power_coefficient
Expand All @@ -752,7 +752,7 @@ def compute(self, inputs, outputs):
try:
CTT[kdx], run_flag = _unint(
# thrust coeff at baseline point for kdx
Blade_angle_table[kdx], CT_Angle_table[idx_blade][kdx][:ang_len], BLL[kdx])
Blade_angle_table[kdx][:ang_len], CT_Angle_table[idx_blade][kdx][:ang_len], BLL[kdx])
except IndexError:
raise om.AnalysisError(
"interp failed for CTT (thrust coefficient) in hamilton_standard.py")
Expand Down Expand Up @@ -810,7 +810,7 @@ def compute(self, inputs, outputs):
XFFT[kl], run_flag = _biquad(comp_mach_CT_arr, 1, DMN, CTE2)
CL_tab_idx = CL_tab_idx + 1
if (CL_tab_idx_flg != 1):
cli = inputs[Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICIENT]
cli = inputs[Aircraft.Engine.PROPELLER_INTEGRATED_LIFT_COEFFICIENT][0]
TCLII, run_flag = _unint(
CL_arr[CL_tab_idx_begin:CL_tab_idx_begin+4], TXCLI[CL_tab_idx_begin:CL_tab_idx_begin+4], cli)
xft, run_flag = _unint(
Expand Down
200 changes: 170 additions & 30 deletions aviary/subsystems/propulsion/propeller/propeller_performance.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,40 @@
from aviary.variable_info.variables import Aircraft, Dynamic


def smooth_min(x, b, alpha=100.0):
"""
Smooth approximation of the min function using the log-sum-exp trick.
Parameters:
x (float or array-like): First value.
b (float or array-like): Second value.
alpha (float): The smoothing factor. Higher values make it closer to the true minimum. Try between 75 and 275.
Returns:
float or array-like: The smooth approximation of min(x, b).
"""
sum_log_exp = np.log(np.exp(np.multiply(-alpha, x)) + np.exp(np.multiply(-alpha, b)))
rv = -(1 / alpha) * sum_log_exp
return rv


def d_smooth_min(x, b, alpha=100.0):
"""
Derivative of function smooth_min(x)
Parameters:
x (float or array-like): First value.
b (float or array-like): Second value.
alpha (float): The smoothing factor. Higher values make it closer to the true minimum. Try between 75 and 275.
Returns:
float or array-like: The smooth approximation of derivative of min(x, b).
"""
d_sum_log_exp = np.exp(np.multiply(-alpha, x)) / \
(np.exp(np.multiply(-alpha, x)) + np.exp(np.multiply(-alpha, b)))
return d_sum_log_exp


class TipSpeedLimit(om.ExplicitComponent):
"""
Computation of propeller tip speed.
Expand Down Expand Up @@ -305,6 +339,136 @@ def compute_partials(self, inputs, J):
np.sqrt(helical_mach * helical_mach - mach * mach)


class AreaSquareRatio(om.ExplicitComponent):
"""
Compute the area ratio nacelle and propeller with a maximum 0.5.
"""

def initialize(self):
self.options.declare("num_nodes", types=int)
self.options.declare('smooth_sqa', default=True, types=bool)
self.options.declare('alpha', default=100.0, types=float)

def setup(self):
nn = self.options["num_nodes"]
arange = np.arange(self.options["num_nodes"])
self.add_input("DiamNac", val=0.0, units='ft')
self.add_input("DiamProp", val=0.0, units='ft')

self.add_output('sqa_array', val=np.zeros(nn), units='unitless')

self.declare_partials("sqa_array",
[
"DiamNac",
"DiamProp",
],
rows=arange, cols=np.zeros(nn))

def compute(self, inputs, outputs):
nn = self.options["num_nodes"]
diamNac = inputs["DiamNac"]
diamProp = inputs["DiamProp"]
sqa = diamNac**2 / diamProp**2

smooth = self.options["smooth_sqa"]
if smooth:
alpha = self.options['alpha']
sqa = smooth_min(sqa, 0.50, alpha)
else:
sqa = np.minimum(sqa, 0.50)
outputs["sqa_array"] = np.ones(nn) * sqa

def compute_partials(self, inputs, partials):
diamNac = inputs["DiamNac"]
diamProp = inputs["DiamProp"]
sqa = diamNac**2 / diamProp**2

dSQA_dNacDiam = 2 * diamNac / diamProp**2
dSQA_dPropDiam = -2 * diamNac**2 / diamProp**3

smooth = self.options["smooth_sqa"]
if smooth:
alpha = self.options['alpha']
dSQA_dNacDiam = d_smooth_min(sqa, 0.50, alpha) * dSQA_dNacDiam
dSQA_dPropDiam = d_smooth_min(sqa, 0.50, alpha) * dSQA_dPropDiam
else:
dSQA_dNacDiam = np.piecewise(
sqa, [sqa < 0.5, sqa >= 0.5], [1, 0]) * dSQA_dNacDiam
dSQA_dPropDiam = np.piecewise(
sqa, [sqa < 0.5, sqa >= 0.5], [1, 0]) * dSQA_dPropDiam
partials['sqa_array', "DiamNac"] = dSQA_dNacDiam
partials['sqa_array', "DiamProp"] = dSQA_dPropDiam


class AdvanceRatio(om.ExplicitComponent):
"""
Compute the advance ratio jze with a maximum 5.0.
"""

def initialize(self):
self.options.declare(
'num_nodes', types=int, default=1,
desc='Number of nodes to be evaluated in the RHS')
self.options.declare('smooth_zje', default=True, types=bool)
self.options.declare('alpha', default=100.0, types=float)

def setup(self):
nn = self.options['num_nodes']
range = np.arange(nn)
self.add_input("vktas", val=np.zeros(nn), units='knot')
self.add_input("tipspd", val=np.zeros(nn), units='ft/s')
self.add_input("sqa_array", val=np.zeros(nn), units='unitless')
self.add_output("equiv_adv_ratio", val=np.zeros(nn), units='unitless')

self.declare_partials("equiv_adv_ratio",
["vktas", "tipspd"],
rows=range, cols=range)

self.declare_partials("equiv_adv_ratio",
["sqa_array"],
rows=range, cols=range)

def compute(self, inputs, outputs):
nn = self.options['num_nodes']
vktas = inputs["vktas"]
tipspd = inputs["tipspd"]
sqa_array = inputs["sqa_array"]
equiv_adv_ratio = (1.0 - 0.254 * sqa_array) * 5.309 * vktas / tipspd

smooth = self.options["smooth_zje"]
if smooth:
alpha = self.options['alpha']
jze = smooth_min(equiv_adv_ratio, np.ones(nn) * 5.0, alpha)
else:
jze = np.minimum(equiv_adv_ratio, np.ones(nn) * 5.0)
outputs["equiv_adv_ratio"] = jze

def compute_partials(self, inputs, partials):
nn = self.options['num_nodes']
vktas = inputs["vktas"]
tipspd = inputs["tipspd"]
sqa_array = inputs["sqa_array"]
jze = (1.0 - 0.254 * sqa_array) * 5.309 * vktas / tipspd

djze_dsqa = -0.254 * 5.309 * vktas / tipspd
djze_dvktas = (1.0 - 0.254 * sqa_array) * 5.309 / tipspd
djze_dtipspd = -(1.0 - 0.254 * sqa_array) * 5.309 * vktas / tipspd**2

smooth = self.options["smooth_zje"]
if smooth:
alpha = self.options["alpha"]
djze_dsqa = d_smooth_min(jze, np.ones(nn) * 5.0, alpha) * djze_dsqa
djze_dvktas = d_smooth_min(jze, np.ones(nn) * 5.0, alpha) * djze_dvktas
djze_dtipspd = d_smooth_min(jze, np.ones(nn) * 5.0, alpha) * djze_dtipspd
else:
djze_dsqa = np.piecewise(jze, [jze < 5, jze >= 5], [1, 0]) * djze_dsqa
djze_dvktas = np.piecewise(jze, [jze < 5, jze >= 5], [1, 0]) * djze_dvktas
djze_dtipspd = np.piecewise(jze, [jze < 5, jze >= 5], [1, 0]) * djze_dtipspd
partials["equiv_adv_ratio", "sqa_array"] = djze_dsqa
partials["equiv_adv_ratio", "vktas"] = djze_dvktas
partials["equiv_adv_ratio", "tipspd"] = djze_dtipspd


class InstallLoss(om.Group):
"""
Compute installation loss
Expand All @@ -322,42 +486,18 @@ def setup(self):
nn = self.options['num_nodes']
self.add_subsystem(
name='sqa_comp',
subsys=om.ExecComp(
'sqa = minimum(DiamNac**2/DiamProp**2, 0.50)',
DiamNac={'val': 0, 'units': 'ft'},
DiamProp={'val': 0, 'units': 'ft'},
sqa={'units': 'unitless'},
has_diag_partials=True,
),
subsys=AreaSquareRatio(num_nodes=nn, smooth_sqa=True),
promotes_inputs=[("DiamNac", Aircraft.Nacelle.AVG_DIAMETER),
("DiamProp", Aircraft.Engine.PROPELLER_DIAMETER)],
promotes_outputs=["sqa"],
promotes_outputs=["sqa_array"],
)

# We should update these minimum calls to use a smooth minimum so that the
# gradient information is C1 continuous.
self.add_subsystem(
name='zje_comp', subsys=om.ExecComp(
'equiv_adv_ratio = minimum((1.0 - 0.254 * sqa) * 5.309 * vktas/tipspd, 5.0)',
vktas={'units': 'knot', 'val': np.zeros(nn)},
tipspd={'units': 'ft/s', 'val': np.zeros(nn)},
sqa={'units': 'unitless'},
equiv_adv_ratio={'units': 'unitless', 'val': np.zeros(nn)},
has_diag_partials=True,),
promotes_inputs=["sqa", ("vktas", Dynamic.Mission.VELOCITY),
name='zje_comp',
subsys=AdvanceRatio(num_nodes=nn, smooth_zje=True),
promotes_inputs=["sqa_array", ("vktas", Dynamic.Mission.VELOCITY),
("tipspd", Dynamic.Mission.PROPELLER_TIP_SPEED)],
promotes_outputs=["equiv_adv_ratio"],)

self.add_subsystem(
'convert_sqa',
om.ExecComp(
'sqa_array = sqa',
sqa={'units': 'unitless'},
sqa_array={'units': 'unitless', 'shape': (nn,)},
has_diag_partials=True,
),
promotes_inputs=["sqa"],
promotes_outputs=["sqa_array"],
promotes_outputs=["equiv_adv_ratio"],
)

self.blockage_factor_interp = self.add_subsystem(
Expand Down
Loading

0 comments on commit 990527d

Please sign in to comment.