Source code for gemseo.problems.sobieski.core.utils

# 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
"""Sobieski's SSBJ base class."""

from __future__ import annotations

import cmath
import math
from typing import TYPE_CHECKING

from numpy import array
from numpy import atleast_2d
from numpy import clip
from numpy import complex128
from numpy import concatenate
from numpy import float64
from numpy import ndarray
from strenum import StrEnum

if TYPE_CHECKING:
    from collections.abc import Sequence
    from types import ModuleType


DEG_TO_RAD = math.pi / 180.0
# DO NOT CHANGE THE DEFAULT DESIGN VALUE: IT IS USED FOR POLYNOMIAL APPROXIMATION
_DEFAULT_DESIGN = array([0.25, 1.0, 1.0, 0.5, 0.05, 45000.0, 1.6, 5.5, 55.0, 1000.0])
_MINIMUM_FUEL_WEIGHT = 2000.0
_MISC_WEIGHT = 25000.0
_MAXIMUM_LOAD_FACTOR = 6.0
_REFERENCE_ENGINE_WEIGHT = 4360.0
_MINIMUM_DRAG_COEFFICIENT = 0.01375
_CONSTANTS = array([
    _MINIMUM_FUEL_WEIGHT,
    _MISC_WEIGHT,
    _MAXIMUM_LOAD_FACTOR,
    _REFERENCE_ENGINE_WEIGHT,
    _MINIMUM_DRAG_COEFFICIENT,
])
_DESIGN_BOUNDS = array([
    (0.1, 0.4),
    (0.75, 1.25),
    (0.75, 1.25),
    (0.1, 1),
    (0.01, 0.09),
    (30000.0, 60000.0),
    (1.4, 1.8),
    (2.5, 8.5),
    (40.0, 70.0),
    (500.0, 1500.0),
])
_NAMES_TO_BOUNDS = {
    "x_1": _DESIGN_BOUNDS[0:2],
    "x_2": _DESIGN_BOUNDS[2],
    "x_3": _DESIGN_BOUNDS[3],
    "x_shared": _DESIGN_BOUNDS[4:10],
    "y_14": array([(4.97e04 * 0.5, 5.14e04 * 1.5), (-1.54e04 * 0.5, 3.0e04 * 1.5)]),
    "y_12": array([(4.97e04 * 0.5, 5.15e04 * 1.5), (0.9 * 0.5, 1.00)]),
    "y_21": array([(4.97e04 * 0.5, 5.15e04 * 1.5)]),
    "y_23": array([(6.73e03 * 0.5, 1.76e04 * 1.5)]),
    "y_24": array([(0.88 * 0.5, 7.42 * 1.5)]),
    "y_31": array([(5.92e03 * 0.5, 6.79e03 * 1.5)]),
    "y_32": array([(0.47 * 0.5, 0.53 * 1.5)]),
    "y_34": array([(0.88 * 0.5, 1.32 * 1.5)]),
}


