Source code for gemseo.problems.scalable.data_driven.problem

# -*- 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:  Matthias De Lozzo
#    OTHER AUTHORS   - MACROSCOPIC CHANGES
"""
Scalable MDO problem
====================

This module implements the concept of scalable problem by means of the
:class:`.ScalableProblem` class.

Given

- a MDO scenario based on a set of sampled disciplines
  with a particular problem dimension,
- a new problem dimension (= number of inputs and outputs),

a scalable problem:

1. makes each discipline scalable based on the new problem dimension,
2. creates the corresponding MDO scenario.

Then, this MDO scenario can be executed and post-processed.

We can repeat this tasks for different sizes of variables
and compare the scalability, which is the dependence of the scenario results
on the problem dimension.

.. seealso:: :class:`.MDODiscipline`, :class:`.ScalableDiscipline`
   and :class:`.Scenario`
"""
from __future__ import division, unicode_literals

import logging
import os
from copy import deepcopy

from numpy import array, ones, random, where, zeros

from gemseo.algos.design_space import DesignSpace
from gemseo.api import (
    create_design_space,
    create_scenario,
    generate_coupling_graph,
    generate_n2_plot,
    get_all_inputs,
)
from gemseo.core.coupling_structure import MDOCouplingStructure
from gemseo.mda.mda_factory import MDAFactory
from gemseo.problems.scalable.data_driven.discipline import ScalableDiscipline
from gemseo.utils.string_tools import MultiLineString

LOGGER = logging.getLogger(__name__)


