diff --git a/changelog/fragments/1538.added.rst b/changelog/fragments/1538.added.rst new file mode 100755 index 0000000000000000000000000000000000000000..61e057a40dc52f073c3eaf8fe2ec798fe9af45d9 --- /dev/null +++ b/changelog/fragments/1538.added.rst @@ -0,0 +1 @@ +"BiLevelBCD" is a formulation derived from "BiLevel" that uses the Block Coordinate Descent algorithm to solve the lower-optimization problem. diff --git a/src/gemseo/formulations/bilevel.py b/src/gemseo/formulations/bilevel.py index 370829aaca57021634b1ffc810367036a6b46f3a..cee704e5bbb3fa929335c482b2baaf9077237464 100644 --- a/src/gemseo/formulations/bilevel.py +++ b/src/gemseo/formulations/bilevel.py @@ -68,13 +68,33 @@ class BiLevel(BaseMDOFormulation): 1. a first MDA to compute the coupling variables, 2. several disciplinary optimizations on the local design variables in parallel, 3. a second MDA to update the coupling variables. + + + The residual norm of MDA1 and MDA2 can be captured into scenario + observables thanks to different namespaces :attr:`.BiLevel.MDA1_RESIDUAL_NAMESPACE` + and :attr:`.BiLevel.MDA2_RESIDUAL_NAMESPACE`. """ DEFAULT_SCENARIO_RESULT_CLASS_NAME: ClassVar[str] = BiLevelScenarioResult.__name__ + """The default name of the scenario results.""" + + SYSTEM_LEVEL: ClassVar[str] = "system" + """The name of the system level.""" + + SUBSCENARIOS_LEVEL: ClassVar[str] = "sub-scenarios" + """The name of the sub-scenarios level.""" - SYSTEM_LEVEL = "system" - SUBSCENARIOS_LEVEL = "sub-scenarios" LEVELS = (SYSTEM_LEVEL, SUBSCENARIOS_LEVEL) + """The collection of levels.""" + + CHAIN_NAME: ClassVar[str] = "bilevel_chain" + """The name of the internal chain.""" + + MDA1_RESIDUAL_NAMESPACE: ClassVar[str] = "MDA1" + """The name of the namespace for the MDA1 residuals.""" + + MDA2_RESIDUAL_NAMESPACE: ClassVar[str] = "MDA2" + """The name of the namespace for the MDA2 residuals.""" Settings: ClassVar[type[BiLevel_Settings]] = BiLevel_Settings @@ -98,8 +118,7 @@ class BiLevel(BaseMDOFormulation): settings_model=settings_model, **settings, ) - self._shared_dv = design_space.variable_names - self.scenario_adapters = [] + self._scenario_adapters = [] self.coupling_structure = CouplingStructure( get_sub_disciplines(self.disciplines) ) @@ -126,10 +145,9 @@ class BiLevel(BaseMDOFormulation): def _build_scenario_adapters( self, output_functions: bool = False, - use_non_shared_vars: bool = False, adapter_class: type[MDOScenarioAdapter] = MDOScenarioAdapter, **adapter_options, - ) -> list[MDOScenarioAdapter]: + ) -> None: """Build the MDOScenarioAdapter required for each sub scenario. This is used to build the self.chain. @@ -137,20 +155,14 @@ class BiLevel(BaseMDOFormulation): Args: output_functions: Whether to add the optimization functions in the adapter outputs. - use_non_shared_vars: Whether the non-shared design variables are inputs - of the scenarios adapters. adapter_class: The class of the adapters. **adapter_options: The options for the adapters' initialization. - - Returns: - The adapters for the sub-scenarios. """ - adapters = [] scenario_log_level = adapter_options.pop( "scenario_log_level", self._settings.sub_scenarios_log_level ) for scenario in self.get_sub_scenarios(): - adapter_inputs = self._compute_adapter_inputs(scenario, use_non_shared_vars) + adapter_inputs = self._compute_adapter_inputs(scenario) adapter_outputs = self._compute_adapter_outputs(scenario, output_functions) adapter = adapter_class( scenario, @@ -159,8 +171,7 @@ class BiLevel(BaseMDOFormulation): scenario_log_level=scenario_log_level, **adapter_options, ) - adapters.append(adapter) - return adapters + self._scenario_adapters.append(adapter) def _compute_adapter_outputs( self, @@ -197,35 +208,39 @@ class BiLevel(BaseMDOFormulation): def _compute_adapter_inputs( self, scenario: BaseScenario, - use_non_shared_vars: bool, ) -> list[str]: """Compute the scenario adapter inputs. Args: scenario: A sub-scenario. - use_non_shared_vars: Whether to add the non-shared variables - as inputs of the adapter. Returns: The input variables of the adapter. """ - shared_dv = set(self._shared_dv) + local_dv = [ + var + for scn in self.get_sub_scenarios() + for var in scn.design_space.variable_names + ] + shared_dv = set(self.optimization_problem.design_space.variable_names) couplings = self.coupling_structure.all_couplings mda1_outputs = self._get_mda1_outputs() top_disc = scenario.formulation.get_top_level_disciplines() top_inputs = [inpt for disc in top_disc for inpt in disc.io.input_grammar] + nonshared_var = set(scenario.design_space.variable_names) + # All couplings of the scenarios are taken from the MDA - adapter_inputs = list( + return list( # Add shared variables from system scenario driver - set(top_inputs) & (set(couplings) | shared_dv | set(mda1_outputs)) - ) - if use_non_shared_vars: - adapter_inputs = list( - set(adapter_inputs) - | set(top_inputs).intersection(scenario.formulation.design_space) + set(top_inputs) + & ( + set(couplings) + | shared_dv + | set(mda1_outputs) + | set(local_dv) - nonshared_var ) - return adapter_inputs + ) def _get_mda1_outputs(self) -> list[str]: """Return the MDA1 outputs. @@ -300,6 +315,7 @@ class BiLevel(BaseMDOFormulation): settings_model=self._settings.main_mda_settings, ) mda1.settings.warm_start = True + else: LOGGER.warning( "No strongly coupled disciplines detected, " @@ -312,6 +328,7 @@ class BiLevel(BaseMDOFormulation): settings_model=self._settings.main_mda_settings, ) mda2.settings.warm_start = False + return mda1, mda2 def _build_chain_dis_sub_opts( @@ -322,7 +339,7 @@ class BiLevel(BaseMDOFormulation): Returns: The first MDA (if exists) and the sub-scenarios. """ - return [] if self._mda1 is None else [self._mda1], self.scenario_adapters + return [] if self._mda1 is None else [self._mda1], self._scenario_adapters def _create_multidisciplinary_chain(self) -> MDOChain: """Create the multidisciplinary chain. @@ -333,7 +350,7 @@ class BiLevel(BaseMDOFormulation): The multidisciplinary chain. """ # Build the scenario adapters to be chained with MDAs - self.scenario_adapters = self._build_scenario_adapters( + self._build_scenario_adapters( reset_x0_before_opt=self._settings.reset_x0_before_opt, keep_opt_history=True, ) @@ -354,11 +371,11 @@ class BiLevel(BaseMDOFormulation): chain_dis += [self._mda2] if self._settings.reset_x0_before_opt: - return MDOChain(chain_dis, name="bilevel_chain") + return MDOChain(chain_dis, name=self.CHAIN_NAME) return MDOWarmStartedChain( chain_dis, - name="bilevel_chain", + name=self.CHAIN_NAME, variable_names_to_warm_start=self._get_variable_names_to_warm_start(), ) @@ -372,7 +389,7 @@ class BiLevel(BaseMDOFormulation): """ variable_names = [ name - for adapter in self.scenario_adapters + for adapter in self._scenario_adapters for name in adapter.io.output_grammar ] if self._mda1: diff --git a/src/gemseo/formulations/bilevel_bcd.py b/src/gemseo/formulations/bilevel_bcd.py new file mode 100644 index 0000000000000000000000000000000000000000..a91df88ed3e34844447e5f425c67925988eedfd9 --- /dev/null +++ b/src/gemseo/formulations/bilevel_bcd.py @@ -0,0 +1,87 @@ +# 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 - API and implementation and/or documentation +# :author: Francois Gallard +# OTHER AUTHORS - MACROSCOPIC CHANGES +"""A Bi-level formulation using the block coordinate descent algorithm.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING +from typing import ClassVar + +from gemseo.core.chains.warm_started_chain import MDOWarmStartedChain +from gemseo.formulations.bilevel import BiLevel +from gemseo.formulations.bilevel_bcd_settings import BiLevel_BCD_Settings +from gemseo.mda.factory import MDAFactory + +if TYPE_CHECKING: + from gemseo.core.chains.chain import MDOChain + + +class BiLevelBCD(BiLevel): + """Block Coordinate Descent bi-level formulation. + + This formulation draws an optimization architecture + that involves multiple optimization problems to be solved + in order to obtain the solution of the MDO problem. + + At each iteration on the global design variables, + the bi-level-BCD MDO formulation implementation performs: + + 1. a first MDA to compute the coupling variables, + 2. a loop on several disciplinary optimizations on the + local design variables using an :class:`.MDAGaussSeidel`, + 3. an MDA to compute precisely the system optimization criteria. + """ + + Settings: ClassVar[type[BiLevel_BCD_Settings]] = BiLevel_BCD_Settings + + _settings: BiLevel_BCD_Settings + + __mda_factory: ClassVar[MDAFactory] = MDAFactory() + """The MDA factory.""" + + def _create_multidisciplinary_chain(self) -> MDOChain: + """Build the chain on top of which all functions are built. + + This chain is: MDA -> MDAGaussSeidel(MDOScenarios) -> MDA. + + Returns: + The multidisciplinary chain. + """ + # Build the scenario adapters to be chained with MDAs + self._build_scenario_adapters( + reset_x0_before_opt=self._settings.reset_x0_before_opt, + keep_opt_history=True, + ) + chain_disc, sub_opts = self._build_chain_dis_sub_opts() + + bcd_mda = self.__mda_factory.create( + "MDAGaussSeidel", + sub_opts, + settings_model=self._settings.bcd_mda_settings, + ) + + chain_disc += [bcd_mda] + if self._mda2: + chain_disc += [self._mda2] + + return MDOWarmStartedChain( + chain_disc, + name=self.CHAIN_NAME, + variable_names_to_warm_start=self._get_variable_names_to_warm_start(), + ) diff --git a/src/gemseo/formulations/bilevel_bcd_settings.py b/src/gemseo/formulations/bilevel_bcd_settings.py new file mode 100644 index 0000000000000000000000000000000000000000..c2a852d050adda6b57f4817ce422d884fc52d40e --- /dev/null +++ b/src/gemseo/formulations/bilevel_bcd_settings.py @@ -0,0 +1,52 @@ +# 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. +"""Settings of the BiLevel BCD formulation .""" + +from __future__ import annotations + +from functools import partial +from typing import TYPE_CHECKING + +from pydantic import Field +from pydantic import model_validator + +from gemseo.formulations.bilevel_settings import BiLevel_Settings +from gemseo.mda.base_mda_settings import BaseMDASettings +from gemseo.mda.gauss_seidel_settings import MDAGaussSeidel_Settings +from gemseo.utils.pydantic import copy_field + +if TYPE_CHECKING: + from typing_extensions import Self + + +copy_field_opt = partial(copy_field, model=BaseMDASettings) + + +class BiLevel_BCD_Settings(BiLevel_Settings): # noqa: N801 + """Settings of the :class:`.BiLevel` formulation.""" + + _TARGET_CLASS_NAME = "BiLevel_BCD" + + bcd_mda_settings: MDAGaussSeidel_Settings = Field( + default=MDAGaussSeidel_Settings(warm_start=True), + description="The settings for the MDA used in the BCD method.", + ) + + @model_validator(mode="after") + def __force_bcd_mda_warm_start(self) -> Self: + """Validates the state of the BCD MDA warm start setting as True.""" + if not self.bcd_mda_settings.warm_start: + self.bcd_mda_settings.warm_start = True + return self diff --git a/src/gemseo/scenarios/scenario_results/bilevel_scenario_result.py b/src/gemseo/scenarios/scenario_results/bilevel_scenario_result.py index 5827ab0845975f809f9c4c2e70b2934fa173fc25..126b6817aadb9c25739e9269a55b26ad45c69ed2 100644 --- a/src/gemseo/scenarios/scenario_results/bilevel_scenario_result.py +++ b/src/gemseo/scenarios/scenario_results/bilevel_scenario_result.py @@ -46,7 +46,7 @@ class BiLevelScenarioResult(ScenarioResult): main_problem = formulation.optimization_problem x_shared_opt = main_problem.solution.x_opt i_opt = main_problem.database.get_iteration(x_shared_opt) - 1 - scenario_adapters = formulation.scenario_adapters + scenario_adapters = formulation._scenario_adapters self.__n_sub_problems = len(scenario_adapters) for index, scenario_adapter in enumerate(scenario_adapters): sub_problem = scenario_adapter.scenario.formulation.optimization_problem diff --git a/src/gemseo/settings/formulations.py b/src/gemseo/settings/formulations.py index 2833c09bbe2b56ebed5ddc2812450dd9aadd2e1f..9e1cf4244c588b1699d2ef1ce25db700eb8f16be 100644 --- a/src/gemseo/settings/formulations.py +++ b/src/gemseo/settings/formulations.py @@ -16,6 +16,7 @@ from __future__ import annotations +from gemseo.formulations.bilevel_bcd_settings import BiLevel_BCD_Settings # noqa:F401 from gemseo.formulations.bilevel_settings import BiLevel_Settings # noqa: F401 from gemseo.formulations.disciplinary_opt_settings import ( # noqa: F401 DisciplinaryOpt_Settings, diff --git a/src/gemseo/utils/testing/bilevel_test_helper.py b/src/gemseo/utils/testing/bilevel_test_helper.py index be4536d4a49814b6e966da6cfb01bcb76d644e24..8189ac15bafb8ca4ad543da8142e4462fbc43dca 100644 --- a/src/gemseo/utils/testing/bilevel_test_helper.py +++ b/src/gemseo/utils/testing/bilevel_test_helper.py @@ -19,6 +19,14 @@ from __future__ import annotations from typing import Any from typing import Callable +from tests.core.test_dependency_graph import create_disciplines_from_desc + +from gemseo import create_design_space +from gemseo import create_discipline +from gemseo import create_scenario +from gemseo.problems.mdo.aerostructure.aerostructure_design_space import ( + AerostructureDesignSpace, +) from gemseo.problems.mdo.sobieski.disciplines import SobieskiAerodynamics from gemseo.problems.mdo.sobieski.disciplines import SobieskiMission from gemseo.problems.mdo.sobieski.disciplines import SobieskiProblem @@ -67,7 +75,11 @@ def create_sobieski_bilevel_scenario( def create_sobieski_sub_scenarios() -> tuple[MDOScenario, MDOScenario, MDOScenario]: - """Return the sub-scenarios of Sobieski's SuperSonic Business Jet.""" + """Creates the sub-scenarios for the Sobieski's SSBJ problem. + + Returns: + The sub-scenarios of Sobieski's SuperSonic Business Jet. + """ design_space = SobieskiProblem().design_space propulsion = MDOScenario( [SobieskiPropulsion()], @@ -99,3 +111,214 @@ def create_sobieski_sub_scenarios() -> tuple[MDOScenario, MDOScenario, MDOScenar ) return structure, aerodynamics, propulsion + + +def create_sobieski_bilevel_bcd_scenario() -> Callable[[dict[str, Any]], MDOScenario]: + """Create a function to generate a Sobieski Scenario. + + Returns: + A function which generates a Sobieski scenario with specific options. + """ + + def func(**settings): + """Create a Sobieski BiLevel scenario. + + Args: + **settings: The settings of the system scenario. + + Returns: + A Sobieski BiLevel Scenario. + """ + propulsion = SobieskiPropulsion() + aerodynamics = SobieskiAerodynamics() + struct = SobieskiStructure() + mission = SobieskiMission() + sub_disciplines = [struct, propulsion, aerodynamics, mission] + + ds = SobieskiProblem().design_space + + def create_block(design_var, name="MDOScenario"): + scenario = MDOScenario( + sub_disciplines, + "y_4", + ds.filter([design_var], copy=True), + formulation_name="MDF", + main_mda_name="MDAGaussSeidel", + maximize_objective=True, + name=name, + ) + scenario.set_algorithm(max_iter=50, algo_name="SLSQP") + scenario.formulation.optimization_problem.objective *= 0.001 + return scenario + + sc_prop = create_block("x_3", "PropulsionScenario") + sc_prop.add_constraint("g_3", constraint_type="ineq") + + sc_aero = create_block("x_2", "AerodynamicsScenario") + sc_aero.add_constraint("g_2", constraint_type="ineq") + + sc_str = create_block("x_1", "StructureScenario") + sc_str.add_constraint("g_1", constraint_type="ineq") + + # Gather the sub-scenarios and mission for objective computation + sub_scenarios = [sc_aero, sc_str, sc_prop, sub_disciplines[-1]] + + sc_system = MDOScenario( + sub_scenarios, + "y_4", + ds.filter(["x_shared"], copy=True), + formulation_name="BiLevelBCD", + maximize_objective=True, + bcd_mda_settings={"tolerance": 1e-5, "max_mda_iter": 10}, + **settings, + ) + sc_system.formulation.optimization_problem.objective *= 0.001 + sc_system.set_differentiation_method("finite_differences") + return sc_system + + return func + + +def create_dummy_bilevel_scenario(formulation_name: str) -> MDOScenario: + """Create a dummy BiLevel scenario. + + It has to be noted that there is no strongly coupled discipline in this example. + It implies that MDA1 will not be created. Yet, MDA2 will be created, + as it is built with all the sub-disciplines passed to the BiLevel formulation. + + Args: + formulation_name: The name of the BiLevel formulation to be used. + + Returns: + A dummy BiLevel MDOScenario. + """ + disc_expressions = { + "disc_1": (["x_1"], ["a"]), + "disc_2": (["a", "x_2"], ["b"]), + "disc_3": (["x", "x_3", "b"], ["obj"]), + } + discipline_1, discipline_2, discipline_3 = create_disciplines_from_desc( + disc_expressions + ) + + system_design_space = create_design_space() + system_design_space.add_variable("x_3") + + sub_design_space_1 = create_design_space() + sub_design_space_1.add_variable("x_1") + sub_scenario_1 = create_scenario( + [discipline_1, discipline_3] + + ([discipline_2] if formulation_name == "BiLevelBCD" else []), + "obj", + sub_design_space_1, + formulation_name="MDF", + ) + + sub_design_space_2 = create_design_space() + sub_design_space_2.add_variable("x_2") + sub_scenario_2 = create_scenario( + [discipline_2, discipline_3] + + ([discipline_1] if formulation_name == "BiLevelBCD" else []), + "obj", + sub_design_space_2, + formulation_name="MDF", + ) + + return create_scenario( + [sub_scenario_1, sub_scenario_2], + "obj", + design_space=system_design_space, + formulation_name=formulation_name, + ) + + +def create_aerostructure_scenario(formulation_name: str): + """Create an Aerostructure scenario. + + Args: + formulation_name: The name of the BiLevel formulation to be used. + + Returns: + An executed Aerostructure scenario. + """ + algo_settings = { + "xtol_rel": 1e-8, + "xtol_abs": 1e-8, + "ftol_rel": 1e-8, + "ftol_abs": 1e-8, + "ineq_tolerance": 1e-5, + "eq_tolerance": 1e-3, + } + + aero_formulas = { + "drag": "0.1*((sweep/360)**2 + 200 + thick_airfoils**2 - thick_airfoils - " + "4*displ)", + "forces": "10*sweep + 0.2*thick_airfoils - 0.2*displ", + "lift": "(sweep + 0.2*thick_airfoils - 2.*displ)/3000.", + } + aerodynamics = create_discipline( + "AnalyticDiscipline", name="Aerodynamics", expressions=aero_formulas + ) + struc_formulas = { + "mass": "4000*(sweep/360)**3 + 200000 + 100*thick_panels + 200.0*forces", + "reserve_fact": "-3*sweep - 6*thick_panels + 0.1*forces + 55", + "displ": "2*sweep + 3*thick_panels - 2.*forces", + } + structure = create_discipline( + "AnalyticDiscipline", name="Structure", expressions=struc_formulas + ) + mission_formulas = {"range": "8e11*lift/(mass*drag)"} + mission = create_discipline( + "AnalyticDiscipline", name="Mission", expressions=mission_formulas + ) + sub_scenario_options = { + "max_iter": 2, + "algo_name": "NLOPT_SLSQP", + **algo_settings, + } + design_space_ref = AerostructureDesignSpace() + + design_space_aero = design_space_ref.filter(["thick_airfoils"], copy=True) + + aero_scenario = create_scenario( + [aerodynamics, mission] + + ([structure] if formulation_name == "BiLevelBCD" else []), + "range", + design_space=design_space_aero, + formulation_name=( + "DisciplinaryOpt" if formulation_name == "BiLevel" else "MDF" + ), + maximize_objective=True, + ) + + aero_scenario.set_algorithm(**sub_scenario_options) + + design_space_struct = design_space_ref.filter(["thick_panels"], copy=True) + struct_scenario = create_scenario( + [structure, mission] + + ([aerodynamics] if formulation_name == "BiLevelBCD" else []), + "range", + design_space=design_space_struct, + formulation_name=( + "DisciplinaryOpt" if formulation_name == "BiLevel" else "MDF" + ), + maximize_objective=True, + ) + struct_scenario.set_algorithm(**sub_scenario_options) + + design_space_system = design_space_ref.filter(["sweep"], copy=True) + system_scenario = create_scenario( + [aero_scenario, struct_scenario, mission], + "range", + design_space=design_space_system, + formulation_name=formulation_name, + maximize_objective=True, + main_mda_name="MDAJacobi", + main_mda_settings={"tolerance": 1e-8}, + ) + + system_scenario.add_constraint("reserve_fact", constraint_type="ineq", value=0.5) + system_scenario.add_constraint("lift", value=0.5) + system_scenario.execute(algo_name="NLOPT_COBYLA", max_iter=5, **algo_settings) + + return system_scenario diff --git a/tests/core/test_factory.py b/tests/core/test_factory.py index 97f93a917f878ad30cedeaefd3605491b5bee9c6..e589ab14e1909718e9d523f2d93b25e74fe31a24 100644 --- a/tests/core/test_factory.py +++ b/tests/core/test_factory.py @@ -79,7 +79,7 @@ def test_print_configuration(reset_factory) -> None: assert re.findall(pattern, line) # check table body - formulations = ["BiLevel", "DisciplinaryOpt", "IDF", "MDF"] + formulations = ["BiLevel", "BiLevelBCD", "DisciplinaryOpt", "IDF", "MDF"] for formulation in formulations: pattern = f"\\|\\s+{formulation}\\s+\\|\\s+Yes\\s+\\|.+\\|" diff --git a/tests/formulations/test_bilevel.py b/tests/formulations/test_bilevel.py index 899fe4c6d2a3c160ca2c9647f910cb5e06418d65..0c6454a7d0e50c4e163fd54abffb417cbce8e421 100644 --- a/tests/formulations/test_bilevel.py +++ b/tests/formulations/test_bilevel.py @@ -22,27 +22,33 @@ import logging import pytest -from gemseo import create_design_space from gemseo import create_discipline -from gemseo import create_scenario from gemseo.algos.design_space import DesignSpace +from gemseo.algos.optimization_result import OptimizationResult from gemseo.core.chains.warm_started_chain import MDOWarmStartedChain from gemseo.core.discipline import Discipline from gemseo.disciplines.analytic import AnalyticDiscipline from gemseo.disciplines.utils import get_sub_disciplines from gemseo.formulations.bilevel import BiLevel -from gemseo.problems.mdo.aerostructure.aerostructure_design_space import ( - AerostructureDesignSpace, -) +from gemseo.formulations.bilevel_bcd import BiLevelBCD from gemseo.problems.mdo.sobieski.core.problem import SobieskiProblem from gemseo.problems.mdo.sobieski.disciplines import SobieskiAerodynamics from gemseo.problems.mdo.sobieski.disciplines import SobieskiMission from gemseo.problems.mdo.sobieski.disciplines import SobieskiPropulsion from gemseo.problems.mdo.sobieski.disciplines import SobieskiStructure from gemseo.scenarios.mdo_scenario import MDOScenario +from gemseo.utils.testing.bilevel_test_helper import create_aerostructure_scenario +from gemseo.utils.testing.bilevel_test_helper import create_dummy_bilevel_scenario +from gemseo.utils.testing.bilevel_test_helper import ( + create_sobieski_bilevel_bcd_scenario, +) from gemseo.utils.testing.bilevel_test_helper import create_sobieski_bilevel_scenario from gemseo.utils.testing.bilevel_test_helper import create_sobieski_sub_scenarios -from tests.core.test_dependency_graph import create_disciplines_from_desc + + +@pytest.fixture(params=["BiLevel", "BiLevelBCD"]) +def scenario_formulation_name(request): + return request.param @pytest.fixture @@ -52,54 +58,29 @@ def sobieski_bilevel_scenario(): @pytest.fixture -def dummy_bilevel_scenario() -> MDOScenario: - """Create a dummy BiLevel scenario. +def sobieski_bilevel_bcd_scenario(): + """Fixture from an existing function.""" + return create_sobieski_bilevel_bcd_scenario() - It has to be noted that there is no strongly coupled discipline in this example. - It implies that MDA1 will not be created. Yet, MDA2 will be created, - as it is built with all the sub-disciplines passed to the BiLevel formulation. - Returns: A dummy BiLevel MDOScenario. - """ - disc_expressions = { - "disc_1": (["x_1"], ["a"]), - "disc_2": (["a", "x_2"], ["b"]), - "disc_3": (["x", "x_3", "b"], ["obj"]), - } - discipline_1, discipline_2, discipline_3 = create_disciplines_from_desc( - disc_expressions - ) +@pytest.fixture +def dummy_bilevel_scenario(scenario_formulation_name) -> MDOScenario: + """Fixture from an existing function.""" + return create_dummy_bilevel_scenario(scenario_formulation_name) - system_design_space = create_design_space() - system_design_space.add_variable("x_3") - sub_design_space_1 = create_design_space() - sub_design_space_1.add_variable("x_1") - sub_scenario_1 = create_scenario( - [discipline_1, discipline_3], - "obj", - sub_design_space_1, - formulation_name="MDF", - ) +@pytest.fixture +def aerostructure_scenario(scenario_formulation_name) -> MDOScenario: + return create_aerostructure_scenario(scenario_formulation_name) - sub_design_space_2 = create_design_space() - sub_design_space_2.add_variable("x_2") - sub_scenario_2 = create_scenario( - [discipline_2, discipline_3], - "obj", - sub_design_space_2, - formulation_name="MDF", - ) - return create_scenario( - [sub_scenario_1, sub_scenario_2], - "obj", - system_design_space, - formulation_name="BiLevel", - ) +@pytest.fixture +def sobieski_sub_scenarios() -> tuple[MDOScenario, MDOScenario, MDOScenario]: + """Fixture from an existing function.""" + return create_sobieski_sub_scenarios() -def test_execute(sobieski_bilevel_scenario) -> None: +def test_constraint_not_in_sub_scenario(sobieski_bilevel_scenario) -> None: """Test the execution of the Sobieski BiLevel Scenario.""" scenario = sobieski_bilevel_scenario( apply_cstr_tosub_scenarios=True, apply_cstr_to_system=False @@ -136,73 +117,12 @@ def test_get_sub_options_grammar() -> None: assert "acceleration_method" in sub_option_values -def test_bilevel_aerostructure() -> None: +def test_bilevel_aerostructure(aerostructure_scenario) -> None: """Test the Bi-level formulation on the aero-structure problem.""" - algo_options = { - "xtol_rel": 1e-8, - "xtol_abs": 1e-8, - "ftol_rel": 1e-8, - "ftol_abs": 1e-8, - "ineq_tolerance": 1e-5, - "eq_tolerance": 1e-3, - } + scenario = aerostructure_scenario - aero_formulas = { - "drag": "0.1*((sweep/360)**2 + 200 + thick_airfoils**2 - thick_airfoils - " - "4*displ)", - "forces": "10*sweep + 0.2*thick_airfoils - 0.2*displ", - "lift": "(sweep + 0.2*thick_airfoils - 2.*displ)/3000.", - } - aerodynamics = create_discipline( - "AnalyticDiscipline", name="Aerodynamics", expressions=aero_formulas - ) - struc_formulas = { - "mass": "4000*(sweep/360)**3 + 200000 + 100*thick_panels + 200.0*forces", - "reserve_fact": "-3*sweep - 6*thick_panels + 0.1*forces + 55", - "displ": "2*sweep + 3*thick_panels - 2.*forces", - } - structure = create_discipline( - "AnalyticDiscipline", name="Structure", expressions=struc_formulas - ) - mission_formulas = {"range": "8e11*lift/(mass*drag)"} - mission = create_discipline( - "AnalyticDiscipline", name="Mission", expressions=mission_formulas - ) - design_space_ref = AerostructureDesignSpace() - - design_space_aero = design_space_ref.filter(["thick_airfoils"], copy=True) - aero_scenario = create_scenario( - [aerodynamics, mission], - "range", - design_space_aero, - formulation_name="DisciplinaryOpt", - maximize_objective=True, - ) - aero_scenario.set_algorithm(algo_name="NLOPT_SLSQP", max_iter=2, **algo_options) - - design_space_struct = design_space_ref.filter(["thick_panels"], copy=True) - struct_scenario = create_scenario( - [structure, mission], - "range", - design_space_struct, - formulation_name="DisciplinaryOpt", - maximize_objective=True, - ) - struct_scenario.set_algorithm(algo_name="NLOPT_SLSQP", max_iter=2, **algo_options) - - design_space_system = design_space_ref.filter(["sweep"], copy=True) - system_scenario = create_scenario( - [aero_scenario, struct_scenario, mission], - "range", - design_space_system, - formulation_name="BiLevel", - maximize_objective=True, - main_mda_name="MDAJacobi", - main_mda_settings={"tolerance": 1e-8}, - ) - system_scenario.add_constraint("reserve_fact", constraint_type="ineq", value=0.5) - system_scenario.add_constraint("lift", value=0.5) - system_scenario.execute(algo_name="NLOPT_COBYLA", max_iter=5, **algo_options) + assert isinstance(scenario.optimization_result, OptimizationResult) + assert scenario.formulation.optimization_problem.database.n_iterations == 5 def test_bilevel_weak_couplings(dummy_bilevel_scenario) -> None: @@ -216,18 +136,23 @@ def test_bilevel_weak_couplings(dummy_bilevel_scenario) -> None: # a and b are weak couplings of all the disciplines, # and they are in the top-level outputs of the first adapter disciplines = dummy_bilevel_scenario.formulation.chain.disciplines - assert "b" in disciplines[0].io.input_grammar assert "a" in disciplines[0].io.output_grammar + assert "b" in disciplines[1].io.output_grammar - # a is a weak coupling of all the disciplines, - # and it is in the top-level inputs of the second adapter - assert "a" in disciplines[1].io.input_grammar + if dummy_bilevel_scenario.formulation == BiLevel: + assert "b" in disciplines[0].io.input_grammar + assert "a" in disciplines[1].io.input_grammar - # a is a weak coupling of all the disciplines, - # and is in the top-level inputs of the second adapter - assert "b" in disciplines[1].io.output_grammar + if dummy_bilevel_scenario.formulation == BiLevelBCD: + assert "b" in disciplines[0].io.output_grammar + assert "a" in disciplines[1].io.output_grammar + assert "a" not in disciplines[0].io.input_grammar + assert "b" not in disciplines[0].io.input_grammar + assert "a" not in disciplines[1].io.input_grammar + assert "b" not in disciplines[1].io.input_grammar +@pytest.mark.parametrize("scenario_formulation_name", ["BiLevel"], indirect=True) def test_bilevel_mda_getter(dummy_bilevel_scenario) -> None: """Test that the user can access the MDA1 and MDA2.""" # In the Dummy scenario, there's not strongly coupled disciplines -> No MDA1 @@ -235,6 +160,7 @@ def test_bilevel_mda_getter(dummy_bilevel_scenario) -> None: assert "obj" in dummy_bilevel_scenario.formulation.mda2.io.output_grammar +@pytest.mark.parametrize("scenario_formulation_name", ["BiLevel"], indirect=True) def test_bilevel_mda_setter(dummy_bilevel_scenario) -> None: """Test that the user cannot modify the MDA1 and MDA2 after instantiation.""" discipline = create_discipline("SellarSystem") @@ -244,13 +170,16 @@ def test_bilevel_mda_setter(dummy_bilevel_scenario) -> None: dummy_bilevel_scenario.formulation.mda2 = discipline -def test_get_sub_disciplines(sobieski_bilevel_scenario) -> None: +@pytest.mark.parametrize( + "scenario", [sobieski_bilevel_scenario, sobieski_bilevel_bcd_scenario] +) +def test_get_sub_disciplines(scenario, request) -> None: """Test the get_sub_disciplines method with the BiLevel formulation. Args: sobieski_bilevel_scenario: Fixture to instantiate a Sobieski BiLevel Scenario. """ - scenario = sobieski_bilevel_scenario() + scenario = request.getfixturevalue(scenario.__name__)() classes = [ discipline.__class__ for discipline in get_sub_disciplines(scenario.formulation.disciplines) @@ -264,13 +193,16 @@ def test_get_sub_disciplines(sobieski_bilevel_scenario) -> None: } -def test_bilevel_warm_start(sobieski_bilevel_scenario) -> None: +@pytest.mark.parametrize( + "scenario", [sobieski_bilevel_scenario, sobieski_bilevel_bcd_scenario] +) +def test_bilevel_warm_start(scenario, request) -> None: """Test the warm start of the BiLevel chain. Args: sobieski_bilevel_scenario: Fixture to instantiate a Sobieski BiLevel Scenario. """ - scenario = sobieski_bilevel_scenario() + scenario = request.getfixturevalue(scenario.__name__)() scenario.formulation.chain.set_cache(Discipline.CacheType.MEMORY_FULL) bilevel_chain_cache = scenario.formulation.chain.cache scenario.formulation.chain.disciplines[0].set_cache( @@ -285,19 +217,21 @@ def test_bilevel_warm_start(sobieski_bilevel_scenario) -> None: assert (mda1_inputs[1]["y_12"] == chain_outputs[0]["y_12"]).all() assert mda1_inputs[2]["y_21"] == chain_outputs[1]["y_21"] assert (mda1_inputs[2]["y_12"] == chain_outputs[1]["y_12"]).all() + assert mda1_inputs[1]["y_32"] == chain_outputs[0]["y_32"] + assert (mda1_inputs[1]["y_23"] == chain_outputs[0]["y_23"]).all() + assert mda1_inputs[2]["y_32"] == chain_outputs[1]["y_32"] + assert (mda1_inputs[2]["y_23"] == chain_outputs[1]["y_23"]).all() +@pytest.mark.parametrize("scenario_formulation_name", ["BiLevel"], indirect=True) def test_bilevel_warm_start_no_mda1(dummy_bilevel_scenario) -> None: """Test that a warm start chain is built even if the process does not include any MDA1. - - Args: - dummy_bilevel_scenario: Fixture to instantiate a dummy weakly - coupled scenario. """ assert isinstance(dummy_bilevel_scenario.formulation.chain, MDOWarmStartedChain) +@pytest.mark.parametrize("scenario_formulation_name", ["BiLevel"], indirect=True) def test_bilevel_get_variable_names_to_warm_start_without_mdas( dummy_bilevel_scenario, monkeypatch ) -> None: @@ -307,13 +241,12 @@ def test_bilevel_get_variable_names_to_warm_start_without_mdas( def _no_mda2(*args, **kwargs): return None - scenario = dummy_bilevel_scenario - monkeypatch.setattr(scenario.formulation, "_mda2", _no_mda2()) + monkeypatch.setattr(dummy_bilevel_scenario.formulation, "_mda2", _no_mda2()) variables = [] - for adapter in scenario.formulation.scenario_adapters: + for adapter in dummy_bilevel_scenario.formulation._scenario_adapters: variables.extend(adapter.io.output_grammar) assert sorted(set(variables)) == sorted( - scenario.formulation._get_variable_names_to_warm_start() + dummy_bilevel_scenario.formulation._get_variable_names_to_warm_start() ) @@ -330,15 +263,21 @@ def test_test_bilevel_get_variable_names_to_warm_start_from_mdas( @pytest.mark.parametrize( - "options", + ("settings", "sub_scenario_formulation", "scenario_formulation"), [ - {}, - {"sub_scenarios_log_level": None}, - {"sub_scenarios_log_level": logging.INFO}, - {"sub_scenarios_log_level": logging.WARNING}, + ({}, "MDF", "BiLevelBCD"), + ({"sub_scenarios_log_level": None}, "MDF", "BiLevelBCD"), + ({"sub_scenarios_log_level": logging.INFO}, "MDF", "BiLevelBCD"), + ({"sub_scenarios_log_level": logging.WARNING}, "MDF", "BiLevelBCD"), + ({}, "DisciplinaryOpt", "BiLevel"), + ({"sub_scenarios_log_level": None}, "DisciplinaryOpt", "BiLevel"), + ({"sub_scenarios_log_level": logging.INFO}, "DisciplinaryOpt", "BiLevel"), + ({"sub_scenarios_log_level": logging.WARNING}, "DisciplinaryOpt", "BiLevel"), ], ) -def test_scenario_log_level(caplog, options) -> None: +def test_scenario_log_level( + caplog, settings, sub_scenario_formulation, scenario_formulation +) -> None: """Check scenario_log_level.""" design_space = DesignSpace() design_space.add_variable("x", lower_bound=0.0, upper_bound=1.0, value=0.5) @@ -347,7 +286,7 @@ def test_scenario_log_level(caplog, options) -> None: [AnalyticDiscipline({"z": "(x+y)**2"})], "z", design_space.filter(["y"], copy=True), - formulation_name="DisciplinaryOpt", + formulation_name=sub_scenario_formulation, name="FooScenario", ) sub_scenario.set_algorithm(algo_name="NLOPT_COBYLA", max_iter=2) @@ -355,23 +294,17 @@ def test_scenario_log_level(caplog, options) -> None: [sub_scenario], "z", design_space.filter(["x"]), - formulation_name="BiLevel", - **options, + formulation_name=scenario_formulation, + **settings, ) scenario.execute(algo_name="NLOPT_COBYLA", max_iter=2) - sub_scenarios_log_level = options.get("sub_scenarios_log_level") + sub_scenarios_log_level = settings.get("sub_scenarios_log_level") if sub_scenarios_log_level == logging.WARNING: assert "Start FooScenario execution" not in caplog.text else: assert "Start FooScenario execution" in caplog.text -@pytest.fixture -def sobieski_sub_scenarios() -> tuple[MDOScenario, MDOScenario, MDOScenario]: - """Fixture from an existing function.""" - return create_sobieski_sub_scenarios() - - def test_remove_couplings_from_ds(sobieski_sub_scenarios) -> None: """Check the removal of strong couplings for the design space.""" formulation = BiLevel( @@ -381,3 +314,97 @@ def test_remove_couplings_from_ds(sobieski_sub_scenarios) -> None: ) for strong_coupling in ["y_12", "y_21", "y_23", "y_31", "y_32"]: assert strong_coupling not in formulation.design_space + + +@pytest.mark.parametrize( + ("scenario", "subscenario"), + [ + ( + sobieski_bilevel_bcd_scenario, + { + "StructureScenario_adapter": { + "ssbj_local_variables": {"x_1", "x_2", "x_3"}, + "ssbj_input_couplings": {"y_31", "y_21", "y_23"}, + }, + "PropulsionScenario_adapter": { + "ssbj_local_variables": {"x_1", "x_2", "x_3"}, + "ssbj_input_couplings": {"y_31", "y_21", "y_23"}, + }, + "AerodynamicsScenario_adapter": { + "ssbj_local_variables": {"x_1", "x_2", "x_3"}, + "ssbj_input_couplings": {"y_31", "y_21", "y_23"}, + }, + }, + ), + ( + sobieski_bilevel_scenario, + { + "StructureScenario_adapter": { + "ssbj_local_variables": {"x_1"}, + "ssbj_input_couplings": {"y_31", "y_21"}, + }, + "PropulsionScenario_adapter": { + "ssbj_local_variables": {"x_3"}, + "ssbj_input_couplings": {"y_23"}, + }, + "AerodynamicsScenario_adapter": { + "ssbj_local_variables": {"x_2"}, + "ssbj_input_couplings": {"y_32", "y_12"}, + }, + }, + ), + ], +) +def test_adapters_inputs_outputs(scenario, subscenario, request) -> None: + """Test that the ScenarioAdapters within the BCD loop have the right inputs and + outputs. + """ + scenario = request.getfixturevalue(scenario.__name__)() + all_ssbj_couplings = { + "y_14", + "y_12", + "y_32", + "y_31", + "y_34", + "y_21", + "y_23", + "y_24", + } + + # Necessary couplings as inputs, + # depends on the order of the disciplines within the block MDAs. + ssbj_shared_variables = {"x_shared"} + for scenario_adapter in scenario.formulation._scenario_adapters: + ssbj_local_variables = subscenario[scenario_adapter.name][ + "ssbj_local_variables" + ] + ssbj_input_couplings = subscenario[scenario_adapter.name][ + "ssbj_input_couplings" + ] + adapter = scenario_adapter + + design_variable = set( + adapter.scenario.formulation.optimization_problem.design_space.variable_names + ) + other_local = ssbj_local_variables.difference(design_variable) + # Check the inputs + inputs = set(adapter.io.input_grammar) + # Check the outputs + outputs = set(adapter.io.output_grammar) + # Shared variables should always be present. + assert ssbj_shared_variables.issubset(inputs) + # Only necessary couplings should always be present. + assert ssbj_input_couplings.issubset(inputs) + # All local variables, excepted the optimized ones, should be present. + assert other_local.issubset(inputs) + assert not design_variable.issubset(inputs) + + # Shared variables should never be present + assert not ssbj_shared_variables.issubset(outputs) + # Only the optimized local variables should be present. + assert design_variable.issubset(outputs) + + if isinstance(scenario.formulation, BiLevelBCD): + # All couplings should always be present + assert all_ssbj_couplings.issubset(outputs) + assert not other_local.issubset(outputs) diff --git a/tests/formulations/test_formulations_factory.py b/tests/formulations/test_formulations_factory.py index 9108d54a6d44556db8b7e836b2d9d8f000c21f10..1c8f66b57a32de8b4c85d373e0823382b82892c5 100644 --- a/tests/formulations/test_formulations_factory.py +++ b/tests/formulations/test_formulations_factory.py @@ -55,7 +55,7 @@ def test_create_with_wrong_formulation_name(factory) -> None: ImportError, match=re.escape( "The class foo is not available; " - "the available ones are: BiLevel, DisciplinaryOpt, IDF, MDF." + "the available ones are: BiLevel, BiLevelBCD, DisciplinaryOpt, IDF, MDF." ), ): factory.create("foo", None, None, None)