[docs] class SobieskiBase: """Utilities for Sobieski's SSBJ use case.""" # TODO: API: remove these unused class attributes. DTYPE_COMPLEX = "complex128" DTYPE_DOUBLE = "float64" DTYPE_DEFAULT = DTYPE_COMPLEX
[docs] class DataType(StrEnum): """A NumPy data type.""" COMPLEX = "complex128" FLOAT = "float64"
dtype: DataType """The NumPy data type.""" math: ModuleType """The library of mathematical functions.""" def __init__( self, dtype: DataType, ) -> None: """ Args: dtype: The NumPy data type. """ # noqa: D205 D212 if dtype == self.DataType.COMPLEX: self.math = cmath self.dtype = complex128 elif dtype == self.DataType.FLOAT: self.math = math self.dtype = float64 else: msg = f"Unknown dtype: {dtype}." raise ValueError(msg) self.__constants = _CONSTANTS.astype(self.dtype) self.__coeff_mtrix = array( [ [0.2736, 0.3970, 0.8152, 0.9230, 0.1108], [0.4252, 0.4415, 0.6357, 0.7435, 0.1138], [0.0329, 0.8856, 0.8390, 0.3657, 0.0019], [0.0878, 0.7248, 0.1978, 0.0200, 0.0169], [0.8955, 0.4568, 0.8075, 0.9239, 0.2525], ], dtype=self.dtype, ) self.__mtx_shifted = array( [ [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], ], dtype=self.dtype, ) self.__f_bound = array( [[1.0], [1.0], [1.0]], dtype=self.dtype, ) self.__initial_design = _DEFAULT_DESIGN.astype(self.dtype) @property def initial_design(self) -> ndarray: """The initial design.""" return self.__initial_design @property def constants(self) -> ndarray: """The default constants.""" return self.__constants @property def design_bounds(self) -> tuple[ndarray, ndarray]: """The lower and upper bounds of the design variables.""" return _DESIGN_BOUNDS[:, 0], _DESIGN_BOUNDS[:, 1]
[docs] @classmethod def get_bounds_by_name( cls, variable_names: str | Sequence[str], ) -> tuple[ndarray, ndarray]: """Return the bounds of the design and coupling variables. Args: variable_names: The names of the variables. Returns: The lower and upper bounds of these variables. """ if isinstance(variable_names, str): variable_names = [variable_names] bounds = atleast_2d( concatenate( [_NAMES_TO_BOUNDS[variable_name] for variable_name in variable_names], axis=0, ) ) return bounds[:, 0], bounds[:, 1]
def __compute_mtx_shifted( self, s_bound: ndarray, index: int, ) -> None: """Compute a matrix of shifted values of the design variables. Args: s_bound: The bounds used to control the slope of the polynomial function. index: The index of the design variable in the polynomial function. """ self.__mtx_shifted[0, 1] = -s_bound[index] self.__mtx_shifted[2, 1] = s_bound[index] self.__mtx_shifted[0, 2] = self.__mtx_shifted[2, 2] = s_bound[index] ** 2 @staticmethod def __compute_a( mtx_shifted: ndarray, f_bound: ndarray, ao_coeff: ndarray, ai_coeff: ndarray, aij_coeff: ndarray, index: int, ) -> None: """Compute the interpolation terms. Args: mtx_shifted: The shift matrix. f_bound: The f bound. ao_coeff: The a0 term. ai_coeff: The ai terms aij_coeff: The aij term index: The index. """ ao_coeff[:] = f_bound[1] ai_coeff[index] = -(f_bound[2] - f_bound[0]) / (2 * mtx_shifted[0, 1]) aij_coeff[index, index] = ( f_bound[0] - f_bound[1] + (f_bound[2] - f_bound[0]) / 2 ) / mtx_shifted[0, 2] def __update_aij( self, aij_coeff: ndarray, imax: int, ) -> ndarray: """Update the quadratic interpolation terms. Args: aij_coeff: The quadratic terms. imax: The size of the interpolation matrix. Returns: The modified quadratic terms. """ for i in range(imax): for j in range(i + 1, imax): aij_coeff[i, j] = aij_coeff[i, i] * self.__coeff_mtrix[i, j] aij_coeff[j, i] = aij_coeff[i, j] return aij_coeff def __compute_fbound( self, flag: ndarray, s_shifted: ndarray, a_coeff: ndarray, b_coeff: ndarray, index: int, ) -> ndarray: """Compute right-hand side of polynomial function system. Args: flag: The functional relationship between variables. s_shifted: The normalized values of the independent variables shifted around the origin. a_coeff: The a coefficient. b_coeff: The b coefficient. index: The index. Returns: The right-hand side of polynomial function system. """ if flag[index] == 5: self.__f_bound[0, 0] = 1.0 + 0.25 * a_coeff * a_coeff self.__f_bound[2, 0] = 1.0 + 0.25 * b_coeff * b_coeff else: if flag[index] == 0: s_shifted = 0.0 elif flag[index] == 3: a_coeff = -a_coeff b_coeff = a_coeff elif flag[index] == 2: b_coeff = 2 * a_coeff elif flag[index] == 4: a_coeff = -a_coeff b_coeff = 2 * a_coeff self.__f_bound[0, 0] = 1.0 - 0.5 * a_coeff self.__f_bound[2, 0] = 1.0 + 0.5 * b_coeff return s_shifted @staticmethod def _compute_normalization( s_ref: ndarray, s_new: ndarray, ) -> ndarray: """Normalization of input variables for polynomial approximation. Args: s_ref: The initial values of the independent variables (5 variables at max). s_new: The current values of the independent variables. Returns: The normalized and centered value. """ return clip(s_new / s_ref, 0.75, 1.25) - 1.0
[docs] @staticmethod def derive_normalization( s_ref: float | ndarray, s_new: float | ndarray, ) -> float | ndarray: """Derivation of normalization of input variables. For use of polynomial approximation. Args: s_ref: The initial values of the independent variables (5 variables at max). s_new: The current values of the independent variables. Returns: Derivatives of normalized value. """ s_norm = s_new.real / s_ref.real if 0.75 <= s_norm <= 1.25: return 1.0 / s_ref return 0.0
[docs] def derive_polynomial_approximation( self, s_ref: ndarray, s_new: ndarray, flag: int, s_bound: ndarray, a0_coeff: ndarray, ai_coeff: ndarray, aij_coeff: ndarray, ) -> tuple[ndarray, ndarray, ndarray, ndarray]: """Compute the polynomial coefficients for both evaluation and linearization. These coefficients characterize the behavior of certain synthetic variables and function modifiers. Args: s_ref: The initial values of the independent variables (5 variables at max). s_new: The current values of the independent variables. flag: The functional relationship between the variables: - flag = 1: linear >0, - flag = 2: nonlinear >0, - flag = 3: linear < 0, - flag = 4: nonlinear <0, - flag = 5: parabolic. s_bound: The offset value for normalization. Returns: The value of the synthetic variables or function modifiers for both evaluation and linearization. """ imax = s_ref.shape[0] # Normalize new S with initial s_shifted = self._compute_normalization(s_ref, s_new) a0_coeff[...] = 0 ai_coeff[...] = 0 aij_coeff[...] = 0 for i in range(imax): a_coeff = 0.1 b_coeff = a_coeff self.__compute_mtx_shifted(s_bound, i) s_shifted = self.__compute_fbound(flag, s_shifted, a_coeff, b_coeff, i) self.__compute_a( self.__mtx_shifted, self.__f_bound, a0_coeff, ai_coeff, aij_coeff, i ) aij_coeff = self.__update_aij(aij_coeff, imax) poly_value = ( a0_coeff + ai_coeff @ s_shifted + 0.5 * s_shifted.T @ (aij_coeff[:imax, :imax] @ s_shifted) ) return poly_value[0], ai_coeff, aij_coeff[:imax, :imax], s_shifted
[docs] def compute_polynomial_approximation( self, s_ref: ndarray, s_new: ndarray, flag: int, s_bound: ndarray, a0_coeff: ndarray, ai_coeff: ndarray, aij_coeff: ndarray, ) -> ndarray: """Compute the polynomial coefficients. These coefficients characterize the behavior of certain synthetic variables and function modifiers. Args: s_ref: The initial values of the independent variables (5 variables at max). s_new: The current values of the independent variables. flag: The functional relationship between the variables: - flag = 1: linear >0, - flag = 2: nonlinear >0, - flag = 3: linear < 0, - flag = 4: nonlinear <0, - flag = 5: parabolic. s_bound: The offset value for normalization. Returns: The value of the synthetic variables or function modifiers. """ imax = s_ref.shape[0] # Normalize new S with initial s_shifted = self._compute_normalization(s_ref, s_new) # - ai_coeff: vector of coeff for 2nd term # - aij_coeff: matrix of coeff for 3rd term # - ao_coeff: scalar coeff a0_coeff[...] = 0 ai_coeff[...] = 0 aij_coeff[...] = 0 for i in range(imax): a_coeff = 0.1 b_coeff = a_coeff self.__compute_mtx_shifted(s_bound, i) s_shifted = self.__compute_fbound(flag, s_shifted, a_coeff, b_coeff, i) self.__compute_a( self.__mtx_shifted, self.__f_bound, a0_coeff, ai_coeff, aij_coeff, i ) aij_coeff = self.__update_aij(aij_coeff, imax) poly_value = ( a0_coeff + ai_coeff @ s_shifted + 0.5 * s_shifted.T @ (aij_coeff[:imax, :imax] @ s_shifted) ) return poly_value[0]
[docs] def compute_half_span( self, aspect_ratio: float, wing_surface_area: float, ) -> float: """Compute the half-span from the wing surface and aspect ratio. Args: aspect_ratio: The aspect ratio. wing_surface_area: The wing_surface_area. Returns: The half-span. """ return self.math.sqrt(aspect_ratio * wing_surface_area) * 0.5
[docs] def compute_thickness( self, aspect_ratio: float, thickness_to_chord_ratio: float, wing_surface_area: float, ) -> float: """Compute a wing thickness. Args: aspect_ratio: The aspect ratio. thickness_to_chord_ratio: The thickness to chord ratio. wing_surface_area: The wing_surface_area. Returns: The wing thickness. """ return ( thickness_to_chord_ratio * wing_surface_area / (self.math.sqrt(aspect_ratio * wing_surface_area)) )
[docs] @staticmethod def compute_aero_center( wing_taper_ratio: float, ) -> float: """Computes the aerodynamic center. Args: wing_taper_ratio: The wing taper ratio. Returns: The aerodynamic center. """ return (1.0 + 2.0 * wing_taper_ratio) / (3.0 * (1 + wing_taper_ratio))
[docs] def get_initial_values(self) -> tuple[float]: """Return the initial values used by the polynomial functions. Returns: The initial values used by the polynomial functions. """ x_initial = self.initial_design[1] tc_initial = self.initial_design[4] half_span_initial = ( self.math.sqrt(self.initial_design[7] * self.initial_design[9]) * 0.5 ) aero_center_initial = (1 + 2 * self.initial_design[0]) / ( 3 * (1 + self.initial_design[0]) ) cf_initial = self.initial_design[2] mach_initial = self.initial_design[6] h_initial = self.initial_design[5] throttle_initial = self.initial_design[3] lift_initial = self.dtype(1.0) twist_initial = self.dtype(1.0) esf_initial = self.dtype(1.0) return ( x_initial, tc_initial, half_span_initial, aero_center_initial, cf_initial, mach_initial, h_initial, throttle_initial, lift_initial, twist_initial, esf_initial, )