# 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
"""The Sobieski's SSBJ problem."""
from __future__ import annotations
import cmath
import logging
import math
import random
from collections import namedtuple
from pathlib import Path
from random import uniform
from typing import Iterable
from typing import Sequence
from numpy import array
from numpy import complex128
from numpy import concatenate
from numpy import float64
from numpy import ndarray
from numpy import ones
from gemseo.algos.design_space import DesignSpace
from gemseo.problems.sobieski.core.aerodynamics import SobieskiAerodynamics
from gemseo.problems.sobieski.core.mission import SobieskiMission
from gemseo.problems.sobieski.core.propulsion import SobieskiPropulsion
from gemseo.problems.sobieski.core.structure import SobieskiStructure
from gemseo.problems.sobieski.core.utils import SobieskiBase
LOGGER = logging.getLogger(__name__)
_Disciplines = namedtuple(
"Disciplines",
("aerodynamics", "structure", "propulsion", "mission"),
)
[docs]class SobieskiProblem:
r"""The Sobieski's SSBJ problem.
This problem seeks to maximize the range of a super-sonic business jet (SSBJ)
over a bounded design space whilst satisfying inequality constraints.
The objective and constraint functions come from a system of four disciplines:
1. the structure,
computing the total aircraft mass :math:`y_{1,4,0}` and :math:`y_{1,2,0}`,
the fuel mass :math:`y_{1,4,1}`,
the wing twist :math:`y_{1,2,1}`
and the five stress constraints :math:`g_{1,0},\ldots,g_{4,0}`.
2. the aerodynamics discipline,
computing the lift :math:`y_{2,1}`,
drag :math:`y_{2,3}`,
lift-to-drag ratio
and pressure gradient constraint :math:`g_2`,
3. the propulsion,
computing the engine mass :math:`y_{3,1,0}`,
the engine scale factor :math:`y_{3,2,0}`,
which is also the constraint :math:`g_{3,0}`,
the specific fuel consumption :math:`y_{3,4,0}`,
the engine temperature constraint :math:`g_{3,1}`,
and the throttle setting constraint :math:`g_{3,2}`
4. the mission,
computing the range :math:`y_{4,0}`.
Notes:
- The structure, aerodynamics, propulsion and mission disciplines
are numbered from 1 to 4.
- The variable :math:`y_{i,j,k}` is a coupling variable
from the discipline :math:`i` to the discipline :math:`j`.
- The aerodynamics, structure and propulsion disciplines are strongly coupled,
i.e. each of them depends directly or indirectly on the others,
and provide inputs to the mission discipline.
The design variables can be classified into four groups:
- the design variables which are inputs to at least two disciplines,
- :math:`x_{0,0}`: the thickness-to-chord ratio,
- :math:`x_{0,1}`: the altitude (ft),
- :math:`x_{0,2}`: the Mach number,
- :math:`x_{0,3}`: the aspect ratio,
- :math:`x_{0,4}`: the wing sweep (deg),
- :math:`x_{0,5}`: the wing surface area (ft :sup:`2`\),
- the design variables which are inputs of the structure discipline only:
- :math:`x_{1,0}`: the wing taper ratio,
- :math:`x_{1,1}`: the wingbox x-sectional area
- the design variables which are inputs of the aerodynamics discipline only:
- :math:`x_{2,0}`: the skin friction coefficient,
- the design variables which are inputs of the propulsion discipline only:
- :math:`x_{3,0}`: the throttle setting (engin mass flow).
Lastly,
this problem is based on five constants:
- :math:`c_0`: the minimum fuel weight,
- :math:`c_1`: the miscellaneous weight,
- :math:`c_2`: the maximum load factor,
- :math:`c_3`: the reference engine weight,
- :math:`c_4`: the minimum drag coefficient.
"""
CONTRAINTS_NAMES_INEQUALITY = (
"c_Stress_x1",
"c_Stress_x2",
"c_Stress_x3",
"c_Stress_x4",
"c_Stress_x5",
"c_Twist_upper",
"c_Twist_lower",
"c_Pgrad",
"c_ESF_upper",
"c_ESF_lower",
"c_Throttle",
"c_Temperature",
)
CONTRAINTS_NAMES = (
"Stress_x1",
"Stress_x2",
"Stress_x3",
"Stress_x4",
"Stress_x5",
"Twist",
"Pgrad",
"ESF",
"Temperature",
"Throttle",
)
DV_NAMES_NORMALIZED = (
"x_TaperRatio",
"x_SectionalArea",
"x_Cf",
"x_Throttle_setting",
"x_eta",
"x_h",
"x_Mach",
"x_AR",
"x_Phi",
"x_sref",
)
DV_NAMES = (
"TaperRatio",
"SectionalArea",
"Cf",
"Throttle_setting",
"eta",
"h",
"Mach",
"AR",
"Phi",
"sref",
)
COUPLING_VARIABLES_NAMES = (
"Total weight",
"Fuel weight",
"Wing twist",
"Lift",
"Drag",
"Lift/Drag",
"SFC",
"Engine weight",
)
STRESS_LIMIT = SobieskiStructure.STRESS_LIMIT
TWIST_UPPER_LIMIT = SobieskiStructure.TWIST_UPPER_LIMIT
TWIST_LOWER_LIMIT = SobieskiStructure.TWIST_LOWER_LIMIT
PRESSURE_GRADIENT_LIMIT = SobieskiAerodynamics.PRESSURE_GRADIENT_LIMIT
ESF_UPPER_LIMIT = SobieskiPropulsion.ESF_UPPER_LIMIT
ESF_LOWER_LIMIT = SobieskiPropulsion.ESF_LOWER_LIMIT
TEMPERATURE_LIMIT = SobieskiPropulsion.TEMPERATURE_LIMIT
def __init__(
self,
dtype: str = SobieskiBase.DTYPE_DOUBLE,
) -> None:
"""
Args:
dtype: The data type for the NumPy arrays, either "float64" or "complex128".
"""
if dtype == SobieskiBase.DTYPE_COMPLEX:
self.__dtype = complex128
self.__math = cmath
elif dtype == SobieskiBase.DTYPE_DOUBLE:
self.__math = math
self.__dtype = float64
else:
raise ValueError(f"Unknown dtype: {dtype}.")
self.__base = SobieskiBase(dtype=self.__dtype)
self.disciplines = _Disciplines(
aerodynamics=SobieskiAerodynamics(self.__base),
mission=SobieskiMission(self.__base),
propulsion=SobieskiPropulsion(self.__base),
structure=SobieskiStructure(self.__base),
)
(
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.__l_b, self.__u_b = self.__base.design_bounds
self.__design_space = None
self.__names_to_feasible_values = {
"x_1": array([0.14951, 7.5e-01], dtype=self.__dtype),
"x_2": array([7.5e-01], dtype=self.__dtype),
"x_3": array([0.1675], dtype=self.__dtype),
"x_shared": array(
[6.0e-02, 5.4e04, 1.4e00, 4.4e00, 6.6e01, 1.2e03], dtype=self.__dtype
),
}
@property
def initial_design(self) -> ndarray:
"""The initial design :math:`x`."""
return self.__base.initial_design.copy()
@property
def constants(self) -> ndarray:
"""The constant vector."""
return self.__base.constants
@property
def aerodynamics(self) -> SobieskiAerodynamics:
"""The aerodynamics discipline."""
return self.disciplines.aerodynamics
@property
def mission(self) -> SobieskiMission:
"""The mission discipline."""
return self.disciplines.mission
@property
def propulsion(self) -> SobieskiPropulsion:
"""The propulsion discipline."""
return self.disciplines.propulsion
@property
def structure(self) -> SobieskiStructure:
"""The structure discipline."""
return self.disciplines.structure
[docs] def get_bounds_by_name(
self,
variables_names: Sequence[str],
) -> tuple[ndarray, ndarray]:
"""Return the lower and upper bounds of variables.
Args:
variables_names: The names of the variables.
Returns:
The lower and upper bounds of the variables;
the array components keep the order of the variables.
"""
return self.__base.get_bounds_by_name(variables_names)
def __set_indata(
self,
design_vector: ndarray,
names: str | Iterable[str] | None,
) -> dict[str, ndarray]:
"""Return the default values of some variables of the problem.
Args:
design_vector: The design vector to be used.
names: The names of the variables of interest.
If ``None``, use all the variables of the problem.
Returns:
The default values of some variables of the problem.
"""
dtype = self.__dtype
constants = self.__base.constants
names_to_default_values = {
"x_shared": design_vector[4:10],
"y_12": array([50606.9742, 0.95], dtype=dtype),
"y_14": array((50606.9741711, 7306.20262124), dtype=dtype),
"y_21": array([50606.9741711], dtype=dtype),
"y_23": array([12562.01206488], dtype=dtype),
"y_24": array([4.15006276], dtype=dtype),
"y_31": array([6354.32430691], dtype=dtype),
"y_32": array([0.50279625], dtype=dtype),
"y_34": array([1.10754577], dtype=dtype),
"y_1": ones(3, dtype=dtype),
"y_2": ones(3, dtype=dtype),
"y_3": ones(3, dtype=dtype),
"g_1": ones(6, dtype=dtype),
"g_2": ones(1, dtype=dtype),
"g_3": ones(3, dtype=dtype),
"x_1": design_vector[:2],
"x_2": array([design_vector[2]], dtype=dtype),
"x_3": array([design_vector[3]], dtype=dtype),
"c_0": array([constants[0]], dtype=dtype),
"c_1": array([constants[1]], dtype=dtype),
"c_2": array([constants[2]], dtype=dtype),
"c_3": array([constants[3]], dtype=dtype),
"c_4": array([constants[4]], dtype=dtype),
}
if names is None:
return names_to_default_values
if isinstance(names, str):
names = [names]
return {
names: names_to_default_values[names]
for names in names
if names in names_to_default_values
}
[docs] def get_x0_feasible(self, names: str | Iterable[str] | None = None) -> ndarray:
"""Return a feasible starting point related to some input variables.
Args:
names: The names of the variables.
If ``None``, use all the input variables.
Returns:
The feasible starting point.
"""
if isinstance(names, str):
names = [names]
if names is None:
names = ["x_1", "x_2", "x_3", "x_shared"]
return concatenate([self.__names_to_feasible_values[name] for name in names])
@property
def design_bounds(self) -> tuple[ndarray, ndarray]:
"""The lower and upper bounds of the design variables."""
return self.__base.design_bounds
@property
def optimum_design(self) -> ndarray:
"""The optimal design vector found by Sobieski with BLISS."""
return array(
(0.38757, 0.75, 0.75, 0.15624, 0.06, 60000.0, 1.4, 2.5, 70.0, 1500.0),
dtype=self.__dtype,
)
@property
def optimum_range(self) -> ndarray:
"""The optimal range found by Sobieski with BLISS."""
return array([3963.98], dtype=self.__dtype)
[docs] def get_sobieski_constraints(
self,
g_1: ndarray,
g_2: ndarray,
g_3: ndarray,
true_cstr: ndarray = False,
) -> ndarray:
"""Return either the value of the constraints or the distance to the thresholds.
Args:
g_1: The constraints from the structure discipline:
``g_1[0]`` to ``g_1[4]`` are the stresses on wing
and ``g_1[5]`` is the wing twist.
g_2: The constraint (pressure gradient) from the aerodynamics discipline.
g_3: The constraints from the propulsion discipline:
``g_3[0]`` is the engine scale factor,
``g_3[1]`` is the engine temperature
and ``g_3[2]`` is the throttle setting constraint.
true_cstr: If ``True``,
return the value of the outputs.
Otherwise,
return the distance to the corresponding constraint thresholds.
Returns:
The constraints or the distance to the thresholds,
according to ``true_cstr``.
"""
if true_cstr:
return concatenate(
(g_1[0:5], array((g_1[5], g_2[0], g_3[0], g_3[2], g_3[1])))
)
return concatenate(
(
g_1[0:5] - self.STRESS_LIMIT,
array(
(
g_1[5] - self.TWIST_UPPER_LIMIT,
self.TWIST_LOWER_LIMIT - g_1[5],
g_2[0] - self.PRESSURE_GRADIENT_LIMIT,
g_3[0] - self.ESF_UPPER_LIMIT,
self.ESF_LOWER_LIMIT - g_3[0],
g_3[2],
g_3[1] - self.TEMPERATURE_LIMIT,
)
),
)
)
def __compute_mda(
self,
design_vector: tuple[ndarray],
true_cstr: bool = False,
accuracy: float = 1e-3,
) -> tuple[ndarray]:
"""Compute the output variables at equilibrium with the Gauss-Seidel algorithm.
Args:
design_vector: The design vector.
true_cstr: If ``True``,
return the value of the constraint outputs.
Otherwise,
return the distance to the corresponding constraint thresholds.
accuracy: The system resolution accuracy.
Returns:
The output variables at equilibrium:
- ``y_1``: The weight analysis outputs:
- ``y_1[0]``: The total aircraft weight,
- ``y_1[1]``: The fuel weight,
- ``y_1[2]``: The wing twist,
- ``y_2``: The outputs of the aerodynamics analysis:
- ``y_2[0]``: the lift,
- ``y_2[1]``: the drag,
- ``y_2[2]``: the lift/drag ratio,
- ``y_3``: The outputs of the propulsion analysis, namely
- ``y_3[0]``: the specific fuel consumption (SFC),
- ``y_3[1]``: the engine weight (EW),
- ``y_3[2]``: the engine scale factor (ESF),
- ``y_4``: The A/C range,
- ``y_12``: The coupling variables for the aerodynamics discipline:
- ``y_12[0]``: The total aircraft weight,
- ``y_12[1]``: The wing twist,
- ``y_14``: The coupling variables for the mission discipline:
- ``y_14[0]``: The total aircraft weight
- ``y_14[1]``: The fuel weight
- ``y_21``: The coupling variable (lift) for the structure discipline,
- ``y_23``: The coupling variable (drag) for the propulsion discipline,
- ``y_24``: The coupling variable (lift/drag ratio)
for the mission discipline,
- ``y_34``: The couping variable SFC for the mission discipline,
- ``y_31``: The coupling variable (EW) for the structure discipline,
- ``y_32``: The coupling variable (ESG) for the aerodynamics discipline,
- ``g_1``: the outputs to be constrained,
where ``g_1[0]`` to ``g_1[4]`` are stresses on the wing
and ``g_1[5]`` is the wing twist,
- ``g_2``: The pressure gradient to be constrained.
- ``g_3``: The propulsion constraints,
where ``g_3[0]`` is the engine scale factor constraint,
`g_3[1]``is the engine temperature,
and ``g_3[2]`` is the throttle setting constraint.
"""
dtype = self.__dtype
x_1 = design_vector[:2]
x_2 = array([design_vector[2]], dtype=dtype)
x_3 = array([design_vector[3]], dtype=dtype)
x_shared = design_vector[4:]
y_12 = ones(2, dtype=dtype)
y_21 = ones(1, dtype=dtype)
y_31 = ones(1, dtype=dtype)
y_32 = ones(1, dtype=dtype)
# Execute Gauss Seidel iteration on system to find Y variables
lift_convergence = y_21[0] + 10.0
eng_weight_conv = y_31[0] + 10.0
esf_convergence = y_32[0] + 10.0
# loop_index = 1
while (abs(lift_convergence - y_21[0]) > (y_21[0] * accuracy)) or (
(abs(eng_weight_conv - y_31[0]) > (y_31[0] * accuracy))
or (abs(esf_convergence - y_32[0]) > (y_32[0] * accuracy))
):
lift_convergence = y_21[0]
eng_weight_conv = y_31[0]
esf_convergence = y_32[0]
# Call Black Boxes
y_1, _, y_12, y_14, g_1 = self.structure.execute(
x_shared, y_21, y_31, x_1, true_cstr=true_cstr
)
y_2, y_21, y_23, y_24, g_2 = self.aerodynamics.execute(
x_shared, y_12, y_32, x_2, true_cstr=true_cstr
)
y_3, y_34, y_31, y_32, g_3 = self.propulsion.execute(
x_shared, y_23, x_3, true_cstr=true_cstr
)
y_4 = self.mission.execute(x_shared, y_14, y_24, y_34)
return (
y_1,
y_2,
y_3,
y_4,
y_12,
y_14,
y_21,
y_23,
y_24,
y_31,
y_32,
y_34,
g_1,
g_2,
g_3,
)
[docs] def get_constraints(
self,
design_vector: ndarray,
true_cstr: bool = False,
) -> tuple[ndarray, ndarray, ndarray]:
"""Compute all the constraints.
Args:
design_vector: The design vector.
true_cstr: If ``True``,
return the value of the constraint outputs.
Otherwise,
return the distance to the corresponding constraint thresholds.
Returns:
The value of the constraints :math:`g_1`, :math:`g_2` and :math:`g_3`.
"""
outputs = self.__compute_mda(design_vector, true_cstr=true_cstr)
return outputs[-3], outputs[-2], outputs[-1]
@property
def design_space(self) -> DesignSpace:
"""The design space."""
if self.__design_space is None:
input_file = Path(__file__).parent / "sobieski_design_space.txt"
self.__design_space = DesignSpace.read_from_txt(input_file)
if self.__dtype == complex128:
current_x = self.__design_space.get_current_value(as_dict=True)
for variable_name, variable_value in current_x.items():
current_x[variable_name] = array(variable_value, dtype=complex128)
self.__design_space.set_current_value(current_x)
return self.__design_space
@property
def design_space_with_physical_naming(self) -> DesignSpace:
"""The design space with physical naming."""
if self.__design_space is None:
input_file = Path(__file__).parent / "sobieski_design_space_pn.txt"
self.__design_space = DesignSpace.read_from_txt(input_file)
if self.__dtype == complex128:
current_x = self.__design_space.get_current_value(as_dict=True)
for variable_name, current_value in current_x.items():
current_x[variable_name] = array(current_value, dtype=complex128)
self.__design_space.set_current_value(current_x)
return self.__design_space
