# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# Contributors:
# INITIAL AUTHORS - initial API and implementation
# and/or initial documentation
# :author: Sobieski, Agte, and Sandusky
# OTHER AUTHORS - MACROSCOPIC CHANGES
# :author: Damien Guenot
# :author: Francois Gallard
# From NASA/TM-1998-208715
# Bi-Level Integrated System Synthesis (BLISS)
# Sobieski, Agte, and Sandusky
"""
SSBJ Propulsion computations
****************************
"""
from __future__ import absolute_import, division, print_function, unicode_literals
from math import pi
from future import standard_library
from numpy import append, array, atleast_2d, zeros
standard_library.install_aliases()
from gemseo import LOGGER
DEG_TO_RAD = pi / 180.0
[docs]class SobieskiPropulsion(object):
"""Class defining propulsion analysis for Sobieski problem and
related method to the propulsion problem
such as disciplines computation, constraints, reference optimum
"""
DTYPE_COMPLEX = "complex128"
DTYPE_DOUBLE = "float64"
ESF_UPPER_LIMIT = 1.5
ESF_LOWER_LIMIT = 0.5
TEMPERATURE_LIMIT = 1.02
def __init__(self, sobieski_base):
"""
Constructor
"""
self.base = sobieski_base
self.constants = self.base.default_constants()
(
self.x_initial,
self.tc_initial,
self.half_span_initial,
self.aero_center_initial,
self.cf_initial,
self.mach_initial,
self.h_initial,
self.throttle_initial,
self.lift_initial,
self.twist_initial,
self.esf_initial,
) = self.base.get_initial_values()
self.dtype = self.base.dtype
self.math = self.base.math
# Surface fit to engine deck with least square method
# Polynomial coefficients for SFC computation
self.sfc_coeff = array(
[
1.13238425638512,
1.53436586044561,
-0.00003295564466,
-0.00016378694115,
-0.31623315541888,
0.00000410691343,
-0.00005248000590,
-0.00000000008574,
0.00000000190214,
0.00000001059951,
],
dtype=self.dtype,
)
# Polynomial coefficients for throttle constraint
self.thua_coeff = array(
[
11483.7822254806,
10856.2163466548,
-0.5080237941,
3200.157926969,
-0.1466251679,
0.0000068572,
],
dtype=self.dtype,
)
self.throttle_coeff = self.dtype(16168.6)
def __set_coeff_temp(self, x_shared, x_3):
"""
Prepare settings of polynomial function
for temperature
:param x_shared: global design variables
:type x_shared: numpy array
:param x_3: local design variables
:type x_3: numpy array
:returns: s_initial, s_new, flag, bound
:rtype: numpy array
"""
s_initial = array(
[self.mach_initial, self.h_initial, self.throttle_initial], dtype=self.dtype
)
s_new = array([x_shared[2], x_shared[1], x_3[0]], dtype=self.dtype)
flag = array([2, 4, 2], dtype=self.dtype)
bound = array([0.25, 0.25, 0.25], dtype=self.dtype)
return s_initial, s_new, flag, bound
[docs] def compute_dim_throttle(self, adim_throttle):
"""Compute a dimensioned value of throttle from adim value
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: throttle
:rtype: numpy array
"""
return adim_throttle * self.throttle_coeff
[docs] def compute_sfc(self, x_shared, adim_throttle):
"""Compute Specific Fuel Consumption (SFC) from global design variables
(M, h,) and local design variable (throttle) by polynomial function
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: sfc
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
sfc = (
self.sfc_coeff[0]
+ self.sfc_coeff[1] * x_shared[2]
+ self.sfc_coeff[2] * x_shared[1]
+ self.sfc_coeff[3] * throttle
+ self.sfc_coeff[4] * x_shared[2] ** 2
+ 2 * x_shared[1] * x_shared[2] * self.sfc_coeff[5]
+ 2 * throttle * x_shared[2] * self.sfc_coeff[6]
+ self.sfc_coeff[7] * x_shared[1] ** 2
+ 2 * throttle * x_shared[1] * self.sfc_coeff[8]
+ self.sfc_coeff[9] * throttle ** 2
)
return sfc
[docs] def compute_throttle_ua(self, x_shared):
"""Compute throttle upper limit from global design variables
(M, h,)
:param x_shared: global design vector
:type x_shared: numpy array
:returns: throttle_uA
:rtype: numpy array
"""
throttle_ua = (
self.thua_coeff[0]
+ self.thua_coeff[1] * x_shared[2]
+ self.thua_coeff[2] * x_shared[1]
+ self.thua_coeff[3] * x_shared[2] * x_shared[2]
+ 2 * self.thua_coeff[4] * x_shared[2] * x_shared[1]
+ self.thua_coeff[5] * x_shared[1] * x_shared[1]
)
return throttle_ua
[docs] def compute_throttle_constraint(self, x_shared, adim_throttle):
"""Compute throttle constraint
(M, h,) and local design variable (throttle) by polynomial function
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: throttle constraint as throttle / throttle_ua - 1.0
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
throttle_ua = self.compute_throttle_ua(x_shared)
throttle_constraint = throttle / throttle_ua - 1.0 # throttle setting
return throttle_constraint
[docs] def compute_dthrconst_dthrottle(self, x_shared):
"""Compute derivative of throttle constraint wrt throttle
:param x_shared: global design vector
:type x_shared: numpy array
:returns: dthrottle_constraint_dthrottle
:rtype: numpy array
"""
throttle_ua = self.compute_throttle_ua(x_shared)
return self.throttle_coeff / throttle_ua
[docs] def compute_dthrcons_dh(self, x_shared, adim_throttle):
"""Compute derivative of throttle constraint wrt altitude
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: dthrottle_constraint_dh
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
throttle_ua = self.compute_throttle_ua(x_shared)
dthrottle_ua_dh = (
self.thua_coeff[2]
+ 2 * self.thua_coeff[4] * x_shared[2]
+ 2.0 * self.thua_coeff[5] * x_shared[1]
)
return -throttle * dthrottle_ua_dh / (throttle_ua * throttle_ua)
[docs] def compute_dthrconst_dmach(self, x_shared, adim_throttle):
"""Compute derivative of throttle constraint wrt Mach number
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: dthrottle_constraint_dmach
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
throttle_ua = self.compute_throttle_ua(x_shared)
dthrottle_ua_dmach = (
self.thua_coeff[1]
+ 2 * self.thua_coeff[3] * x_shared[2]
+ 2.0 * self.thua_coeff[4] * x_shared[1]
)
return -throttle * dthrottle_ua_dmach / (throttle_ua * throttle_ua)
[docs] def compute_esf(self, drag, adim_throttle):
"""Compute engine scale factor
:param drag: drag
:type drag: float
:param adim_throttle: throttle (not dimensioned)
:type adim_throttle: float
:returns: Engine Scale Factor
:rtype: numpy array
"""
return drag / (3.0 * self.compute_dim_throttle(adim_throttle))
[docs] def compute_desf_ddrag(self, adim_throttle):
"""Compute derivative of ESF wrt aero drag
:param adim_throttle: local design variables
:type adim_throttle: numpy array
:returns: derivative of ESF wrt drag
:rtype: numpy array
"""
return 1.0 / (3 * self.compute_dim_throttle(adim_throttle))
[docs] def compute_desf_dthrottle(self, drag, adim_throttle):
"""Compute derivative of ESF wrt aero drag
:param drag: coupling design variables
:type drag: numpy array
:param adim_throttle: local design variables
:type adim_throttle: numpy array
:returns: derivative of ESF wrt drag
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
return -self.throttle_coeff * drag / (3.0 * throttle ** 2)
[docs] def compute_temp(self, x_shared, x_3):
"""Compute engine temperature
:param x_3: local design vector
:type x_3: numpy array
:param x_shared: global design vector
:type x_shared: numpy array
:returns: engine temperature
:rtype: numpy array
"""
s_initial, s_new, flag, bound = self.__set_coeff_temp(x_shared, x_3)
return self.base.poly_approx(s_initial, s_new, flag, bound)
[docs] def compute_engine_weight(self, esf):
"""Compute engine weight
:param esf: engine scale factor
:returns: engine weight
:rtype: numpy array
"""
return self.constants[3] * (esf ** 1.05) * 3
[docs] def blackbox_propulsion(self, x_shared, y_23, x_3, true_cstr=False):
"""This function calculates fuel comsumption, engine weight and
engine scale factor
:param x_shared: shared design variable vector:
- x_shared[0]: thickness/chord ratio
- x_shared[1]: altitude
- x_shared[2]: Mach
- x_shared[3]: aspect ratio
- x_shared[4]: wing sweep
- x_shared[5]: wing surface area
:type x_shared: numpy array
:param y_23: shared variables coming from blackbox_aerodynamics (drag)
:type y_23: numpy array
:param x_3: power/propulsion design variable (throttle setting)
:type x_3: numpy array
:param true_cstr: Default value = False)
:returns: y_3, y_34, y_31, y_32, g_3:
- y_3: output variables for propulsion analysis
- y_3[0]: SFC
- y_3[1]: engine weight
- y_3[2]: engine scale factor
- y_34: shared variable for BlackBoxMission (SFC)
- y_31: shared variable for blackbox_structure (engine weight)
- y_32: shared variable for blackbox_aerodynamics (esf)
- g_3: propulsion constraints
- g_3[0]: engine scale factor constraint
- g_3[1]: engine temperature
- g_3[2]: throttle setting constraint
:rtype: numpy array, numpy array, numpy array, numpy array, numpy array
"""
y_3 = zeros((3), dtype=self.dtype)
g_3 = zeros((3), dtype=self.dtype)
y_31 = zeros((1), dtype=self.dtype)
y_32 = zeros((1), dtype=self.dtype)
y_34 = zeros((1), dtype=self.dtype)
esf = self.compute_esf(y_23[0], x_3[0])
y_3[2] = esf
y_3[1] = self.compute_engine_weight(esf)
y_3[0] = self.compute_sfc(x_shared, x_3[0])
y_31[0] = y_3[1]
y_34[0] = y_3[0]
y_32[0] = y_3[2]
# THIS SECTION COMPUTES SFC, esf, AND ENGINE WEIGHT
g_3[0] = y_3[2] # engine scale factor
# engine temperature
temp = self.compute_temp(x_shared, x_3)
g_3[1] = temp
g_3[2] = self.compute_throttle_constraint(x_shared, x_3[0])
if not true_cstr:
g_3 = append(
g_3[0] - self.ESF_UPPER_LIMIT,
(
self.ESF_LOWER_LIMIT - g_3[0],
g_3[2],
g_3[1] - self.TEMPERATURE_LIMIT,
),
)
return y_3, y_34, y_31, y_32, g_3
[docs] def compute_dengineweight_dvar(self, esf, desf_dx):
"""Computes derivative of engine weight wrt to
a variable(drag or throttle)
:param esf: ESF
:type esf: float
:param desf_dx: partial derivative of ESF wrt to input variable
:type desf_dx: numpy array
:returns: derivative of engine weight wrt a variable (drag or throttle)
:rtype: numpy array
"""
return 3 * self.constants[3] * 1.05 * desf_dx * esf ** 0.05
[docs] def compute_dsfc_dthrottle(self, x_shared, adim_throttle):
"""Compute derivative of sfc constraint wrt throttle
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: dsfc_dthrottle
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
dsfc_dthrottle = (
self.sfc_coeff[3]
+ 2 * x_shared[2] * self.sfc_coeff[6]
+ 2 * x_shared[1] * self.sfc_coeff[8]
+ 2 * self.sfc_coeff[9] * throttle
)
dsfc_dthrottle *= self.throttle_coeff
return dsfc_dthrottle
[docs] def compute_dsfc_dh(self, x_shared, adim_throttle):
"""Compute derivative of sfc constraint wrt altitude
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: dsfc_dh
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
dsfc_dh = (
self.sfc_coeff[2]
+ 2 * x_shared[2] * self.sfc_coeff[5]
+ +2 * self.sfc_coeff[7] * x_shared[1]
+ 2 * throttle * self.sfc_coeff[8]
)
return dsfc_dh
[docs] def compute_dsfc_dmach(self, x_shared, adim_throttle):
"""Compute derivative of sfc constraint wrt Mach number
:param x_shared: global design vector
:type x_shared: numpy array
:param adim_throttle: local design vector
:type adim_throttle: numpy array
:returns: dsfc_dmach
:rtype: numpy array
"""
throttle = self.compute_dim_throttle(adim_throttle)
dsfc_dmach = (
self.sfc_coeff[1]
+ 2 * self.sfc_coeff[4] * x_shared[2]
+ 2 * x_shared[1] * self.sfc_coeff[5]
+ 2 * throttle * self.sfc_coeff[6]
)
return dsfc_dmach
def __dadimthrottle_dthrottle(self, x_3):
"""
Compute partial derivative of adim throttle of polynomial
function wrt throttle
:param x_3: local design variables
:type x_3: numpy array
:returns: dadimthrottle_dthrottle
:rtype: numpy array
"""
return self.base.derive_normalize_s(self.throttle_initial, x_3[0])
def __compute_dadimh_dh(self, x_shared):
"""
Compute partial derivative of adim throttle of polynomial
function wrt altitude
:param x_shared: global design variables
:type x_shared: numpy array
:returns: dadimh_dh
:rtype: numpy array
"""
return self.base.derive_normalize_s(self.h_initial, x_shared[1])
def __compute_dadimmach_dmach(self, x_shared):
"""
Compute partial derivative of adim throttle of polynomial
function wrt Mach number
:param x_shared: global design variables
:type x_shared: numpy array
:returns: dadimthrottle_dthrottle
:rtype: numpy array
"""
return self.base.derive_normalize_s(self.mach_initial, x_shared[2])
def __initialize_jacobian(self, true_cstr):
"""
Initialization of jacobian matrix
:param true_cstr:
:type true_cstr: logical
:returns: jacobian
:rtype: dict of dict of numpy array
"""
# Jacobian matrix as a dictionary
jacobian = {}
jacobian["y_3"] = {}
jacobian["g_3"] = {}
jacobian["y_31"] = {}
jacobian["y_32"] = {}
jacobian["y_34"] = {}
jacobian["y_3"]["x_3"] = zeros((3, 1), dtype=self.dtype)
jacobian["y_3"]["x_shared"] = zeros((3, 6), dtype=self.dtype)
jacobian["y_3"]["y_23"] = zeros((3, 1), dtype=self.dtype)
if not true_cstr:
n_constraints = 4
else:
n_constraints = 3
jacobian["g_3"]["x_3"] = zeros((n_constraints, 1), dtype=self.dtype)
jacobian["g_3"]["x_shared"] = zeros((n_constraints, 6), dtype=self.dtype)
jacobian["g_3"]["y_23"] = zeros((n_constraints, 1), dtype=self.dtype)
return jacobian
[docs] def derive_blackbox_propulsion(self, x_shared, y_23, x_3, true_cstr=False):
"""Compute jacobian matrix of propulsion analysis
:param x_shared: shared design variable vector:
- x_shared[0]: thickness/chord ratio
- x_shared[1]: altitude
- x_shared[2]: Mach
- x_shared[3]: aspect ratio
- x_shared[4]: wing sweep
- x_shared[5]: wing surface area
:type x_shared: numpy array
:param y_23: shared variables coming from blackbox_aerodynamics (drag)
:type y_23: numpy array
:param x_3: power/propulsion design variable (throttle setting)
:type x_3: numpy array
:param true_cstr: Default value = False)
:returns: jacobian : Jacobian matrix
:rtype: dict(dict(ndarray))
"""
# Jacobian matrix as a dictionary
jacobian = self.__initialize_jacobian(true_cstr)
dg_3_dx_3 = zeros((3, 1), dtype=self.dtype)
dg_3_dxs = zeros((3, 6), dtype=self.dtype)
dg_3_dy_23 = zeros((3, 1), dtype=self.dtype)
drag = y_23[0]
esf = self.compute_esf(y_23[0], x_3[0])
# dSFC_dthrottle
jacobian["y_3"]["x_3"][0, 0] = self.compute_dsfc_dthrottle(x_shared, x_3[0])
# dESF_dthrottle
jacobian["y_3"]["x_3"][2, 0] = self.compute_desf_dthrottle(drag, x_3[0])
# dengineweight_dthrottle
jacobian["y_3"]["x_3"][1, :] = self.compute_dengineweight_dvar(
esf, jacobian["y_3"]["x_3"][2, 0]
)
# dSFC_d(t/c) = 0
# dESF_d(t/c) = 0
# dengineweight_d(t/c) = 0
# dSFC_dh
jacobian["y_3"]["x_shared"][0, 1] = self.compute_dsfc_dh(x_shared, x_3[0])
# dESF_dh= 0.0
# dengineweight_dh= 0.0
# dSFC_dM
jacobian["y_3"]["x_shared"][0, 2] = self.compute_dsfc_dmach(x_shared, x_3[0])
# dESF_dM= 0.0
# dengineweight_dM= 0.0
# jacobian['y_3']['x_shared'][:, 2:] = 0.0
# dSFC_ddrag
# jacobian['y_3']['y_23'][0, 0] = 0.0
# dESF_ddrag
jacobian["y_3"]["y_23"][2, 0] = self.compute_desf_ddrag(x_3[0])
# dengineweight_ddrag
jacobian["y_3"]["y_23"][1, :] = self.compute_dengineweight_dvar(
esf, jacobian["y_3"]["y_23"][2, 0]
)
# dtemp_ddrag
# jacobian['g_3']['y_23'][0, 0] = 0.0
s_initial, s_new, flag, bound = self.__set_coeff_temp(x_shared, x_3)
_, ai_coeff, aij_coeff, s_shifted = self.base.derive_poly_approx(
s_initial, s_new, flag, bound
)
# dadimtemp_dtemp = self.__compute_dadimtemp_dtemp(x_shared)
# g_3[0, :] = ESF
dg_3_dx_3[0, :] = jacobian["y_3"]["x_3"][2, :]
# dtemp_dthrottle
dg_3_dx_3[1, 0] = self.__dadimthrottle_dthrottle(x_3) * (
ai_coeff[2]
+ aij_coeff[2, 0] * s_shifted[0]
+ aij_coeff[2, 1] * s_shifted[1]
+ aij_coeff[2, 2] * s_shifted[2]
)
# d(throttle-throttle_ua)_dthrottle
dg_3_dx_3[2, 0] = self.compute_dthrconst_dthrottle(x_shared)
# g_3[0, :] = ESF
dg_3_dxs[0, :] = jacobian["y_3"]["x_shared"][2, :]
# dtemp_d(t/c)= 0.0
# d(throttle-throttle_ua)_d(t/c)= 0.0
# dtemp_dh
dg_3_dxs[1, 1] = self.__compute_dadimh_dh(x_shared) * (
ai_coeff[1]
+ aij_coeff[1, 0] * s_shifted[0]
+ aij_coeff[1, 1] * s_shifted[1]
+ aij_coeff[1, 2] * s_shifted[2]
)
# d(throttle-throttle_ua)_dh
dg_3_dxs[2, 1] = self.compute_dthrcons_dh(x_shared, x_3[0])
# dtemp_dM
dg_3_dxs[1, 2] = self.__compute_dadimmach_dmach(x_shared) * (
ai_coeff[0]
+ aij_coeff[0, 0] * s_shifted[0]
+ aij_coeff[0, 1] * s_shifted[1]
+ aij_coeff[0, 2] * s_shifted[2]
)
# d(throttle-throttle_ua)_dM
dg_3_dxs[2, 2] = self.compute_dthrconst_dmach(x_shared, x_3[0])
jacobian = self.__set_coupling_jacobian(jacobian)
# THIS SECTION COMPUTES SFC, ESF, AND ENGINE WEIGHT
dg_3_dxs[0, :] = jacobian["y_3"]["x_shared"][2, :]
dg_3_dy_23[0, :] = jacobian["y_3"]["y_23"][2, :]
if true_cstr:
jacobian["g_3"]["x_3"] = dg_3_dx_3
jacobian["g_3"]["x_shared"] = dg_3_dxs
jacobian["g_3"]["y_23"] = dg_3_dy_23
else:
for i_g_jac, i_orig in enumerate([0, 0, 2, 1]):
jacobian["g_3"]["x_shared"][i_g_jac, :] = dg_3_dxs[i_orig, :]
jacobian["g_3"]["y_23"][i_g_jac, :] = dg_3_dy_23[i_orig, :]
jacobian["g_3"]["x_3"][i_g_jac, :] = dg_3_dx_3[i_orig, :]
if i_g_jac == 1:
jacobian["g_3"]["x_shared"][i_g_jac, :] *= -1
jacobian["g_3"]["y_23"][i_g_jac, :] *= -1
jacobian["g_3"]["x_3"][i_g_jac, :] *= -1
return jacobian
@staticmethod
def __set_coupling_jacobian(jacobian):
"""
Set jacobian of coupling variables
"""
jacobian["y_31"]["x_3"] = atleast_2d(jacobian["y_3"]["x_3"][1, :])
jacobian["y_31"]["x_shared"] = atleast_2d(jacobian["y_3"]["x_shared"][1, :])
jacobian["y_31"]["y_23"] = atleast_2d(jacobian["y_3"]["y_23"][1, :])
jacobian["y_32"]["x_3"] = atleast_2d(jacobian["y_3"]["x_3"][2, :])
jacobian["y_32"]["x_shared"] = atleast_2d(jacobian["y_3"]["x_shared"][2, :])
jacobian["y_32"]["y_23"] = atleast_2d(jacobian["y_3"]["y_23"][2, :])
jacobian["y_34"]["x_3"] = atleast_2d(jacobian["y_3"]["x_3"][0, :])
jacobian["y_34"]["x_shared"] = atleast_2d(jacobian["y_3"]["x_shared"][0, :])
jacobian["y_34"]["y_23"] = atleast_2d(jacobian["y_3"]["y_23"][0, :])
return jacobian