diff --git a/bark/benchmarks/__init__.py b/bark/benchmarks/__init__.py index bdf92ad..e9ccc8c 100644 --- a/bark/benchmarks/__init__.py +++ b/bark/benchmarks/__init__.py @@ -13,6 +13,7 @@ ) from .xgboost_mnist import XGBoostMNIST from .svr_bench import SVRBench +from .welded_beam import WeldedBeam BENCHMARK_MAP: dict[str, Type[Benchmark]] = { # unconstrained spaces @@ -38,6 +39,7 @@ # multi-fidelity # "currin": CurrinExp2D, "SVRBench": SVRBench, + "WeldedBeam": WeldedBeam, } diff --git a/bark/benchmarks/welded_beam.py b/bark/benchmarks/welded_beam.py new file mode 100644 index 0000000..2660d40 --- /dev/null +++ b/bark/benchmarks/welded_beam.py @@ -0,0 +1,164 @@ +import numpy as np +import pandas as pd +from math import sqrt +from bofire.benchmarks.api import Benchmark +from bofire.data_models.domain.api import Domain, Inputs, Outputs, Constraints +from bofire.data_models.features.api import ( + ContinuousInput, + CategoricalInput, + ContinuousOutput, +) +from bofire.data_models.objectives.api import MinimizeObjective +from bofire.data_models.constraints.api import ( + LinearInequalityConstraint, +) +from bark.bofire_utils.constraints import FunctionalInequalityConstraint +from bark.bofire_utils.domain import build_integer_input + + +class WeldedBeam(Benchmark): + """ + Welded Beam problem from https://link.springer.com/content/pdf/10.1007/s00158-018-2182-1.pdf + + Only implementing for continuous inputs. + """ + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + # Constants + self.F = 6000 + self.delta_max = 0.25 + self.L = 14 + + # Material parameters (C1, C2, sigma_d, E, G) + self.material_params = { + 0: (0.1047, 0.0481, 3e4, 3e7, 12e6), # steel + 1: (0.0489, 0.0224, 8e3, 14e6, 6e6), # cast iron + 2: (0.5235, 0.2405, 5e3, 1e7, 4e6), # aluminum + 3: (0.5584, 0.2566, 8e3, 16e6, 1e7), # brass + } + + # Define domain + self._domain = Domain( + inputs=Inputs( + features=[ + build_integer_input(key="w", bounds=[0, 1]), # binary, welding type + CategoricalInput(key="m", categories=[0, 1, 2, 3]), # material type + ContinuousInput(key="h", bounds=(0.0625, 2)), + ContinuousInput(key="l", bounds=(0.1, 10)), + ContinuousInput(key="t", bounds=(2, 20)), + ContinuousInput(key="b", bounds=(0.0625, 2)), + ] + ), + outputs=Outputs( + features=[ContinuousOutput(key="manufacturing_cost", objective=MinimizeObjective())] + ), + constraints=Constraints( + constraints=[ + FunctionalInequalityConstraint( + func=self._constraint_g1, + rhs=0.0 + ), + FunctionalInequalityConstraint( + func=self._constraint_g2, + rhs=0.0 + ), + FunctionalInequalityConstraint( + func=self._constraint_g3, + rhs=0.0 + ), + FunctionalInequalityConstraint( + func=self._constraint_g4, + rhs=0.0 + ), + FunctionalInequalityConstraint( + func=self._constraint_g5, + rhs=0.0 + ), + ] + ), + ) + + def unpack_inputs(self, x): + vars = {k: x[i] for i, k in enumerate(self.domain.inputs.get_keys())} + return vars["w"], vars["m"], vars["h"], vars["l"], vars["t"], vars["b"] + + def _get_beam_params(self, x): + """Calculate beam parameters for a single configuration""" + w, m, h, l, t, b = self.unpack_inputs(x) + + if w == 0: + # TODO: think if we want to use sqrt from corelib or sth else? + A = sqrt(2) * h * l + J = A * ((h + t)**2 / 4 + l**2 / 12) + R = 0.5 * sqrt(l**2 + (h + t)**2) + else: + A = sqrt(2) * h * (t + l) + J = ( + sqrt(2) * h * l * ((h + t)**2 / 4 + l**2 / 12) + + sqrt(2) * h * t * ((h + l)**2 / 4 + t**2 / 12) + ) + R = 0.5 * max( + sqrt(l**2 + (h + t)**2), + sqrt(t**2 + (h + l)**2) + ) + + cos_theta = l / (2 * R) + return A, J, R, cos_theta + + def _f(self, X: pd.DataFrame) -> pd.DataFrame: + """Calculate objective value (cost)""" + x = X[self.domain.inputs.get_keys()].to_numpy() + + costs = [] + for xi in x: + w, m, h, l, t, b = self.unpack_inputs(xi) + C1, C2, _, _, _ = self.material_params[int(m)] + cost = (1 + C1) * (w * t + l) * h**2 + C2 * t * b * (self.L + l) + costs.append(cost) + + return pd.DataFrame(data=np.array(costs).reshape(-1, 1), + columns=self.domain.outputs.get_keys()) + + def _constraint_g1(self, x, model_core=None): + """Shear stress constraint""" + w, m, h, l, t, b = self.unpack_inputs(x) + A, J, R, cos_theta = self._get_beam_params(x) + _, _, sigma_d, _, _ = self.material_params[int(m)] + + tau_prime = sqrt(self.F / A) + tau_double_prime = self.F * (self.L + 0.5 * l) * R / J + tau = sqrt(tau_prime**2 + tau_double_prime**2 + + 2 * tau_prime * tau_double_prime * cos_theta) + return tau - 0.577 * sigma_d + + def _constraint_g2(self, x, model_core=None): + """Normal stress constraint""" + w, m, h, l, t, b = self.unpack_inputs(x) + _, _, sigma_d, _, _ = self.material_params[int(m)] + + sigma = 6 * self.F * self.L / (t**2 * b) + return sigma - sigma_d + + def _constraint_g3(self, x, model_core=None): + """Width constraint""" + _, _, h, _, _, b = self.unpack_inputs(x) + return h - b + + def _constraint_g4(self, x, model_core=None): + """Buckling load constraint""" + w, m, h, l, t, b = self.unpack_inputs(x) + _, _, _, E, G = self.material_params[int(m)] + + P_c = (4.013 * t * b**3 * sqrt(E * G) / (6 * self.L**2) * + (1 - t / (4 * self.L) * sqrt(E / G))) + return self.F - P_c + + def _constraint_g5(self, x, model_core=None): + """Deflection constraint""" + w, m, h, l, t, b = self.unpack_inputs(x) + _, _, _, E, _ = self.material_params[int(m)] + + delta = 4 * self.F * self.L**3 / (E * t**3 * b) + return delta - self.delta_max \ No newline at end of file