Source code for gemseo.problems.sobieski.base

# -*- 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 base class
***************
"""
from __future__ import division, unicode_literals

import cmath
import logging
import math

from numpy import array, asarray, complex128, concatenate, dot, float64, zeros
from numpy.linalg import lstsq
from scipy import linalg
from six import string_types

LOGGER = logging.getLogger(__name__)


DEG_TO_RAD = math.pi / 180.0


[docs]class SobieskiBase(object): """Class defining Sobieski problem and related method to the problem such as disciplines computation, constraints, reference optimum.""" DTYPE_COMPLEX = "complex128" DTYPE_DOUBLE = "float64" DTYPE_DEFAULT = DTYPE_COMPLEX def __init__(self, dtype): """Constructor. :param dtype: data type :type dtype: str """ self.dtype = dtype if dtype == complex128: self.math = cmath elif dtype == float64: self.math = math elif dtype == self.DTYPE_DOUBLE: self.math = math self.dtype = float64 elif dtype == self.DTYPE_COMPLEX: self.math = cmath self.dtype = complex128 else: raise ValueError("Unknown dtype : " + str(dtype)) self.i_0 = None
[docs] def get_default_x0(self): """Return a default initial value for design variables. :returns: initial design variables :rtype: numpy array warning: ##DO NOT CHANGE VALUE: THEY ARE USED FOR POLYNOMIAL APPROXIMATION ## """ return array( [0.25, 1.0, 1.0, 0.5, 0.05, 45000.0, 1.6, 5.5, 55.0, 1000.0], dtype=self.dtype, )
[docs] def default_constants(self): """Definition of constants vector constants for Sobieski problem. :returns: constant vector :rtype: numpy array """ constants = zeros(5, dtype=self.dtype) # Constants of problem constants[0] = 2000.0 # minimum fuel weight constants[1] = 25000.0 # miscellaneous weight constants[2] = 6.0 # Maximum load factor constants[3] = 4360.0 # Engine weight reference constants[4] = 0.01375 # Minimum drag coefficient return constants
[docs] @staticmethod def get_sobieski_bounds_tuple(): """Set the input design bounds and return them as a tuple of tuples. :returns: Subtuple is build with lower bound and upper bound. :rtype: tuple(tuple) """ bounds_tuple = ( (0.1, 0.4), (0.75, 1.25), (0.75, 1.25), (0.1, 1), (0.01, 0.09), (30000.0, 60000.0), # (0.01, 0.09) with threshold, (0.038, 0.06) # without threshold (1.4, 1.8), (2.5, 8.5), (40.0, 70.0), (500.0, 1500.0), ) return bounds_tuple
[docs] @classmethod def get_sobieski_bounds(cls): """Set the input design bounds and return them as 2 numpy arrays. :returns: upper and lower bounds :rtype: numpy array, numpy array """ bounds_tuple = cls.get_sobieski_bounds_tuple() bounds_array = asarray(bounds_tuple) upper_bounds = bounds_array[:, 1] lower_bounds = bounds_array[:, 0] return upper_bounds, lower_bounds
[docs] @classmethod def get_bounds_by_name(cls, variables_names): """Return bounds of design variables and coupling variables. :param variables_names: names of variables :type variables_names: str or list(str) :returns: lower bound and upper bound """ if isinstance(variables_names, string_types): variables_names = [variables_names] bounds_tuple = cls.get_sobieski_bounds_tuple() bounds_dict = { "x_1": array(bounds_tuple[0:2]), "x_2": array([bounds_tuple[2]]), "x_3": array([bounds_tuple[3]]), "x_shared": array(bounds_tuple[4:10]), "y_14": array([(4.97e04, 5.14e04), (-1.54e04, 3.0e04)]), "y_12": array([(4.97e04, 5.15e04), (0.9, 1.00)]), "y_21": array([(4.97e04, 5.15e04)]), "y_23": array([(6.73e03, 1.76e04)]), "y_24": array([(0.88, 7.42)]), "y_31": array([(5.92e03, 6.79e03)]), "y_32": array([(0.47, 0.53)]), "y_34": array([(0.88, 1.32)]), } for yuk in ["y_21", "y_14", "y_12", "y_23", "y_24", "y_31", "y_32", "y_34"]: bounds_dict[yuk][:, 0] = bounds_dict[yuk][:, 0] * 0.5 bounds_dict[yuk][:, 1] = bounds_dict[yuk][:, 1] * 1.5 bounds = bounds_dict[variables_names[0]] # LOGGER.debug(15 * " " + "get_bounds_by_name") for var_name in variables_names[1:]: # LOGGER.debug("var_name " + var_name) bounds = concatenate((bounds, bounds_dict[var_name]), axis=0) return bounds[:, 0], bounds[:, 1]
def __compute_mtx_shifted(self, s_bound, index): """Compute a matrix of shifted values of design variables. :param s_bound: vector of bounds used to control slope of polynomial function :type s_bound: numpy array :param index: index of design variable in polynomial function :type index: integer :returns: mtx_shifted :rtype: numpy array """ # s_mid = 0.0 # independent variable mid point # s_lower = s_mid - s_bound[index] # s_upper = s_mid + s_bound[index] s_lower = -s_bound[index] s_upper = s_bound[index] s_lower_2 = s_lower * s_lower # So_2 = s_mid * s_mid s_upper_2 = s_upper * s_upper mtx_shifted = array( [ [1.0, s_lower, s_lower_2], # [1.0, s_mid, s_mid_2], [1.0, 0.0, 0.0], [1.0, s_upper, s_upper_2], ], dtype=self.dtype, ) return mtx_shifted @staticmethod def __compute_a(mtx_shifted, f_bound, ao_coeff, ai_coeff, aij_coeff, index): """Compute the interpolation terms. :returns: ao_coeff, ai_coeff, aij_coeff constant, linear and quadratic terms :rtype: float, 1D array, 2D array (numpy) """ try: a_mat = linalg.solve(mtx_shifted, f_bound) except linalg.LinAlgError: LOGGER.warning( "Sobieski polynomial approximation: " "exact linear system solve failed, using " " approximate least squares method instead." ) a_mat, _, _, _ = lstsq(mtx_shifted, f_bound) ao_coeff = a_mat[0] ai_coeff[index] = a_mat[1][0] aij_coeff[index, index] = a_mat[2][0] return ao_coeff, ai_coeff, aij_coeff def __update_aij(self, aij_coeff, imax): """Update of quadratic interpolation terms. :param aij_coeff: array of quadratic terms :type aij_coeff: numpy array :param imax: size of interpolation matrix :type imax: integer :returns: aij_coeff (modified quadratic terms) :rtype: numpy array """ 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, ) for i in range(imax): for j in range(i + 1, imax): aij_coeff[i, j] = aij_coeff[i, i] * coeff_mtrix[i, j] aij_coeff[j, i] = aij_coeff[i, j] return aij_coeff def __compute_fbound(self, flag, s_shifted, a_coeff, b_coeff, index): """Compute right-hand side of polynomial function system. :param flag: functional relationship between var :type flag: numpy array :param s_shifted: vector of normalized values of independent variables shifted around origin :type s_shifted: numpy array """ if flag[index] == 5: f_bound = array( [ [1.0 + (0.25 * a_coeff * a_coeff)], [1.0], [1 + 0.25 * b_coeff * b_coeff], ], dtype=self.dtype, ) 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 f_bound = array( [[1.0 - (0.5 * a_coeff)], [1.0], [1 + 0.5 * b_coeff]], dtype=self.dtype ) return s_shifted, f_bound @staticmethod def _normalize_s(s_ref, s_new): """Normalization of input variables for use of polynomial approximation. :param s_ref: vector of initial values of independent variables (5 variables at max) :type s_ref: numpy array :param s_new: vector of current values of independent variables :type s_new: numpy array :returns: normalized value and normalized+centered value :rtype: numpy array """ s_norm = s_new / s_ref s_norm[s_norm > 1.25] = 1.25 s_norm[s_norm < 0.75] = 0.75 s_shifted = s_norm - 1.0 # Shift S near origin return s_norm, s_shifted
[docs] def derive_normalize_s(self, s_ref, s_new): """Derivation of normalization of input variables for use of polynomial approximation. :param s_ref: vector of initial values of independent variables (5 variables at max) :type s_ref: numpy array :param s_new: vector of current values of independent variables :type s_new: numpy array :returns: derivatives of normalized value and normalized+centered value :rtype: numpy array """ if self.dtype == complex128: s_norm = s_new.real / s_ref.real else: s_norm = s_new / s_ref if s_norm > 1.25: derive_s_norm = 0.0 elif s_norm < 0.75: derive_s_norm = 0.0 else: derive_s_norm = 1.0 / s_ref # derive_s_norm = s_new / s_ref # derive_s_norm[derive_s_norm > 1.25] = 0.0 # derive_s_norm[derive_s_norm < 0.75] = 0.0 # derive_s_norm = derive_s_norm / s_ref return derive_s_norm
[docs] def derive_poly_approx(self, s_ref, s_new, flag, s_bound): """Compute the polynomial coefficients to characterize the behavior of certain synthetic variables and function modifiers. Compared to poly_approx, also returns polynomial coeff for linearization. :param s_ref: vector of initial values of independent variables (5 variables at max) :type s_ref: numpy array :param s_new: vector of current values of independent variables :type s_new: numpy array :param flag: indicates functional relationship between variables: - flag = 1: linear >0 - flag = 2: nonlinear >0 - flag = 3: linear < 0 - flag = 4: nonlinear <0 - flag = 5: parabolic :type flag: int :param s_bound: offset value for normalization :type s_bound: numpy array :returns: poly_value: value of synthetic variable or modifier :rtype: numpy array """ imax = s_ref.shape[0] # Normalize new S with initial _, s_shifted = self._normalize_s(s_ref, s_new) # derive_S_norm = self._derive_normalize_S(s_ref, s_new) ao_coeff = 0.0 ai_coeff = zeros(imax, dtype=self.dtype) aij_coeff = zeros((imax, imax), dtype=self.dtype) for i in range(imax): a_coeff = 0.1 b_coeff = a_coeff mtx_shifted = self.__compute_mtx_shifted(s_bound, i) s_shifted, f_bound = self.__compute_fbound( flag, s_shifted, a_coeff, b_coeff, i ) ao_coeff, ai_coeff, aij_coeff = self.__compute_a( mtx_shifted, f_bound, ao_coeff, ai_coeff, aij_coeff, i ) aij_coeff = self.__update_aij(aij_coeff, imax) poly_value = ( ao_coeff + dot(ai_coeff, s_shifted) + 0.5 * dot(dot(s_shifted, aij_coeff[:imax, :imax]), s_shifted) ) return poly_value[0], ai_coeff, aij_coeff[:imax, :imax], s_shifted
[docs] def poly_approx(self, s_ref, s_new, flag, s_bound): """This function calculates polynomial coefficients to characterize the behavior of certain synthetic variables and function modifiers. :param s_ref: vector of initial values of independent variables (5 variables at max) :type s_ref: numpy array :param s_new: vector of current values of independent variables :type s_new: numpy array :param flag: indicates functional relationship between variables: - flag = 1: linear >0 - flag = 2: nonlinear >0 - flag = 3: linear < 0 - flag = 4: nonlinear <0 - flag = 5: parabolic :type flag: int :param s_bound: offset value for normalization :type s_bound: numpy array :returns: poly_value: value of synthetic variable or modifier :rtype: numpy array """ imax = s_ref.shape[0] # Normalize new S with initial _, s_shifted = self._normalize_s(s_ref, s_new) # - ai_coeff: vector of coeff for 2nd term # - aij_coeff: matrix of coeff for 3rd term # - ao_coeff: scalar coeff ao_coeff = 0.0 ai_coeff = zeros(imax, dtype=self.dtype) aij_coeff = zeros((imax, imax), dtype=self.dtype) for i in range(imax): a_coeff = 0.1 b_coeff = a_coeff mtx_shifted = self.__compute_mtx_shifted(s_bound, i) s_shifted, f_bound = self.__compute_fbound( flag, s_shifted, a_coeff, b_coeff, i ) ao_coeff, ai_coeff, aij_coeff = self.__compute_a( mtx_shifted, f_bound, ao_coeff, ai_coeff, aij_coeff, i ) aij_coeff = self.__update_aij(aij_coeff, imax) poly_value = ( ao_coeff + dot(ai_coeff, s_shifted) + 0.5 * dot(dot(s_shifted, aij_coeff[:imax, :imax]), s_shifted) ) return poly_value[0]
[docs] def compute_half_span(self, x_shared): """Compute half-span from surface and aspect ratio. :param x_shared: global design variables :type x_shared: numpy array :returns: half-span :rtype: numpy array """ return self.math.sqrt(x_shared[3] * x_shared[5]) * 0.5 # 1/2 span
[docs] def compute_thickness(self, x_shared): """Compute a wing thickness. :param x_shared: global design variables :type x_shared: numpy array :returns: thickness :rtype: numpy array """ return ( x_shared[0] * x_shared[5] / (self.math.sqrt(x_shared[5] * x_shared[3])) ) # thickness
[docs] @staticmethod def compute_aero_center(x_1): """Computes the aerodynamic center. :param x_1: local design variables (structure) :type x_1: numpy array :returns: aerodynamic center :rtype: numpy array """ return (1.0 + 2.0 * x_1[0]) / (3.0 * (1 + x_1[0]))
[docs] def get_initial_values(self): """Get initial values used by polynomial functions. :returns: initial values :rtype: tuple(ndarray) """ self.i_0 = self.get_default_x0() x_initial = self.i_0[1] tc_initial = self.i_0[4] half_span_initial = self.math.sqrt(self.i_0[7] * self.i_0[9]) * 0.5 aero_center_initial = (1 + 2 * self.i_0[0]) / (3 * (1 + self.i_0[0])) cf_initial = self.i_0[2] mach_initial = self.i_0[6] h_initial = self.i_0[5] throttle_initial = self.i_0[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, )