[docs]class ScalableProblem(object): """Scalable problem.""" def __init__( self, datasets, design_variables, objective_function, eq_constraints=None, ineq_constraints=None, maximize_objective=False, sizes=None, **parameters ): """Constructor. :param list(Dataset) datasets: disciplinary datasets. :param list(str) design_variables: list of design variable names :param str objective_function: objective function :param list(str) eq_constraints: equality constraints. Default: None. :param list(str) eq_constraints: inequality constraints. Default: None. :param bool maximize_objective: maximize objective. Default: False. :param dict sizes: sizes of input and output variables. If None, use the original sizes. Default: None. :param parameters: optional parameters for the scalable model. """ self.disciplines = [dataset.name for dataset in datasets] self.data = {dataset.name: dataset for dataset in datasets} self.inputs = { dataset.name: dataset.get_names(dataset.INPUT_GROUP) for dataset in datasets } self.outputs = { dataset.name: dataset.get_names(dataset.OUTPUT_GROUP) for dataset in datasets } self.varsizes = {} for dataset in datasets: self.varsizes.update(dataset.sizes) self.design_variables = design_variables self.objective_function = objective_function self.ineq_constraints = ineq_constraints self.eq_constraints = eq_constraints self.maximize_objective = maximize_objective self.scaled_disciplines = [] self.scaled_sizes = {} self._build_scalable_disciplines(sizes, **parameters) self.scenario = None def __str__(self): """String representation of information about the scalable problem. :return: scalable problem description :rtype: str """ disciplines = ", ".join(self.disciplines) design_variables = None if self.design_variables is not None: design_variables = ", ".join(self.design_variables) ineq_constraints = None if self.ineq_constraints is not None: ineq_constraints = ", ".join(self.ineq_constraints) eq_constraints = None if self.eq_constraints is not None: eq_constraints = ", ".join(self.eq_constraints) sizes = [ name + " ({})".format(size) for name, size in self.scaled_sizes.items() ] sizes = ", ".join(sizes) optimize = "maximize" if self.maximize_objective else "minimize" msg = MultiLineString() msg.add("MDO problem") msg.indent() msg.add("Disciplines: {}", disciplines) msg.add("Design variables: {}", design_variables) msg.add("Objective function: {} (to {})", self.objective_function, optimize) msg.add("Inequality constraints: {}", ineq_constraints) msg.add("Equality constraints: {}", eq_constraints) msg.add("Sizes: {}", sizes) return str(msg)
[docs] def plot_n2_chart(self, save=True, show=False): """Plot a N2 chart. :param bool save: save plot. Default: True. :param bool show: show plot. Default: False. """ generate_n2_plot(self.scaled_disciplines, save=save, show=show)
[docs] def plot_coupling_graph(self): """Plot a coupling graph.""" generate_coupling_graph(self.scaled_disciplines)
[docs] def plot_1d_interpolations( self, save=True, show=False, step=0.01, varnames=None, directory=".", png=False ): """Plot 1d interpolations. :param bool save: save plot. Default: True. :param bool show: show plot. Default: False. :param bool step: Step to evaluate the 1d interpolation function Default: 0.01. :param list(str) varnames: names of the variable to plot; if None, all variables are plotted. Default: None. :param str directory: directory path. Default: '.'. :param bool png: if True, the file format is PNG. Otherwise, use PDF. Default: False. """ if not os.path.exists(directory): os.mkdir(directory) allfnames = [] for scalable_discipline in self.scaled_disciplines: func = scalable_discipline.scalable_model.plot_1d_interpolations fnames = func(save, show, step, varnames, directory, png) allfnames = allfnames + [os.path.join(directory, fname) for fname in fnames] return allfnames
[docs] def plot_dependencies(self, save=True, show=False, directory="."): """Plot dependency matrices. :param bool save: save plot (default: True) :param bool show: show plot (default: False) :param str directory: directory path (default: '.') """ fnames = [] for scalable_discipline in self.scaled_disciplines: scalable_model = scalable_discipline.scalable_model plot_dependency = scalable_model.plot_dependency fname = plot_dependency( add_levels=True, save=save, show=show, directory=directory ) fnames.append(fname) return fnames
def _build_scalable_disciplines(self, sizes=None, **parameters): """Build scalable disciplines. :param dict sizes: dictionary whose keys are variable names and variables sizes. :param parameters: options. """ copied_parameters = deepcopy(parameters) for disc in self.disciplines: varnames = self.inputs[disc] + self.outputs[disc] sizes = sizes or {} new_varsizes = { varname: sizes.get(varname, self.varsizes[varname]) for varname in varnames } if "group_dep" in parameters: copied_parameters["group_dep"] = parameters["group_dep"][disc] if "fill_factor" in parameters: copied_parameters["fill_factor"] = parameters["fill_factor"][disc] self.scaled_disciplines.append( ScalableDiscipline( "ScalableDiagonalModel", self.data[disc], new_varsizes, **copied_parameters ) ) self.scaled_sizes.update(deepcopy(new_varsizes))
[docs] def create_scenario( self, formulation="DisciplinaryOpt", scenario_type="MDO", start_at_equilibrium=False, active_probability=0.1, feasibility_level=0.5, **options ): """Create MDO scenario from the scalable disciplines. :param str formulation: MDO formulation. Default: 'DisciplinaryOpt'. :param str scenario_type: type of scenario ('MDO' or 'DOE'). Default: 'MDO'. :param bool start_at_equilibrium: start at equilibrium using a preliminary MDA. Default: True. :param float active_probability: probability to set the inequality constraints as active at initial step of the optimization. Default: 0.1. :param float feasibility_level: offset of satisfaction for inequality constraints. Default: 0.5. :param options: formulation options. """ equilibrium = None if start_at_equilibrium: equilibrium = self.__get_equilibrium() disciplines = self.scaled_disciplines design_space = self._create_design_space(disciplines, formulation) max_obj = self.maximize_objective if formulation == "BiLevel": self.scenario = self._create_bilevel_scenario(disciplines, **options) else: self.scenario = create_scenario( disciplines, formulation, self.objective_function, deepcopy(design_space), scenario_type=scenario_type, maximize_objective=max_obj, **options ) equilibrium = {} if not isinstance(equilibrium, dict) else equilibrium self.__add_ineq_constraints(active_probability, feasibility_level, equilibrium) self.__add_eq_constraints(equilibrium) return self.scenario
def _create_bilevel_scenario(self, disciplines, **sub_scenario_options): """Create a bilevel scenario from disciplines. :param list(MDODiscipline) disciplines: list of MDODiscipline """ cpl_structure = MDOCouplingStructure(disciplines) st_cpl_disciplines = cpl_structure.strongly_coupled_disciplines() wk_cpl_disciplines = cpl_structure.weakly_coupled_disciplines() obj = self.objective_function max_obj = self.maximize_objective # Construction of the subsystem scenarios sub_scenarios = [] sub_inputs = [] for discipline in st_cpl_disciplines: cplt_disciplines = list(set(disciplines) - {discipline}) sub_disciplines = [discipline] + wk_cpl_disciplines design_space = DesignSpace() inputs = get_all_inputs([discipline]) all_inputs = get_all_inputs(cplt_disciplines) inputs = list(set(inputs) - set(all_inputs)) sub_inputs += inputs for name in inputs: design_space.add_variable( name, self.scaled_sizes[name], "float", 0.0, 1.0, 0.5 ) sub_scenarios.append( create_scenario( disciplines=sub_disciplines, formulation="DisciplinaryOpt", objective_name=obj, design_space=design_space, maximize_objective=max_obj, ) ) sub_scenarios[-1].default_inputs = sub_scenario_options # Construction of the system scenario all_inputs = get_all_inputs(disciplines) inputs = list(set(all_inputs) - set(sub_inputs)) design_space = DesignSpace() for name in inputs: design_space.add_variable( name, self.scaled_sizes[name], "float", 0.0, 1.0, 0.5 ) sub_disciplines = sub_scenarios + wk_cpl_disciplines system_scenario = create_scenario( disciplines=sub_disciplines, formulation="BiLevel", objective_name=obj, design_space=design_space, maximize_objective=max_obj, mda_name="MDAJacobi", tolerance=1e-8, ) return system_scenario def _create_design_space(self, disciplines=None, formulation="DisciplinaryOpt"): """Create a design space into the unit hypercube. :param list(MDODiscipline) disciplines: list of MDODiscipline :param str formulation: MDO formulation (default: 'DisciplinaryOpt') """ design_space = create_design_space() for varname in self.design_variables: size = self.scaled_sizes[varname] l_b = zeros(size) u_b = ones(size) value = 0.5 + zeros(size) design_space.add_variable(varname, size, "float", l_b, u_b, value) if formulation == "IDF": coupling_structure = MDOCouplingStructure(disciplines) all_couplings = set(coupling_structure.get_all_couplings()) for varname in all_couplings: size = self.scaled_sizes[varname] design_space.add_variable( varname, size, "float", zeros(size), ones(size), 0.5 + zeros(size) ) return design_space def __get_equilibrium(self, mda_name="MDAJacobi", **options): """Get the equilibrium point from a MDA method. :param str mda_name: MDA name (default: 'MDAJacobi') :return: equilibrium point :rtype: dict """ LOGGER.info("Build a preliminary MDA to start at equilibrium") factory = MDAFactory() mda = factory.create(mda_name, self.scaled_disciplines, **options) if len(mda.strong_couplings) == 0: mda = factory.create("MDAQuasiNewton", self.scaled_disciplines, **options) return mda.execute() def __add_ineq_constraints( self, active_probability, feasibility_level, equilibrium ): """Add inequality constraints. :param float active_probability: probability to set the inequality constraints as active at initial step of the optimization :param float feasibility_level: offset of satisfaction for inequality constraints :param dict equilibrium: starting point at equilibrium """ if not hasattr(feasibility_level, "__len__"): feasibility_level = { constraint: feasibility_level for constraint in self.ineq_constraints } for constraint, alphai in feasibility_level.items(): if constraint in list(equilibrium.keys()): sample = random.rand(len(equilibrium[constraint])) val = equilibrium[constraint] taui = where( sample < active_probability, val, alphai + (1 - alphai) * val ) else: taui = 0.0 self.scenario.add_constraint(constraint, "ineq", value=taui) def __add_eq_constraints(self, equilibrium): """Add equality constraints. :param dict equilibrium: starting point at equilibrium """ for constraint in self.eq_constraints: cstr_value = equilibrium.get(constraint, array([0.0]))[0] self.scenario.add_constraint(constraint, "eq", value=cstr_value)
[docs] def exec_time(self, do_sum=True): """Get total execution time per discipline. :param bool do_sum: sum over disciplines (default: True) :return: execution time :rtype: list(float) or float """ exec_time = [discipline.exec_time for discipline in self.scenario.disciplines] if do_sum: exec_time = sum(exec_time) return exec_time
@property def n_calls_top_level(self): """Get number of top level disciplinary calls per discipline. :return: number of top level disciplinary calls per discipline :rtype: list(int) or int """ disciplines = self.scenario.formulation.get_top_level_disc() n_calls = {discipline.name: discipline.n_calls for discipline in disciplines} return n_calls @property def n_calls_linearize_top_level(self): """Get number of top level disciplinary calls per discipline. :return: number of top level disciplinary calls per discipline :rtype: list(int) or int """ disciplines = self.scenario.formulation.get_top_level_disc() n_calls = { discipline.name: discipline.n_calls_linearize for discipline in disciplines } return n_calls @property def n_calls(self): """Get number of disciplinary calls per discipline. :return: number of disciplinary calls per discipline :rtype: list(int) or int """ n_calls = { discipline.name: discipline.n_calls for discipline in self.scenario.disciplines } return n_calls @property def n_calls_linearize(self): """Get number of disciplinary calls per discipline. :return: number of disciplinary calls per discipline :rtype: list(int) or int """ tmp = { discipline.name: discipline.n_calls_linearize for discipline in self.scenario.disciplines } return tmp @property def status(self): """Get the status of the scenario.""" return self.scenario.optimization_result.status @property def is_feasible(self): """Get the feasibility property of the scenario.""" return self.scenario.optimization_result.is_feasible