# 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: Matthias De Lozzo
# OTHER AUTHORS - MACROSCOPIC CHANGES
r"""The aerostructure MDO problem.
The **aerostructure** module implements all :class:`.MDODiscipline`
included in the Aerostructure problem:
.. math::
\text{OVERALL AIRCRAFT DESIGN} = \left\{
\begin{aligned}
&\text{minimize }\text{range}(\text{thick\\_airfoils},
\text{thick\\_panels}, \text{sweep}) = 8
\times10^{11}\times\text{lift}\times\text{mass}/\text{drag} \\
&\text{with respect to }\text{thick\\_airfoils},\,\text{thick\\_panels},
\,\text{sweep} \\
&\text{subject to }\\
& \text{rf}-0.5 = 0\\
& \text{lift}-0.5 \leq 0
\end{aligned}\right.
where
.. math::
\text{AERODYNAMICS} = \left\{
\begin{aligned}
&\text{drag}=0.1\times((\text{sweep}/360)^2 + 200 +
\text{thick\\_airfoils}^2 - \text{thick\\_airfoils} -
4\times\text{displ})\\
&\text{forces}=10\times\text{sweep} +
0.2\times\text{thick\\_airfoils}-0.2\times\text{displ}\\
&\text{lift}=(\text{sweep} + 0.2\times\text{thick\\_airfoils}-
2\times\text{displ})/3000
\end{aligned}
\right.
and
.. math::
\text{STRUCTURE} = \left\{
\begin{aligned}
&\text{mass}=4000\times(\text{sweep}/360)^3 + 200000 +
100\times\text{thick\\_panels} + 200\times\text{forces}\\
&\text{rf}=3\times\text{sweep} - 6\times\text{thick\\_panels} +
0.1\times\text{forces} + 55\\
&\text{displ}=2\times\text{sweep} + 3\times\text{thick\\_panels} -
2\times\text{forces}
\end{aligned}
\right.
"""
from __future__ import annotations
from typing import TYPE_CHECKING
from numpy import array
from numpy import atleast_2d
from numpy import complex128
from numpy import ones
from gemseo.core.discipline import MDODiscipline
if TYPE_CHECKING:
from collections.abc import Iterable
[docs]
class Mission(MDODiscipline):
"""The mission discipline of the aerostructure use case.
Compute the objective and the constraints.
"""
def __init__(self, r_val: float = 0.5, lift_val: float = 0.5) -> None:
"""
Args:
r_val: The threshold to compute the reserve factor constraint.
lift_val: The threshold to compute the lift constraint.
""" # noqa: D205 D212
super().__init__(auto_detect_grammar_files=True)
self.default_inputs = get_inputs("lift", "mass", "drag", "reserve_fact")
self.re_exec_policy = self.ReExecutionPolicy.DONE
self.r_val = r_val
self.lift_val = lift_val
def _run(self) -> None:
lift, mass, drag, reserve_fact = self.get_inputs_by_name([
"lift",
"mass",
"drag",
"reserve_fact",
])
obj = array([self.compute_range(lift, mass, drag)], dtype=complex128)
c_lift = array([self.c_lift(lift, self.lift_val)], dtype=complex128)
c_rf = array([self.c_rf(reserve_fact)], dtype=complex128)
self.store_local_data(range=obj, c_lift=c_lift, c_rf=c_rf)
[docs]
@staticmethod
def compute_range(lift, mass, drag) -> float:
"""Compute the objective function: :math:`range=8.10^{11}*lift/(mass*drag)`.
Args:
lift: The lift.
mass: The mass.
drag: The drag.
Returns:
The range.
"""
return 8e11 * lift[0] / (mass[0] * drag[0])
[docs]
@staticmethod
def c_lift(lift, lift_val: float = 0.5):
"""Compute the lift constraint: :math:`lift-0.5`.
Args:
lift: The lift.
lift_val: The threshold for the lift constraint.
Returns:
The value of the lift constraint.
"""
return lift[0] - lift_val
[docs]
@staticmethod
def c_rf(reserve_fact, rf_val: float = 0.5):
"""Compute the reserve factor constraint: :math:`rf-0.5`.
Args:
reserve_fact: The reserve factor.
rf_val: The threshold for the reserve factor constraint.
Returns:
The value of the reserve factor constraint.
"""
return reserve_fact[0] - rf_val
def _compute_jacobian(
self, inputs: Iterable[str] | None = None, outputs: Iterable[str] | None = None
) -> None:
# Initialize all matrices to zeros
self._init_jacobian(inputs, outputs)
drag, lift, mass = self.get_inputs_by_name(["drag", "lift", "mass"])
self.jac["c_lift"]["lift"] = ones((1, 1))
self.jac["c_rf"]["reserve_fact"] = ones((1, 1))
self.jac["range"]["lift"] = atleast_2d(array([8e11 / (mass[0] * drag[0])]))
self.jac["range"]["mass"] = atleast_2d(
array([-8e11 * lift[0] / (mass[0] ** 2 * drag[0])])
)
self.jac["range"]["drag"] = atleast_2d(
array([-8e11 * lift[0] / (mass[0] * drag[0] ** 2)])
)
[docs]
class Aerodynamics(MDODiscipline):
"""The aerodynamics discipline of the aerostructure use case.
Evaluate: ``[drag, forces, lift] = f(sweep, thick_airfoils, displ)``.
"""
def __init__(self) -> None: # noqa: D107
super().__init__(auto_detect_grammar_files=True)
self.default_inputs = get_inputs("sweep", "thick_airfoils", "displ")
self.re_exec_policy = self.ReExecutionPolicy.DONE
def _run(self) -> None:
sweep, thick_airfoils, displ = self.get_inputs_by_name([
"sweep",
"thick_airfoils",
"displ",
])
drag_out = array([self.compute_drag(sweep, thick_airfoils, displ)])
lift_out = array([self.compute_lift(sweep, thick_airfoils, displ)])
forces_out = array([self.compute_forces(sweep, thick_airfoils, displ)])
self.store_local_data(drag=drag_out, forces=forces_out, lift=lift_out)
[docs]
@staticmethod
def compute_drag(sweep, thick_airfoils, displ) -> float:
r"""Compute the coupling.
:math:`drag=0.1*((sweep/360)^2 + 200 + thick\\_airfoils^2
- thick\\_airfoils - 4*displ)`
Args:
sweep: The sweep.
thick_airfoils: The thickness of the airfoils.
displ: The displacement.
Returns:
The drag.
"""
return 0.1 * (
(sweep[0] / 360) ** 2
+ 200
+ thick_airfoils[0] ** 2
- thick_airfoils[0]
- 4 * displ[0]
)
[docs]
@staticmethod
def compute_forces(sweep, thick_airfoils, displ):
r"""Compute the coupling forces.
:math:`forces=10*sweep + 0.2*thick\\_airfoils-0.2*displ`
Args:
sweep: The sweep.
thick_airfoils: The thickness of the airfoils.
displ: The displacement.
Returns:
The forces.
"""
return 10 * sweep[0] + 0.2 * thick_airfoils[0] - 0.2 * displ[0]
[docs]
@staticmethod
def compute_lift(sweep, thick_airfoils, displ):
r"""Compute the lift.
:math:`lift=(sweep + 0.2*thick\\_airfoils-2.*displ)/3000.`.
Args:
sweep: The sweep.
thick_airfoils: The thickness of the airfoils.
displ: The displacement.
Returns:
The lift.
"""
return (sweep[0] + 0.2 * thick_airfoils[0] - 2.0 * displ[0]) / 3000.0
def _compute_jacobian(
self, inputs: Iterable[str] | None = None, outputs: Iterable[str] | None = None
) -> None:
# Initialize all matrices to zeros
self._init_jacobian(inputs, outputs)
sweep, thick_airfoils = self.get_inputs_by_name(["sweep", "thick_airfoils"])
self.jac["drag"]["sweep"] = atleast_2d(
array([0.1 * 2.0 * sweep[0] / 360.0**2.0])
)
self.jac["drag"]["thick_airfoils"] = atleast_2d(
array([0.1 * (2.0 * thick_airfoils[0] - 1.0)])
)
self.jac["drag"]["displ"] = atleast_2d(0.1 * array([-4.0]))
self.jac["forces"]["sweep"] = atleast_2d(array([10.0]))
self.jac["forces"]["thick_airfoils"] = atleast_2d(array([0.2]))
self.jac["forces"]["displ"] = atleast_2d(array([-0.2]))
self.jac["lift"]["sweep"] = atleast_2d(array([1.0 / 3000.0]))
self.jac["lift"]["thick_airfoils"] = atleast_2d(array([0.2 / 3000.0]))
self.jac["lift"]["displ"] = atleast_2d(array([-2.0 / 3000.0]))
[docs]
class Structure(MDODiscipline):
"""The structure discipline of the aerostructure use case.
Evaluate: ``[mass, rf, displ] = f(sweep, thick_panels, forces)``.
"""
def __init__(self) -> None: # noqa: D107
super().__init__(auto_detect_grammar_files=True)
self.default_inputs = get_inputs("sweep", "forces", "thick_panels")
self.re_exec_policy = self.ReExecutionPolicy.DONE
def _run(self) -> None:
sweep, thick_panels, forces = self.get_inputs_by_name([
"sweep",
"thick_panels",
"forces",
])
mass_out = array([self.compute_mass(sweep, thick_panels, forces)])
rf_out = array([self.compute_rf(sweep, thick_panels, forces)])
displ_out = array([self.compute_displ(sweep, thick_panels, forces)])
self.store_local_data(mass=mass_out, reserve_fact=rf_out, displ=displ_out)
[docs]
@staticmethod
def compute_mass(sweep, thick_panels, forces):
r"""Compute the mass.
:math:`mass=4000*(sweep/360)^3 + 200000 + 100*thick\\_panels + 200.0*forces`.
Args:
sweep: The sweep.
thick_panels: The thickness of the panels.
forces: The forces.
Returns:
The mass.
"""
return (
4000 * (sweep[0] / 360) ** 3
+ 200000
+ 100 * thick_panels[0]
+ 200.0 * forces[0]
)
[docs]
@staticmethod
def compute_rf(sweep, thick_panels, forces):
r"""Compute the coupling.
:math:`rf=-3*sweep - 6*thick\\_panels + 0.1*forces + 55`
Args:
sweep: The sweep.
thick_panels: The thickness of the panels.
forces: The forces.
Returns:
The reserve factor.
"""
return -3 * sweep[0] - 6 * thick_panels[0] + 0.1 * forces[0] + 55
[docs]
@staticmethod
def compute_displ(sweep, thick_panels, forces):
r"""Compute the coupling.
:math:`displ=2*sweep + 3*thick\\_panels - 2.*forces`
Args:
sweep: The sweep.
thick_panels: The thickness of the panels.
forces: The forces.
Returns:
The displacement.
"""
return 2 * sweep[0] + 3 * thick_panels[0] - 2.0 * forces[0]
def _compute_jacobian(
self, inputs: Iterable[str] | None = None, outputs: Iterable[str] | None = None
) -> None:
# Initialize all matrices to zeros
self._init_jacobian(inputs, outputs)
sweep = self.get_inputs_by_name("sweep")
jac = self.jac["mass"]
jac["sweep"] = atleast_2d(array([4000.0 * 3.0 * sweep[0] ** 2 / 360.0**3]))
jac["thick_panels"] = atleast_2d(array([100.0]))
jac["forces"] = atleast_2d(array([200.0]))
jac = self.jac["reserve_fact"]
jac["sweep"] = atleast_2d(array([-3.0]))
jac["thick_panels"] = atleast_2d(array([-6.0]))
jac["forces"] = atleast_2d(array([0.1]))
jac = self.jac["displ"]
jac["sweep"] = atleast_2d(array([2.0]))
jac["thick_panels"] = atleast_2d([3.0])
jac["forces"] = atleast_2d(array([-2.0]))
