diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d3b3f40..0fc6bcc 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -37,7 +37,7 @@ jobs: python -m pip install --upgrade pip setuptools pip install -r requirements-tests.txt pip install -r requirements-assignment.txt - pip install -e .[arch_sbo,rocket] + pip install -e .[arch_sbo,rocket,uncertainty] - name: Test with pytest run: pytest -v sb_arch_opt --durations=10 diff --git a/docs/algo/stochastic.md b/docs/algo/stochastic.md new file mode 100644 index 0000000..9a0e833 --- /dev/null +++ b/docs/algo/stochastic.md @@ -0,0 +1,159 @@ +# Stochastic Optimization + +Architecture evaluations often depend on parameters that are not chosen by the optimizer: manufacturing tolerances, +operating conditions, material properties. Optimizing against their nominal values yields a design that is optimal on +paper, but that might perform poorly in reality. Stochastic optimization +instead optimizes a statistic of the response distribution, and therefore looks for a design that is insensitive to +the uncertainties present in the system. + +SBArchOpt supports this by evaluating each design point for many realizations of the uncertain parameters, and reducing +the resulting response distributions to the single values the optimizer sees. +The problem is defined with a `StochasticArchOptProblem`, a subclass of`ArchOptProblemBase`, so +a stochastic problem can be solved with pymoo, ArchSBO, or any of the connected optimization frameworks. +Each design points costs additional `n_evaluations` evaluations of the underlying model, +which makes Surrogate-Based Optimization (SBO) especially attractive to reduce the number of design points . + +To facilitate uncertainty quantification, sampling and distribution storage, SBArchOpt uses +the objects defined in [OpenTURNS](https://openturns.github.io/). + +Three things are defined on top of the design variables: + +- A `StochasticParameterSpace`: the joint distribution of the uncertain parameters, assumed to be independent. They + influence the evaluation but are not chosen by the optimizer, and therefore are NOT design variables. +- A `UQMethod` (`MonteCarlo` or `PolynomialChaos`): handles the uncertainty propagation, + and turning the responses into statistics. It is given a budget of `n_evaluations` expensive evaluations per design point. +- A `Scalarization` per objective and constraint (`Mean`, `Margin` or `Quantile`): scalarizes the output distribution to + a single number that optimizer sees. If not specified, responses default to `Mean`. + +## Installation + +``` +pip install sb-arch-opt[uncertainty] +``` + +## Usage + +[API Reference](../api/robust.md) + +A stochastic problem is a subclass of `StochasticArchOptProblem`, which extends `ArchOptProblemBase`. Compared to a +deterministic problem, evaluation is implemented in `_arch_evaluate_sample` instead of `_arch_evaluate`. This method +is called once per realization, receives it in `parameters`, and writes the responses for that +realization. The class automatically handles the loop over the sample realizations, output distribution construction +and the scalarization to `f`, `g` and `h` for the the optimizer. + +Since the result is an `ArchOptProblemBase`, the problem is then handed to an algorithm like any other. See the +[Stochastic Optimization Tutorial](../tutorial_stochastic.ipynb) for a worked example comparing the robust Pareto front +found by NSGA2 and by SBO to the deterministic one. + +### Implementing a Stochastic Problem + +Implement `_arch_evaluate_sample` rather than `_arch_evaluate`; everything else works as described in +[implementing an architecture optimization problem](../index.md). Uncertain parameters are stored in the parameter space, +and their samples are reused for all design vectors in `x`. + +```python +import openturns as ot +from pymoo.core.variable import Real +from sb_arch_opt.stochastic_problem import StochasticArchOptProblem +from sb_arch_opt.uncertainty import MonteCarlo, Mean, StochasticParameter, \ + StochasticParameterSpace + + +class MyStochasticProblem(StochasticArchOptProblem): + + def __init__(self): + # The uncertain parameters: not design variables + param_space = StochasticParameterSpace([ + StochasticParameter('u', ot.Normal(1., .05)), + ]) + + super().__init__( + [Real(bounds=(-2., 2.))], + param_space=param_space, + uq_method=MonteCarlo(n_evaluations=1000, seed=42), + n_obj=1, obj_scalar=[Mean()], + ) + + def _is_conditionally_active(self): + return [False] + + def _correct_x(self, x, is_active): + pass + + def _arch_evaluate_sample(self, x, is_active, f_out, g_out, h_out, parameters, + *args, **kwargs): + # Called once per realization; parameters contains one realization of the parameter space + f_out[:, 0] = (parameters[0] - x[:, 0])**2 +``` + +All design points share the same sample realizations (common random numbers) during the evaluation, +ensuring the outputs remain directly comparable. For reproducibility during development, +the UQ method can be run with a fixed seed so that identical realizations are generated across independent runs. + +A failed evaluation (NaN) in any realization fails the whole design point: the response has no distribution and the +optimizer sees NaN, which is how SBArchOpt treats hidden constraint violations elsewhere. + +### Uncertainty Quantification Methods + +Both methods draw samples from the parameter space using Latin Hypercube Sampling (LHS) +and use the same number of expensive model evaluations per design point. +They differ in how the statistics are obtained: Monte Carlo computes them directly from the sampled model outputs, +while PCE first fits a polynomial chaos expansion to those outputs and derives the statistics by performing +Monte Carlo on the trained cheap surrogate model. In general Monte Carlo requires more simulation code executions +to compute the statistics compared to PCE. More UQ methods can be implemented by creating a new subclass of `UQMethod`. + +- `MonteCarlo`: the statistics are computed from the `n_evaluations` responses directly. It is the simplest method to + implement, but slowest to converge. +- `PolynomialChaos`: a polynomial chaos expansion of each response as a function of the uncertain parameters is fitted + from the `n_evaluations` responses, and the statistics are then taken from `n_metamodel_samples` cheap evaluations. + Sampling the trained surrogate model is much cheaper than running the real simulation code, + which makes this method attractive for expensive evaluation codes. The PCE surrogate model is trained + with the least-squares method that requires at least `P` samples, but it is commonly recommended to provide + between `2P` and `3P`, where `P = (n+d)! / (n! d!)` is the number of polynomial terms for `n` parameters and polynomial + degree `d` (Blatman, 2009). + +### Scalarization + +Each objective and constraint is reduced to one value by the selected `Scalarization`. +The resulting `float` is passed to the chosen optimization algorithm. + +- `Mean`: the expected value of the response --> `min(E[F(x)])` +- `Margin`: `min(E[F(x)] + k*sigma[F(x)])`, assuming a Gaussian response. The default `k=1.645` is the one-sided 95% + interval of a normal distribution; use `direction=1` for a response that is maximized, so that the margin is + subtracted rather than added. +- `Quantile`: `min(F_q(x))`, the `q` quantile of the response. Makes no assumption about the shape of the distribution. + +Custom scalarizations are implemented by extending `Scalarization` and implementing `scalarize`, which receives the +whole `StochasticOutput` and can therefore use anything the fitted distribution offers. + +### Sampled Responses + +After each evaluation the fitted response distributions are available in the pymoo output dictionary as +`StochasticOutput` objects, in `out['f_stochastic']`, `out['g_stochastic']` and `out['h_stochastic']`. A response that +turns out not to depend on the uncertain parameters is a plain value instead, and a response that failed is NaN. + +```python +import numpy as np + +problem = ... # Subclass of StochasticArchOptProblem + +out = problem.evaluate(np.array([[.5]]), return_as_dictionary=True) +output = out['f_stochastic'][0, 0] + +print(output.mean, output.std, output.quantile(.95)) +``` + +For polynomial chaos each output also carries the fitted expansion in `method_results`, from which Sobol sensitivity +indices are available to see which uncertain parameters drive the spread of a response: + +```python +import openturns as ot + +sobol = ot.FunctionalChaosSobolIndices(output.method_results) +first_order_u0 = sobol.getSobolIndex(0) +``` + +## References +- Blatman, G. (2009). *Adaptive sparse polynomial chaos expansions for uncertainty + propagation and sensitivity analysis*. PhD thesis, Université Blaise Pascal, + Clermont-Ferrand, France. \ No newline at end of file diff --git a/docs/api/robust.md b/docs/api/robust.md new file mode 100644 index 0000000..0cd3152 --- /dev/null +++ b/docs/api/robust.md @@ -0,0 +1,44 @@ +# Stochastic Problem Definition + +Stochastic optimization uses [OpenTURNS](https://openturns.github.io/) for the uncertainty propagation, +which is an optional dependency: + +``` +pip install sb-arch-opt[uncertainty] +``` + +::: sb_arch_opt.stochastic_problem.StochasticArchOptProblem + handler: python + +::: sb_arch_opt.uncertainty.EvaluationOutput + handler: python + +::: sb_arch_opt.uncertainty.Scalarization + handler: python + +::: sb_arch_opt.uncertainty.Mean + handler: python + +::: sb_arch_opt.uncertainty.Margin + handler: python + +::: sb_arch_opt.uncertainty.Quantile + handler: python + +::: sb_arch_opt.uncertainty.StochasticParameter + handler: python + +::: sb_arch_opt.uncertainty.StochasticParameterSpace + handler: python + +::: sb_arch_opt.uncertainty.StochasticOutput + handler: python + +::: sb_arch_opt.uncertainty.UQMethod + handler: python + +::: sb_arch_opt.uncertainty.MonteCarlo + handler: python + +::: sb_arch_opt.uncertainty.PolynomialChaos + handler: python \ No newline at end of file diff --git a/docs/index.md b/docs/index.md index d101b7b..8cc5a6d 100644 --- a/docs/index.md +++ b/docs/index.md @@ -64,6 +64,7 @@ DOI: [10.2514/6.2024-4401](https://arc.aiaa.org/doi/10.2514/6.2024-4401) See also the tutorial(s): - [SBArchOpt Tutorial](tutorial.ipynb): optimization, implementing new problems +- [Stochastic Optimization Tutorial](tutorial_stochastic.ipynb): optimizing a problem with stochastic parameters ### Connecting to Optimization Frameworks @@ -86,6 +87,11 @@ Following proprietary frameworks are also integrated: Optimization framework dependencies (except pymoo) are always optional! Therefore, when installing SBArchOpt, the interfaces to optimization frameworks are probably not working yet. +### Stochastic Optimization + +Problems with uncertain parameters can be optimized for a statistic of the response distribution rather than for its +nominal value, see: [Stochastic Optimization](algo/stochastic.md) + ### Test Problems For an overview of the available test problems and how to use them, see: diff --git a/docs/tutorial_stochastic.ipynb b/docs/tutorial_stochastic.ipynb new file mode 100644 index 0000000..397cccd --- /dev/null +++ b/docs/tutorial_stochastic.ipynb @@ -0,0 +1,199 @@ +{ + "cells": [ + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "# Stochastic Optimization Tutorial\n", + "\n", + "Try it out yourself:\n", + "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/jbussemaker/SBArchOpt/HEAD?labpath=docs%2Ftutorial_stochastic.ipynb)\n", + "\n", + "This tutorial takes `MDMOHimmelblau`, the mixed-discrete multi-objective Himmelblau test problem, and turns it into a stochastic problem by assuming that the continuous design variable cannot be measured accurately. This uncertainty is modeled by adding a perturbation to the variable, represented by a zero-mean Gaussian distribution. The tutorial then compares the robust Pareto fronts found by NSGA2 and SBO with the deterministic one. For more background information, refer to the [documentation](../algo/stochastic).\n", + "\n", + "**Note**: make sure you have installed the `uncertainty` and `arch_sbo` dependencies:\n", + "```\n", + "pip install sb-arch-opt[uncertainty,arch_sbo]\n", + "```" + ], + "id": "cb2dfcb356ad626b" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "import warnings\n", + "warnings.filterwarnings('ignore')" + ], + "id": "853186ccc64b311b" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## The Stochastic Problem Class\n", + "\n", + "Stochastic optimization problems in SBArchOpt are defined by extending the `StochasticArchOptProblem` class, which itself extends `ArchOptProblemBase`. Three things are defined on top of the design variables:\n", + "- A `StochasticParameterSpace`: the joint distribution of the uncertain parameters, assumed to be independent. They\n", + " influence the evaluation but are not chosen by the optimizer, and therefore are NOT design variables.\n", + "- A `UQMethod` (`MonteCarlo` or `PolynomialChaos`): handles the uncertainty propagation,\n", + " and turning the responses into statistics. It is given a budget of `n_evaluations` expensive evaluations per design point.\n", + "- A `Scalarization` per objective and constraint (`Mean`, `Margin` or `Quantile`): scalarizes the output distribution to\n", + " a single number that optimizer sees. If not specified, responses default to `Mean`.\n", + "\n", + "Evaluation is implemented in `_arch_evaluate_sample` instead of `_arch_evaluate`. This function is called once per realization of the uncertain parameters, receives it in `parameters`, and writes the responses for that realization only.\n", + "\n", + "Below, `MDMOHimmelblau` is wrapped in such a class. Its continuous design variable is realized with a normally-distributed deviation from its nominal value, and the underlying problem is then evaluated at the realized design vector. Since the design space of the underlying problem is reused, the design variables and the correction of design vectors are inherited as-is and `_correct_x` does not need to be implemented." + ], + "id": "56d6dab3b1e5f7b" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "import numpy as np\n", + "import openturns as ot\n", + "from sb_arch_opt.problems.md_mo import MDMOHimmelblau\n", + "from sb_arch_opt.stochastic_problem import StochasticArchOptProblem\n", + "from sb_arch_opt.uncertainty import PolynomialChaos, StochasticParameter, StochasticParameterSpace\n", + "\n", + "\n", + "class StochasticMDMOHimmelblau(StochasticArchOptProblem):\n", + " \"\"\"MDMOHimmelblau where the continuous design variable is realized with a manufacturing tolerance\"\"\"\n", + "\n", + " def __init__(self, std=.3, n=10, seed=42, obj_scalar=None):\n", + " self._problem = MDMOHimmelblau()\n", + " self.std = std\n", + "\n", + " # One uncertain parameter per continuous design variable: the deviation from its nominal value\n", + " self._i_cont = np.where(~self._problem.is_discrete_mask)[0]\n", + " param_space = StochasticParameterSpace([\n", + " StochasticParameter(f'dx{i}', ot.Normal(0., std)) for i in self._i_cont])\n", + "\n", + " super().__init__(\n", + " self._problem.design_space, # Same design space, so the same variables and correction\n", + " param_space=param_space,\n", + " uq_method=PolynomialChaos(n_evaluations=n, seed=seed, degree=3, n_metamodel_samples=1000),\n", + " n_obj=self._problem.n_obj, obj_scalar=obj_scalar,\n", + " )\n", + "\n", + " def _arch_evaluate_sample(self, x, is_active, f_out, g_out, h_out, parameters, *args, **kwargs):\n", + " # Called once per realization: apply the deviations and evaluate the underlying problem\n", + " x_realized = x.copy()\n", + " x_realized[:, self._i_cont] += parameters\n", + "\n", + " f_out[:, :] = self._problem.evaluate(x_realized, return_as_dictionary=True)['F']\n", + "\n", + " def __repr__(self):\n", + " return f'{self.__class__.__name__}(std={self.std}, n={self.uq_method.n_evaluations})'" + ], + "id": "7c735da49e170b2b" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "# Instantiate the problem class and print some statistics\n", + "# Sample size for degree d PCE regression: N = 2P to 3P, where\n", + "# P = (n+d)! / (n! d!) is the number of polynomial terms for n parameters.\n", + "problem = StochasticMDMOHimmelblau(std=.3, n=12)\n", + "problem.print_stats()" + ], + "id": "2db9a8d421e3f270" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "The resulting problem has the same design variables and objectives as the deterministic one, only the evaluation has changed. Note that one design point now costs 12 evaluations of the underlying problem to train the PCE surrogate. Then the cheap PCE model is evaluated with 1000 samples to construct the distribution. Note that the uncertainty propagation requires overall more model evaluation which makes Surrogate-Based Optimization (SBO) especially attractive for stochastic problems.", + "id": "d0ecd1d63f400dc5" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Solving with NSGA2\n", + "\n", + "The stochastic problem is an `ArchOptProblemBase` subclass, so the algorithms in SBArchOpt can be used without modification. The code below runs the pre-configured NSGA2 and compares the resulting robust Pareto front to the known deterministic Pareto front of `MDMOHimmelblau`." + ], + "id": "d58823d06acbc6ed" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "from pymoo.optimize import minimize\n", + "from sb_arch_opt.algo.pymoo_interface import plot, get_nsga2\n", + "\n", + "# Create the pre-configured NSGA2 algorithm\n", + "nsga2 = get_nsga2(pop_size=100)\n", + "\n", + "# Run the optimization: 20 generations (20 x 100 = 2000 design points, each with 50 realizations)\n", + "# NOTE: remove \"seed\" argument when using in production!\n", + "result_nsga2 = minimize(problem, nsga2, termination=('n_gen', 20), seed=42)\n", + "\n", + "# Compare the robust Pareto front to the deterministic one\n", + "pf_deterministic = MDMOHimmelblau().pareto_front()\n", + "plot(pf_deterministic, result_nsga2.opt.get('F'),\n", + " labels=['Deterministic Pareto front', 'NSGA2 robust (2000 design points)'])" + ], + "id": "762c1b2ad0127509" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "The robust Pareto front lies above the deterministic one, which clearly shows that the perturbation in design vector has effect on the objectives. The deterministic front is therefore not achievable in practice, and optimizing towards it would select designs that are more sensitive to the perturbations.", + "id": "a4f99e6120985213" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Solving with SBO\n", + "\n", + "Because every design point requires 12 evaluations for each sample realization, the number of design points for evaluation should be kept low. Surrogate-Based Optimization (SBO) in general requires fewer evaluations than genetic algorithms to find the Pareto front, therefore they are well suited for the stochastic optimization problems." + ], + "id": "efe12d81968cda69" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "from sb_arch_opt.algo.arch_sbo.api import get_arch_sbo_gp\n", + "\n", + "# Create the algorithm (with an initial database size of 20)\n", + "n_init = 20\n", + "sbo = get_arch_sbo_gp(problem, n_parallel=4, init_size=n_init)\n", + "\n", + "# Run the optimization: 30 infill points (20 + 30 = 50 design points)\n", + "# NOTE: remove \"seed\" argument when using in production!\n", + "n_infill = 30\n", + "result_sbo = minimize(problem, sbo, termination=('n_eval', n_init+n_infill), seed=42)\n", + "\n", + "# Compare both robust Pareto fronts to the deterministic one\n", + "plot(pf_deterministic, result_nsga2.opt.get('F'), result_sbo.opt.get('F'),\n", + " labels=['Deterministic Pareto front', 'NSGA2 robust (2000 design points)',\n", + " 'SBO robust (50 design points)'])" + ], + "id": "68b1694f096c0818" + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "SBO reaches the same robust Pareto front as NSGA2 using 40 times fewer design points. Given that for each design points 12 evaluations are performed to propagate uncertainty, the computational budget savings are tremendous compared with NSGA-2 (overall 480 less evaluations).", + "id": "d4f0ff16942cc3c4" + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/mkdocs.yml b/mkdocs.yml index 1a08fc5..d98247d 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -26,11 +26,14 @@ nav: - 'HEBO': 'algo/hebo.md' - 'SEGOMOE': 'algo/segomoe.md' - 'SMARTy': 'algo/smarty.md' + - 'Stochastic Optimization': 'algo/stochastic.md' - Tutorials: - SBArchOpt: tutorial.ipynb + - Stochastic Optimization: tutorial_stochastic.ipynb - 'Test Problems': 'test_problems.md' - API Reference: - Problem Definition: 'api/problem.md' + - Stochastic Optimization: 'api/robust.md' - pymoo: 'api/pymoo.md' - ArchSBO: 'api/arch_sbo.md' - 'BoTorch (Ax)': 'api/botorch.md' diff --git a/sb_arch_opt/__init__.py b/sb_arch_opt/__init__.py index 4574cc8..218f431 100644 --- a/sb_arch_opt/__init__.py +++ b/sb_arch_opt/__init__.py @@ -1 +1 @@ -__version__ = '1.6.3' +__version__ = '1.7' diff --git a/sb_arch_opt/stochastic_problem.py b/sb_arch_opt/stochastic_problem.py new file mode 100644 index 0000000..dc31adb --- /dev/null +++ b/sb_arch_opt/stochastic_problem.py @@ -0,0 +1,209 @@ +""" +MIT License + +Copyright: (c) 2024, Deutsches Zentrum fuer Luft- und Raumfahrt e.V. +Contact: jasper.bussemaker@dlr.de + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +""" +import numpy as np +from typing import Union, List, Optional +from pymoo.core.variable import Variable +from sb_arch_opt.design_space import ArchDesignSpace +from sb_arch_opt.problem import ArchOptProblemBase + +try: + from sb_arch_opt.uncertainty import StochasticParameterSpace, UQMethod, Scalarization, Mean + HAS_UNCERTAINTY = True + +except ImportError: + HAS_UNCERTAINTY = False + +__all__ = ['StochasticArchOptProblem', 'HAS_UNCERTAINTY', 'check_dependency'] + + +def check_dependency(): + if not HAS_UNCERTAINTY: + raise ImportError('Looks like SBArchOpt uncertainty package is not installed! Run: pip install sb-arch-opt[uncertainty]') + + +class StochasticArchOptProblem(ArchOptProblemBase): + """ + Base class for a stochastic (robust) architecture optimization problem. It extends the ArchOptProblemBase class + with support for stochastic optimization problems. + + Stochastic problem evaluates design vector for many realizations of the uncertain parameters and + reduces the resulting distribution of each response to the single value the optimizer sees. + Three things are therefore needed on top of the design variables: + + - `param_space`: a `StochasticParameterSpace`, the joint distribution of the uncertain parameters. These are + quantities that influence the evaluation but are not chosen by the optimizer, so they are NOT design + variables. + - `uq_method`: a `UQMethod` (`MonteCarlo`, `PolynomialChaos`) that decides which realizations to evaluate and + how to turn the responses into statistics. It is given a budget of `n_evaluations` expensive evaluations per + design point. + - a `Scalarization` per response, below. + + Requires specifying a type of optimization problem for each objective and constraint by providing a child + instance of Scalarization object: + + - Mean: Minimize the expectation of the objective or constraint function + --> for example min(E[F(x)]) + - Margin: Gaussian output distribution expected. Minimize for the objective or constraint function for a given + confidence interval. + --> min(E[F(x)] + k*sigma[F(x)]); The default k=1.645 is the one-sided 95% interval of a normal distribution. + - Quantile: Makes no distributional assumption. Minimize the q quantile of the objective or constraint + function --> min(F_q(x)), with q=0.95 by default. + + When type of optimization problem is not provided, the response defaults to `Mean()`. + + Implement `_arch_evaluate_sample`: it is called once per realization. All design points in a batch see the same + realizations (common random numbers), so they stay comparable to each other and to a surrogate fitted through + them. + + A failed evaluation (NaN) in any realization fails the whole design point: `scalarize` returns NaN, + which is how SBArchOpt treats hidden-constraint violations elsewhere. + + After each evaluation the fitted response distributions are stored in the pymoo output dictionary as + `StochasticOutput` objects, in `out['f_stochastic']`, `out['g_stochastic']` and `out['h_stochastic']. + For polynomial chaos each output also carries the fitted expansion (`method_results`), from which Sobol + indices can be obtained. If the output turns out to be deterministic after evaluation, a float is stored instead of a + 'StochasticOutput' object. + """ + + def __init__(self, des_vars: Union[List[Variable], ArchDesignSpace], + param_space: StochasticParameterSpace, + uq_method: UQMethod, n_obj=1, n_ieq_constr=0, n_eq_constr=0, + obj_scalar: Optional[List[Scalarization]] = None, + ieq_constr_scalar: Optional[List[Scalarization]] = None, + eq_constr_scalar: Optional[List[Scalarization]] = None, + **kwargs): + + check_dependency() + self.obj_scalar = self.check_scalars(obj_scalar, n_obj) + self.ieq_constr_scalar = self.check_scalars(ieq_constr_scalar, n_ieq_constr) + self.eq_constr_scalar = self.check_scalars(eq_constr_scalar, n_eq_constr) + self.scalars = self.obj_scalar + self.ieq_constr_scalar + self.eq_constr_scalar + + if uq_method is None: + raise ValueError('A UQ method must be specified, e.g. uq_method=MonteCarlo(n_evaluations=100, seed=42)') + if not isinstance(uq_method, UQMethod): + raise ValueError(f'uq_method should be a UQMethod instance, got: {uq_method!r}') + + if param_space is None or param_space.n_parameters == 0: + raise ValueError('param_space must contain valid parameter space') + + self.param_space = param_space + self.uq_method = uq_method + + super().__init__(des_vars, n_obj=n_obj, n_ieq_constr=n_ieq_constr, n_eq_constr=n_eq_constr, **kwargs) + + @staticmethod + def check_scalars(scalars: Optional[List[Scalarization]], n: int) -> List[Scalarization]: + """Default unspecified responses to the expected value, and check the count""" + if scalars is None: + # Mean is the default robust scalar type + return [Mean() for _ in range(n)] + + if len(scalars) != n: + raise ValueError(f'scalars should have {n} Scalarization entries, got {len(scalars)}') + for scalar in scalars: + if not isinstance(scalar, Scalarization): + raise ValueError(f'scalars should contain Scalarization instances, got: {scalar!r}') + return list(scalars) + + def _print_extra_stats(self): + print(f'stochastic : True') + print(f'n_stochastic_params : {self.param_space.n_parameters}') + print(f'uq_method : {self.uq_method}') + print(f'n_uq_samples : {self.uq_method.n_evaluations}') + print(f'obj_scalarize : {self.obj_scalar}') + print(f'ieq_constr_scalarize : {self.ieq_constr_scalar}') + print(f'eq_constr_scalarize : {self.eq_constr_scalar}') + + + def _evaluate(self, x, out, *args, **kwargs): + n = x.shape[0] + f_stoch = np.empty((n, self.n_obj), dtype=object) + g_stoch = np.empty((n, self.n_ieq_constr), dtype=object) + h_stoch = np.empty((n, self.n_eq_constr), dtype=object) + kwargs.update(f_stoch_out=f_stoch, g_stoch_out=g_stoch, h_stoch_out=h_stoch) + + super()._evaluate(x, out, *args, **kwargs) + out['f_stochastic'] = f_stoch + out['g_stochastic'] = g_stoch + out['h_stochastic'] = h_stoch + + + def _arch_evaluate(self, x: np.ndarray, is_active_out: np.ndarray, f_out: np.ndarray, g_out: np.ndarray, + h_out: np.ndarray, *args, f_stoch_out: np.ndarray=None, g_stoch_out: np.ndarray=None, h_stoch_out: np.ndarray=None, **kwargs): + """ + Evaluate architecture for the provided design vectors and samples. + Implement _arch_evaluate_sample to evaluate architecture for single samples realizations. + """ + + if not self.design_space.is_explicit(): + self._correct_x_impute(x, is_active_out) + + # Get samples and include deterministic parameter values for evaluation + samples = self.uq_method.get_samples(self.param_space) + + n_x, n_s = x.shape[0], samples.shape[0] + + f_s = np.zeros((n_x, n_s, self.n_obj))*np.nan + g_s = np.zeros((n_x, n_s, self.n_ieq_constr))*np.nan + h_s = np.zeros((n_x, n_s, self.n_eq_constr))*np.nan + + # Evaluate all design vectors for each realization of the uncertain parameters + for i in range(n_s): + self._arch_evaluate_sample(x, is_active_out, f_s[:, i, :], g_s[:, i, :], h_s[:, i, :], samples[i, :],*args, **kwargs) + + # Evaluate the stochastic result for all the evaluated design vectors and samples + for x_i in range(n_x): + outputs = self.uq_method.process_results(np.concatenate([f_s[x_i], g_s[x_i], h_s[x_i]], axis=1), self.param_space) + + # Reduce the sampled responses of each design point to the values the optimizer sees + n_f, n_g = self.n_obj, self.n_ieq_constr + for f_i, output in enumerate(outputs[:n_f]): + f_stoch_out[x_i, f_i] = output + obj_scalar = self.obj_scalar[f_i] + f_out[x_i, f_i] = obj_scalar.scalarize(output) if not isinstance(output, float) else output + for g_i, output in enumerate(outputs[n_f:n_f+n_g]): + g_stoch_out[x_i, g_i] = output + ieq_constr_scalar = self.ieq_constr_scalar[g_i] + g_out[x_i, g_i] = ieq_constr_scalar.scalarize(output) if not isinstance(output, float) else output + for h_i, output in enumerate(outputs[n_f+n_g:]): + h_stoch_out[x_i, h_i] = output + eq_constr_scalar = self.eq_constr_scalar[h_i] + h_out[x_i, h_i] = eq_constr_scalar.scalarize(output) if not isinstance(output, float) else output + + def _arch_evaluate_sample(self, x: np.ndarray, is_active: np.ndarray, f_out: np.ndarray, g_out: np.ndarray, + h_out: np.ndarray, parameters: np.ndarray, *args, **kwargs): + """ + Implement evaluation and write results in the provided output matrices: + - x (design vectors): discrete variables have integer values, imputed design vectors can be output here (except + if using an explicit design space definition) + - is_active (activeness): vector specifying for each design variable whether it was active or not + - f (objectives): written as a minimization + - g (inequality constraints): written as "<= 0" + - h (equality constraints): written as "= 0" + - sample: realization of single sample + """ + + raise NotImplementedError \ No newline at end of file diff --git a/sb_arch_opt/tests/conftest.py b/sb_arch_opt/tests/conftest.py index e9a6a6f..357f016 100644 --- a/sb_arch_opt/tests/conftest.py +++ b/sb_arch_opt/tests/conftest.py @@ -7,6 +7,23 @@ from pymoo.core.variable import Real, Integer, Choice from pymoo.problems.multi.zdt import ZDT1 +try: + import openturns as ot + from sb_arch_opt.uncertainty import * + from sb_arch_opt.stochastic_problem import * + HAS_UNCERTAINTY = True +except ImportError: + HAS_UNCERTAINTY = False + + class StochasticArchOptProblem: + pass + + +def check_dependency(): + if not HAS_UNCERTAINTY: + raise ImportError( + 'Looks like SBArchOpt uncertainty package is not installed! Run: pip install sb-arch-opt[uncertainty]') + class DummyProblem(ArchOptTestProblemBase): @@ -100,7 +117,125 @@ def failing_problem(): return DummyProblem(fail=True) +def make_space(*distributions) -> 'StochasticParameterSpace': + return StochasticParameterSpace([StochasticParameter(f'u{i}', dist) for i, dist in enumerate(distributions)]) + + +class VectorizedProblem(StochasticArchOptProblem): + """f = (u - x0)^2 + x1^2, evaluated for all design points at once""" + + def __init__(self, n=100, seed=42, uq_method=None, fail=False, **kwargs): + check_dependency() + self.fail = fail + self.seen_parameters = [] + super().__init__([Real(bounds=(-2., 2.)), Real(bounds=(-2., 2.))], + param_space=make_space(ot.Normal(1., .05)), + uq_method=uq_method if uq_method is not None else MonteCarlo(n, seed=seed), + n_obj=1, **kwargs) + + def _is_conditionally_active(self): + return [False, False] + + def _correct_x(self, x, is_active): + pass + + def _arch_evaluate_sample(self, x, is_active, f_out, g_out, h_out, parameters, *args, **kwargs): + self.seen_parameters.append(np.asarray(parameters).copy()) + f_out[:, 0] = (parameters[0] - x[:, 0])**2 + x[:, 1]**2 + if self.fail and parameters[0] > 1.02: # some samples fail to evaluate + f_out[:, 0] = np.nan + + +class HierarchicalProblem(StochasticArchOptProblem): + """Implicit (hierarchical) design space, a constraint, and the per-row evaluation pattern""" + + def __init__(self, n=50, seed=1, **kwargs): + check_dependency() + super().__init__([Choice(options=['electric', 'hybrid']), Real(bounds=(.2, 1.)), Real(bounds=(.1, .4))], + param_space=StochasticParameterSpace([ + StochasticParameter('payload', ot.Normal(2., .3)), + StochasticParameter('bsfc', ot.Normal(.42, .075)), + ]), + uq_method=MonteCarlo(n, seed=seed), n_obj=1, n_ieq_constr=1, **kwargs) + + def _is_conditionally_active(self): + return [False, False, True] + + def _correct_x(self, x, is_active): + is_active[:, 2] = x[:, 0] == 1 # fuel fraction only active for the hybrid branch + + def _arch_evaluate_sample(self, x, is_active, f_out, g_out, h_out, parameters, *args, **kwargs): + payload, bsfc = parameters + for i in range(x.shape[0]): + mass = .5 + x[i, 1] + payload + energy = x[i, 2]*130./(bsfc+.1) if x[i, 0] == 1 else x[i, 1]*40. + f_out[i, 0] = mass**1.5 / energy + g_out[i, 0] = mass - 3.5 + + +class AllResponseKindsProblem(StochasticArchOptProblem): + """One objective, one inequality and one equality constraint, each with its own scalar""" + + def __init__(self, **kwargs): + check_dependency() + kwargs.setdefault('obj_scalar', [Mean()]) + kwargs.setdefault('ieq_constr_scalar', [Mean()]) + kwargs.setdefault('eq_constr_scalar', [Margin(k=3.)]) + super().__init__([Real(bounds=(0., 1.))], param_space=make_space(ot.Normal(1., .2)), + uq_method=MonteCarlo(200, seed=5), + n_obj=1, n_ieq_constr=1, n_eq_constr=1, **kwargs) + + def _is_conditionally_active(self): + return [False] + + def _correct_x(self, x, is_active): + pass + + def _arch_evaluate_sample(self, x, is_active, f_out, g_out, h_out, parameters, *args, **kwargs): + f_out[:, 0] = parameters[0] + x[:, 0] + g_out[:, 0] = parameters[0] - x[:, 0] + h_out[:, 0] = parameters[0] * x[:, 0] + + +class DeterministicResponseProblem(StochasticArchOptProblem): + """Two objectives, of which the second does not depend on the stochastic parameters at all""" + + def __init__(self, uq_method=None, **kwargs): + super().__init__([Real(bounds=(0., 1.))], param_space=make_space(ot.Normal(1., .2)), + uq_method=uq_method or MonteCarlo(50, seed=3), n_obj=2, **kwargs) + + def _is_conditionally_active(self): + return [False] + + def _correct_x(self, x, is_active): + pass + + def _arch_evaluate_sample(self, x, is_active, f_out, g_out, h_out, parameters, *args, **kwargs): + f_out[:, 0] = parameters[0] * x[:, 0] + f_out[:, 1] = 2. * x[:, 0] + + +@pytest.fixture +def stochastic_problem(): + return VectorizedProblem() + + +@pytest.fixture +def hierarchical_problem(): + return HierarchicalProblem() + + +@pytest.fixture +def all_response_kinds_problem(): + return AllResponseKindsProblem() + + +@pytest.fixture +def deterministic_response_problem(): + return DeterministicResponseProblem() + + def pytest_sessionstart(session): from sb_arch_opt.util import _prevent_capture print('PREVENT CAPTURE') - _prevent_capture() + _prevent_capture() \ No newline at end of file diff --git a/sb_arch_opt/tests/test_stochastic_problem.py b/sb_arch_opt/tests/test_stochastic_problem.py new file mode 100644 index 0000000..0917e76 --- /dev/null +++ b/sb_arch_opt/tests/test_stochastic_problem.py @@ -0,0 +1,304 @@ +import pytest +import numpy as np + +try: + import openturns as ot + from sb_arch_opt.uncertainty import * + from sb_arch_opt.stochastic_problem import StochasticArchOptProblem + from sb_arch_opt.tests.conftest import (VectorizedProblem, HierarchicalProblem, + DeterministicResponseProblem, make_space) + +except ImportError: + pytest.skip( + "SBArchOpt uncertainty package is not installed! Run: pip install sb-arch-opt[uncertainty]", + allow_module_level=True, + ) + + +def _output(values): + return StochasticOutput.from_output_samples(ot.Sample(np.array(values, dtype=float).reshape((-1, 1)))) + + +def test_parameter_space(): + space = StochasticParameterSpace([ + StochasticParameter('a', ot.Normal(0., 1.)), + StochasticParameter('b', ot.Uniform(0., 1.)), + ]) + + assert space.n_parameters == 2 + assert space.parameter_names == ['a', 'b'] + assert space.joint_dist.getDimension() == 2 + assert [parameter.name for parameter in space.parameters] == ['a', 'b'] + + for samples in [space.get_random_samples(20), space.get_lhs_samples(20)]: + assert samples.shape == (20, 2) + assert np.all(np.isfinite(samples)) + + +def test_scalars(): + out = _output([1., 2., 3., 4., 5.]) + mean, std = out.mean, out.std + + assert Mean().scalarize(out) == pytest.approx(mean) + assert Quantile(q=.9).scalarize(out) == pytest.approx(out.quantile(.9)) + assert Margin(k=2.).scalarize(out) == pytest.approx(out.margin(2.)) + + # The margin is the conservative value, so which tail that is depends on the direction + assert Margin(k=2.).scalarize(out) == pytest.approx(mean + 2. * std) + assert Margin(k=2., direction=-1).scalarize(out) == pytest.approx(mean + 2. * std) + assert Margin(k=2., direction=1).scalarize(out) == pytest.approx(mean - 2. * std) + + # Each instance carries its own parameters + assert Quantile(q=.5).scalarize(out) != Quantile(q=.95).scalarize(out) + assert Margin(k=1.).scalarize(out) != Margin(k=3.).scalarize(out) + + with pytest.raises(ValueError): + Quantile(q=2.) + + +def test_scalars_penalize_spread(): + wide, narrow = _output([1., 3., 5.]), _output([2.5, 3., 3.5]) + assert wide.mean == pytest.approx(narrow.mean) + + assert Margin(k=2.).scalarize(narrow) < Margin(k=2.).scalarize(wide) # minimized: lower is better + assert Margin(k=2., direction=1).scalarize(narrow) > Margin(k=2., direction=1).scalarize(wide) # maximized + + +@pytest.mark.parametrize('method_name', ['MonteCarlo', 'PolynomialChaos']) +def test_uq_method_samples(method_name): + method_class = {'MonteCarlo': MonteCarlo, 'PolynomialChaos': PolynomialChaos}[method_name] + space = make_space(ot.Normal(0., 1.), ot.Normal(0., 1.)) + method = method_class(n_evaluations=40, seed=42) + + samples = method.get_samples(space) + assert samples.shape == (40, 2) + + # Drawn once and reused, so every design point sees the same realizations (common random numbers) + assert np.all(method.get_samples(space) == samples) + + method.resample() + assert method.get_samples(space).shape == (40, 2) + + with pytest.raises(ValueError): + method.get_samples(None) + + +def test_parameter_realization(): + space = StochasticParameterSpace([ + StochasticParameter('a', ot.Normal(0., 1.)), + StochasticParameter('b', ot.Uniform(2., 4.)), + ]) + samples = MonteCarlo(n_evaluations=5, seed=42).get_samples(space) + + for i in range(5): + parameters = space.param_realization(samples, i) + + assert [parameter.name for parameter in parameters] == ['a', 'b'] + # Column j of the design belongs to parameter j + assert [parameter.sample for parameter in parameters] == pytest.approx(list(samples[i, :])) + assert all(isinstance(parameter.sample, float) for parameter in parameters) + + +def test_process_results(): + space = make_space(ot.Normal(0., 1.)) + method = MonteCarlo(n_evaluations=3, seed=42) + + outputs = method.process_results(np.array([[1., 10., 20.], [2., 11., 21.], [3., 12., 22.]]), space) + + assert len(outputs) == 3 # one output per response column + assert [output.mean for output in outputs] == pytest.approx([2., 11., 21.]) + assert all(output.method_results is None for output in outputs) # Monte Carlo has nothing beyond the samples + +def test_user_only_implements_arch_evaluate_sample(stochastic_problem): + out = stochastic_problem.evaluate(np.array([[1., 0.], [0., 0.]]), return_as_dictionary=True) + + assert out['F'].shape == (2, 1) + assert np.all(np.isfinite(out['F'])) + # E[(u-x0)^2] + x1^2 = (E[u]-x0)^2 + var + x1^2 + assert out['F'][0, 0] == pytest.approx(.05**2, abs=2e-3) + assert out['F'][1, 0] == pytest.approx(1. + .05**2, abs=2e-2) + + +def test_multiple_design_points_at_once(stochastic_problem): + x = np.random.RandomState(0).uniform(-2., 2., size=(7, 2)) + out = stochastic_problem.evaluate(x, return_as_dictionary=True) + + assert out['F'].shape == (7, 1) + assert np.all(np.isfinite(out['F'])) + assert len(np.unique(out['F'])) == 7 + + +def test_evaluation_is_deterministic_and_seeded(): + x = np.array([[.5, .5]]) + f_a = VectorizedProblem(n=50, seed=3).evaluate(x, return_as_dictionary=True)['F'] + f_b = VectorizedProblem(n=50, seed=3).evaluate(x, return_as_dictionary=True)['F'] + f_c = VectorizedProblem(n=50, seed=4).evaluate(x, return_as_dictionary=True)['F'] + + assert np.allclose(f_a, f_b) + assert not np.allclose(f_a, f_c) + + +def test_resample_draws_new_samples(): + problem = VectorizedProblem(n=50, seed=3) + x = np.array([[.5, .5]]) + f_a = problem.evaluate(x, return_as_dictionary=True)['F'].copy() + + problem.uq_method.resample() + problem.uq_method.seed = None + assert not np.allclose(f_a, problem.evaluate(x, return_as_dictionary=True)['F']) + + +def test_hierarchical_problem_with_constraint(hierarchical_problem): + x = np.array([[0, .5, .3], [1, .9, .35]]) + out = hierarchical_problem.evaluate(x, return_as_dictionary=True) + + assert out['F'].shape == (2, 1) + assert out['G'].shape == (2, 1) + assert np.all(np.isfinite(out['F'])) + assert np.all(np.isfinite(out['G'])) + + +def test_response_kinds_use_their_own_scalar(all_response_kinds_problem): + problem = all_response_kinds_problem + out = problem.evaluate(np.array([[.5]]), return_as_dictionary=True) + + # Each response kind is published in its own object array, laid out like F, G and H + assert out['F'][0, 0] == pytest.approx(Mean().scalarize(out['f_stochastic'][0, 0])) + assert out['G'][0, 0] == pytest.approx(Mean().scalarize(out['g_stochastic'][0, 0])) + assert out['H'][0, 0] == pytest.approx(Margin(k=3.).scalarize(out['h_stochastic'][0, 0])) + + # The equality constraint uses a margin, so it is above its own mean + assert out['H'][0, 0] > out['h_stochastic'][0, 0].mean + + +def test_scalar_counts_checked_per_response_kind(): + with pytest.raises(ValueError): + HierarchicalProblem(ieq_constr_scalar=[Mean(), Mean()]) + with pytest.raises(ValueError): + HierarchicalProblem(obj_scalar=[Mean(), Mean()]) + with pytest.raises(ValueError): + HierarchicalProblem(ieq_constr_scalar=['quantile']) + + HierarchicalProblem(ieq_constr_scalar=[Mean()]) # a matching count is accepted + + +def test_statistics_available_per_design_point(stochastic_problem): + out = stochastic_problem.evaluate(np.array([[1., 0.], [0., 0.]]), return_as_dictionary=True) + + assert out['f_stochastic'].shape == (2, 1) + assert all(isinstance(output, StochasticOutput) for output in out['f_stochastic'].ravel()) + assert out['g_stochastic'].shape == (2, 0) + + # The realizations really reach the model, so the response scatters + for output in out['f_stochastic'][:, 0]: + assert output.std > 0. + + +def test_reported_statistics_reproduce_the_reduced_value(): + problem = VectorizedProblem(n=200, obj_scalar=[Margin(k=2.)]) + out = problem.evaluate(np.array([[1., 0.], [.5, .5]]), return_as_dictionary=True) + + for i, output in enumerate(out['f_stochastic'][:, 0]): + assert out['F'][i, 0] == pytest.approx(Margin(k=2.).scalarize(output)) + assert out['F'][i, 0] == pytest.approx(output.mean + 2.*output.std) + + +def test_stochastic_output_statistics(): + out = _output(np.linspace(0., 10., 101)) + + assert out.mean == pytest.approx(5.) + # The kernel widens the fitted distribution slightly compared to the samples it was built from + assert out.std == pytest.approx(np.std(np.linspace(0., 10., 101)), rel=.1) + assert out.quantile(.5) == pytest.approx(5., abs=.1) + assert out.dist.computeCDF(5.) == pytest.approx(.5, abs=.05) + assert str(out) == f'(mean = {out.mean:.4g}, sigma = {out.std:.4g})' + + +@pytest.mark.parametrize('method_name', ['MonteCarlo', 'PolynomialChaos']) +def test_a_response_that_does_not_depend_on_the_parameters(method_name): + uq_method = {'MonteCarlo': MonteCarlo, 'PolynomialChaos': PolynomialChaos}[method_name](50, seed=3) + # A response that is constant over the realizations has no distribution to fit: it is reported as the value + # itself, so no scalarization can add a margin to it + problem = DeterministicResponseProblem(uq_method=uq_method, obj_scalar=[Mean(), Margin(k=2.)]) + out = problem.evaluate(np.array([[.5], [1.]]), return_as_dictionary=True) + + assert np.all(np.isfinite(out['F'])) + assert not problem.get_failed_points(out).any() + + for i, x_i in enumerate([.5, 1.]): + varying, constant = out['f_stochastic'][i, :] + assert isinstance(varying, StochasticOutput) + assert varying.std > 0. + + assert isinstance(constant, float) + assert constant == pytest.approx(2.*x_i) + assert out['F'][i, 1] == pytest.approx(2.*x_i) + + +def test_a_failed_sample_fails_the_design_point(): + # Any non-finite sample discards the design point: with a single failed evaluation the model is not + # trustworthy at that point, and SBArchOpt treats NaN as a hidden-constraint violation + problem = VectorizedProblem(n=50, fail=True) + out = problem.evaluate(np.array([[.5, .5]]), return_as_dictionary=True) + + # There is no distribution to fit, so the response comes back as NaN rather than as an output + assert not np.isfinite(out['f_stochastic'][0, 0]) + assert not np.isfinite(out['F'][0, 0]) + assert problem.get_failed_points(out)[0] + + +def test_pce_is_more_accurate_than_monte_carlo(): + # Both methods draw the same kind of LHS design, so this compares fitting an expansion against averaging + # over it. f = (u-x0)^2 + x1^2 is quadratic in u, so the expansion is essentially exact. + x = np.array([[1., 0.]]) + exact = .05**2 # E[(u-x)^2] at x = E[u] is Var[u] + + f_mc = VectorizedProblem(n=20, obj_scalar=[Mean()]).evaluate(x, return_as_dictionary=True)['F'][0, 0] + f_pce = VectorizedProblem(uq_method=PolynomialChaos(20, seed=42), obj_scalar=[Mean()]) \ + .evaluate(x, return_as_dictionary=True)['F'][0, 0] + + assert abs(f_pce - exact) < abs(f_mc - exact) + assert f_pce == pytest.approx(exact, abs=1e-4) + + +def test_pce_through_the_problem(): + # Regression: the problem called process_results without the parameter space, which PCE requires + problem = VectorizedProblem(uq_method=PolynomialChaos(50, seed=42), obj_scalar=[Mean()]) + x = np.array([[.5, .5]]) + out = problem.evaluate(x, return_as_dictionary=True) + + assert np.all(np.isfinite(out['F'])) + assert np.allclose(out['F'], problem.evaluate(x, return_as_dictionary=True)['F']) + + # Statistics come from the fitted expansion, not from the 50 expensive evaluations + assert isinstance(out['f_stochastic'][0, 0].method_results, ot.FunctionalChaosResult) + + +def test_pce_terms_and_sample_requirement(): + space = make_space(*[ot.Normal(0., 1.)]*3) + + assert PolynomialChaos(200, degree=8).n_terms(space) == 165 # degree 8 in 3 dimensions + assert PolynomialChaos(200, degree=8).n_terms(make_space(ot.Normal(0., 1.))) == 9 # ... in 1 dimension + + with pytest.raises(ValueError): + PolynomialChaos(10, degree=8).get_samples(space) # 165 terms cannot be fitted from 10 samples + + +def test_pce_provides_chaos_result_for_sensitivity_analysis(): + problem = VectorizedProblem(uq_method=PolynomialChaos(50, seed=42), obj_scalar=[Mean()]) + out = problem.evaluate(np.array([[.5, .5]]), return_as_dictionary=True) + + chaos_result = out['f_stochastic'][0, 0].method_results + assert isinstance(chaos_result, ot.FunctionalChaosResult) + + sobol = ot.FunctionalChaosSobolIndices(chaos_result) + assert sobol.getSobolIndex(0) == pytest.approx(1., abs=1e-6) # a single parameter explains everything + + +def test_pce_fails_the_design_point_when_evaluations_fail(): + # A response with failed evaluations cannot be fitted; the design point then fails as usual + problem = VectorizedProblem(uq_method=PolynomialChaos(50, seed=42), obj_scalar=[Mean()], fail=True) + out = problem.evaluate(np.array([[.5, .5]]), return_as_dictionary=True) + + assert not np.isfinite(out['f_stochastic'][0, 0]) + assert not np.isfinite(out['F'][0, 0]) \ No newline at end of file diff --git a/sb_arch_opt/uncertainty.py b/sb_arch_opt/uncertainty.py new file mode 100644 index 0000000..aaa18e0 --- /dev/null +++ b/sb_arch_opt/uncertainty.py @@ -0,0 +1,444 @@ +""" +MIT License + +Copyright: (c) 2024, Deutsches Zentrum fuer Luft- und Raumfahrt e.V. +Contact: jasper.bussemaker@dlr.de + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +""" +from __future__ import annotations + +import copy + +import numpy as np +from typing import List, Optional, Union + +try: + import openturns as ot + HAS_UNCERTAINTY = True +except ImportError: + HAS_UNCERTAINTY = False + + +def check_dependency(): + if not HAS_UNCERTAINTY: + raise ImportError( + 'Looks like SBArchOpt uncertainty package is not installed! Run: pip install sb-arch-opt[uncertainty]') + + +__all__ = ['Scalarization', 'Mean', 'Margin', 'Quantile', 'StochasticParameter', 'StochasticParameterSpace', + 'StochasticOutput', 'UQMethod', 'MonteCarlo', 'PolynomialChaos', 'EvaluationOutput', 'HAS_UNCERTAINTY', 'check_dependency'] + +EvaluationOutput = Union['StochasticOutput', float] +"""Output either distribution or numeric value.""" + +class StochasticParameter: + """ + A stochastic parameter used to evaluate architecture objectives or constraints. + + The parameter carries two distinct things. `dist` defines the parameter distribution defined by OpenTurns object. + `sample` is one realization of it. + + :param name: identifies the parameter + :param dist: OpenTurns distribution object + :param ref: Optional reference object used for matching with the related class + """ + + def __init__(self, name, dist: ot.DistributionImplementation, ref=None): + check_dependency() + self.name = name + self.dist = dist + self.ref = ref + self._sample = None + + @property + def sample(self) -> float: + """ + The parameter realization of it, written by `StochasticParameterSpace.param_realization` + """ + return self._sample + + @sample.setter + def sample(self, value: float): + self._sample = value + + +class StochasticParameterSpace: + """ + The joint distribution of all stochastic parameters of a problem, which are assumed to be independent. + The space does not store the sampled parameters as its attribute. The samples are processed by the chosen UQ method. + + :param parameters: the stochastic parameters + """ + + def __init__(self, parameters: List[StochasticParameter]): + check_dependency() + self._parameters = parameters + + @property + def n_parameters(self) -> int: + """Number of stochastic parameters""" + return len(self._parameters) + + @property + def parameters(self) -> List[StochasticParameter]: + """The stochastic parameters, in the order of the sample columns""" + return list(self._parameters) + + @property + def parameter_names(self) -> List[str]: + """List of stochastic parameter names""" + return [parameter.name for parameter in self._parameters] + + @property + def joint_dist(self) -> ot.JointDistribution: + """Joint distribution of independent parameters""" + return ot.JointDistribution([parameter.dist for parameter in self._parameters], + ot.IndependentCopula(self.n_parameters)) + + def param_realization(self, samples: np.ndarray, i_realization: int) -> List[StochasticParameter]: + """ + Return `i_realization` of `samples` (n_samples x n_parameters) as a list of stochastic parameters. + """ + realization = [] + for j, param in enumerate(self._parameters): + realized = copy.copy(param) + realized.sample = float(samples[i_realization, j]) + realization.append(realized) + return realization + + def get_random_samples(self, n_samples: int) -> np.ndarray: + """Draw n samples of the stochastic parameters; returns an n x n_parameters matrix""" + result = self.joint_dist.getSample(n_samples) + return np.array(result) + + def get_lhs_samples(self, n_samples: int) -> np.ndarray: + """Draw n samples of the stochastic parameters; returns an n x n_parameters matrix built with LHS""" + lhs = ot.LHSExperiment(self.joint_dist, n_samples) + result = lhs.generate() + return np.array(result) + + +class StochasticOutput: + """ + The object representing the stochastic output of a single objective or constraint for a given design point. + + :param dist: the output distribution for the specific objective or constraint. + """ + + def __init__(self, dist: ot.Distribution, method_results=None): + check_dependency() + self.dist = dist + self.method_results = method_results + + @staticmethod + def build_distribution(output_samples: ot.Sample) -> Union[ot.Distribution, float]: + values = np.asarray(output_samples, dtype=float).ravel() + + if values.size == 0 or not np.all(np.isfinite(values)): + return np.nan + + mean, std = float(np.mean(values)), float(np.std(values)) + if std <= 1e-12*max(abs(mean), 1.): + return mean + else: + return ot.KernelSmoothing().build(output_samples) + + @classmethod + def from_output_samples(cls, output_samples: ot.Sample, method_results=None) -> EvaluationOutput: + """ + Construct the object from evaluation output samples of a single objective or constraint. + If the variance of the provided samples is statistically insignificant, + a float value is returned to represent deterministic output. If any + sample is nan, np.nan is returned. + """ + dist = cls.build_distribution(output_samples) + if isinstance(dist, float): + return dist + return cls(dist=dist, method_results=method_results) + + @property + def mean(self) -> float: + return self.dist.getMean()[0] + + @property + def std(self) -> float: + return self.dist.getStandardDeviation()[0] + + def quantile(self, q: float) -> float: + return self.dist.computeQuantile(q)[0] + + def margin(self, k: float = 1.645, direction: int =-1) -> float: + return self.mean - np.sign(direction) * k * self.std + + def __repr__(self): + return f"(mean = {self.mean:.4g}, sigma = {self.std:.4g})" + + def __str__(self): + return repr(self) + + +class Scalarization: + """ + Parent class for all types of stochastic optimization problem. Do not instantiate this class directly, use its subclasses: + + - Mean + - Margin + - Quantile + """ + + def __init__(self): + check_dependency() + + def scalarize(self, output: StochasticOutput) -> float: + """Reduce an (n_samples x 1) sample of one response to a single value that the optimizer sees based on the optimization problem type.""" + raise NotImplementedError + + def __repr__(self): + return self.__class__.__name__ + + def __str__(self): + return repr(self) + + +class Mean(Scalarization): + """Minimize the expectation of the objective or constraint function for example min(E[F(x)])""" + + def scalarize(self, output: StochasticOutput) -> float: + return float(output.mean) + + +class Margin(Scalarization): + """ + Gaussian output distribution expected. Minimize for the objective or constraint function for a given + confidence interval. + --> min(E[F(x)] + k*sigma[F(x)]); The default k=1.645 is the one-sided 95% interval of a normal distribution. + + Requires: + + - :param k: coverage factor, i.e. the number of standard deviations added to the mean (default 1.645 covers 95% of a normal response on one side; 1, 2 and 3 give 84.1%, 97.7% and 99.9%). + - :param direction: -1 for minimizing and +1 for maximizing. + """ + + def __init__(self, k: float = 1.645, direction: int =-1): + super().__init__() + self.k = k + self.direction = direction + + def scalarize(self, output: StochasticOutput) -> float: + return float(output.margin(self.k, self.direction)) + + def __repr__(self): + return f"Margin with k = {self.k}" + + +class Quantile(Scalarization): + """ + Makes no distributional assumption. Minimize the q quantile of the objective or constraint function --> min(F_q(x)), with q=0.95 by default. + + Requires: + + :param q: quantile level in [0, 1], i.e. the probability that the response falls below the returned value (default 0.95); + unlike Margin this assumes nothing about the distribution shape, since the quantile is taken from the samples directly. + + """ + + def __init__(self, q: float = 0.95): + super().__init__() + if not 0. <= q <= 1.: + raise ValueError(f'Quantile should be between 0 and 1: {q}') + self.q = q + + def scalarize(self, output: StochasticOutput) -> float: + return float(output.quantile(self.q)) + + def __repr__(self): + return f"Quantile with q = {self.q}" + + +class UQMethod: + """ + Base class for an uncertainty propagation methods. + + It is responsible for sampling the provided 'param_space' based on the provided number of evaluations and + for performing uncertainty quantification study to return a list of StochasticOutputs or floats for a given design vector. + + Subclass by overriding `_draw_samples` (sampling method) and `process_results` (statistics computation). + + :param n_evaluations: number of expensive evaluations spent per design point + :param seed: optional seed the OpenTURNS generator before drawing, so the design is reproducible + """ + + def __init__(self, n_evaluations: int, seed: int = None): + check_dependency() + if n_evaluations is None or n_evaluations <= 0: + raise ValueError('n_evaluations must be specified: it is the number of expensive evaluations ' + 'per design point') + self.n_evaluations = n_evaluations + self.seed = seed + + self._samples: Optional[np.ndarray] = None + self._samples_space: Optional[StochasticParameterSpace] = None + + def get_samples(self, param_space: StochasticParameterSpace) -> np.ndarray: + """ + Get samples from the joint distribution provided by the param_space. The samples are drawn once and reused + for every design point (common random numbers), until `resample` is called or another space is given. + """ + if param_space is None: + raise ValueError('No parameter space to sample') + + if self._samples is None or self._samples_space is not param_space: + if self.seed is not None: + ot.RandomGenerator.SetSeed(self.seed) + self._samples = self._draw_samples(param_space) + self._samples_space = param_space + return self._samples + + def _draw_samples(self, param_space: StochasticParameterSpace) -> np.ndarray: + """Draw the design of experiments in the parameter space; override to use a different design""" + return param_space.get_random_samples(self.n_evaluations) + + def resample(self): + """Draw a new design on the next evaluation""" + self._samples = None + + def process_results(self, results: np.ndarray, param_space: StochasticParameterSpace) -> list[EvaluationOutput]: + """ + Turn the responses of one design point (an n_samples x (n_obj+n_ieq_constr+n_eq_constr) matrix) into the + list of EvaluationOutput objects. + """ + raise NotImplementedError + + def __repr__(self) -> str: + return self.__class__.__name__ + + def __str__(self): + return repr(self) + + +class MonteCarlo(UQMethod): + """ + Monte Carlo uncertainty propagation + """ + + def __init__(self, n_evaluations: int = 100, seed: int = 42): + super().__init__(n_evaluations, seed) + + def _draw_samples(self, param_space: StochasticParameterSpace) -> np.ndarray: + """Draw samples with LHS method""" + samples = param_space.get_lhs_samples(self.n_evaluations) + return samples + + def process_results(self, results: np.ndarray, param_space: StochasticParameterSpace = None) -> list[EvaluationOutput]: + results = np.asarray(results, dtype=float) + + sample = ot.Sample(results) + outputs = [StochasticOutput.from_output_samples(sample[:, i]) for i in range(results.shape[1])] + return outputs + + def __repr__(self) -> str: + return 'Monte Carlo' + + +class PolynomialChaos(UQMethod): + """ + Polynomial chaos expansion (PCE): a surrogate model is fitted for each response as a function of the uncertain parameters. + The trained model is then sampled with Monte Carlo to compute output statistics. + + The fitted `ot.FunctionalChaosResult` of each response is kept on the `StochasticOutput` (`method_result`) with + Sobol sensitivity indices available through `ot.FunctionalChaosSobolIndices`. + + :param degree: total degree of the expansion + :param n_metamodel_samples: number of cheap evaluations of the trained model used to obtain the statistics + """ + + def __init__(self, n_evaluations: int, seed: int = None, degree: int = 3, n_metamodel_samples: int = 10000): + self.degree = degree + self.n_metamodel_samples = n_metamodel_samples + self._metamodel_input: Optional[ot.Sample] = None + super().__init__(n_evaluations, seed) + + def n_terms(self, param_space: StochasticParameterSpace) -> int: + """Number of terms in the expansion, i.e. the minimum number of samples needed to fit it""" + enumerate_function = ot.LinearEnumerateFunction(param_space.n_parameters) + return int(enumerate_function.getStrataCumulatedCardinal(self.degree)) + + def _validate(self, param_space: StochasticParameterSpace): + """Check the expansion can be fitted from the budget""" + n_terms = self.n_terms(param_space) + if self.n_evaluations < n_terms: + raise ValueError(f'A degree-{self.degree} expansion in {param_space.n_parameters} ' + f'parameters has {n_terms} terms, so it needs at least that many samples to fit: ' + f'n_evaluations = {self.n_evaluations}') + + def _draw_samples(self, param_space: StochasticParameterSpace) -> np.ndarray: + """Draw samples with LHS method""" + self._validate(param_space) + samples = param_space.get_lhs_samples(self.n_evaluations) + self._metamodel_input = param_space.joint_dist.getSample(self.n_metamodel_samples) + return samples + + def _get_metamodel_input(self, param_space: StochasticParameterSpace) -> ot.Sample: + """Get samples used for evaluating the surrogate model""" + self.get_samples(param_space) + return self._metamodel_input + + def _build_algorithm(self, param_space: StochasticParameterSpace, input_sample: ot.Sample, output_sample: ot.Sample) -> ot.FunctionalChaosAlgorithm: + """Train surrogate model with the least-squares method""" + distribution = param_space.joint_dist + dimension = distribution.getDimension() + + polynomials = [ot.StandardDistributionPolynomialFactory(distribution.getMarginal(i)) + for i in range(dimension)] + enumerate_function = ot.LinearEnumerateFunction(dimension) + basis = ot.OrthogonalProductPolynomialFactory(polynomials, enumerate_function) + + adaptive_strategy = ot.FixedStrategy(basis, self.n_terms(param_space)) + projection_strategy = ot.LeastSquaresStrategy() + return ot.FunctionalChaosAlgorithm(input_sample, output_sample, distribution, + adaptive_strategy, projection_strategy) + + def process_results(self, results: np.ndarray, param_space: StochasticParameterSpace) -> list[EvaluationOutput]: + results = np.asarray(results, dtype=float) + + input_sample = ot.Sample(self.get_samples(param_space)) + metamodel_input = self._get_metamodel_input(param_space) + + outputs = [] + for i_out in range(results.shape[1]): + values = results[:, i_out] + + # A response with failed evaluations cannot be fitted: pass the raw samples through instead, so that + # the problem's nan_policy decides what happens to this design point + if not np.all(np.isfinite(values)): + samples = ot.Sample(values.reshape((-1, 1))) + chaos_result = None + + else: + algorithm = self._build_algorithm(param_space, input_sample, ot.Sample(values.reshape((-1, 1)))) + algorithm.run() + chaos_result = algorithm.getResult() + samples = chaos_result.getMetaModel()(metamodel_input) + + outputs.append(StochasticOutput.from_output_samples(samples, chaos_result)) + return outputs + + def __repr__(self) -> str: + return 'Polynomial Chaos Expansion' \ No newline at end of file diff --git a/setup.py b/setup.py index 63caeb8..ad28c3c 100644 --- a/setup.py +++ b/setup.py @@ -67,6 +67,9 @@ def _get_readme(): 'numba', 'scikit-learn', ], + 'uncertainty': [ # Stochastic (robust) problems: sb_arch_opt.uncertainty, sb_arch_opt.stochastic_problem + 'openturns', + ], # 'ota': [ # pip install -r requirements-ota.txt # 'open_turb_arch @ git+https://github.com/jbussemaker/OpenTurbofanArchitecting@pymoo_optional#egg=open_turb_arch', # ],