fix(thermoelastic)!: enforce volume fraction as a constraint, not an objective

docs/problems/thermoelastic2d.md

@@ -19,6 +19,7 @@ The design space is then defined by a 2D array representing density values (para
 The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material.
 Total compliance is defined as a linear combination of thermal compliance and structural compliance.
 The weight configuration parameter defines the relative importance of thermal compliance and structural compliance in the total compliance calculation, where a weight of 1 corresponds to a purely structural problem and a weight of 0 corresponds to a purely thermal problem.
+The target volume fraction `volfrac` is enforced as a design constraint rather than reported as an objective, following the Beams2D convention.
 
 ## Conditions
 
@@ -32,7 +33,7 @@ The optimization process itself operates by calculating the sensitivities of the
 The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold.
 
 ## Dataset
-The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_2d_v1).
+The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_2d_v2).
 This dataset contains a set of 1000 optimized thermoelastic designs in a 64x64 domain, where each design is optimized for a unique set of conditions.
 Each datapoint's conditions are randomly generated by arbitrarily placing: a single loaded element along the bottom boundary, two fixed elements (fixed in both the x and y direction) along the left and top boundary, and heatsink elements along the right boundary.
 Furthermore, values for the volume fraction constraint are randomly selected in the range $[0.2, 0.5]$.
@@ -45,11 +46,10 @@ Relevant datapoint fields include:
 - `force_elements_x`: Encodes a binary NxN matrix specifying elements that have a structural load in the x-direction.
 - `force_elements_y`: Encodes a binary NxN matrix specifying elements that have a structural load in the y-direction.
 - `heatsink_elements`: Encodes a binary NxN matrix specifying elements that have a heat sink.
-- `volume_fraction_error`: The volume fraction error with respect to the target volume fraction
 - `structural_compliance`: The structural compliance of the optimized design
 - `thermal_compliance`: The thermal compliance of the optimized design
 - `nelx`: The number of elements in the x-direction
 - `nely`: The number of elements in the y-direction
-- `volume_fraction_target`: The volume fraction target of the optimized design
+- `volfrac`: The volume fraction target of the optimized design
 - `rmin`: The filter size used in the optimization routine
 - `weight`: The domain weighting used in the optimization routine

docs/problems/thermoelastic3d.md

@@ -20,6 +20,7 @@ The design space is then defined by a 3D array representing density values (para
 The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material.
 Total compliance is defined as a linear combination of thermal compliance and structural compliance.
 The weight configuration parameter defines the relative importance of thermal compliance and structural compliance in the total compliance calculation, where a weight of 1 corresponds to a purely structural problem and a weight of 0 corresponds to a purely thermal problem.
+The target volume fraction `volfrac` is enforced as a design constraint rather than reported as an objective, following the Beams2D convention.
 
 ## Conditions
 
@@ -33,7 +34,7 @@ The optimization process itself operates by calculating the sensitivities of the
 The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold.
 
 ## Dataset
-The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_3d_v0).
+The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_3d_v1).
 This dataset contains a set of 100 optimized thermoelastic designs in a 16x16x16 domain, where each design is optimized for a unique set of conditions.
 Each datapoint's conditions are randomly generated by arbitrarily placing: a single loaded element along the bottom boundary, two fixed elements (fixed in the x, y, and z direction) along the left and top boundary, and heatsink elements along the right boundary.
 Furthermore, values for the volume fraction constraint are randomly selected in the range $[0.2, 0.5]$.
@@ -45,7 +46,6 @@ Relevant datapoint fields include:
 - `force_elements_y`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the y-direction.
 - `force_elements_z`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the z-direction.
 - `heatsink_elements`: Encodes a binary NxNxN matrix specifying elements that have a heat sink.
-- `volume_fraction`: The volume fraction value of the optimized design
 - `structural_compliance`: The structural compliance of the optimized design
 - `thermal_compliance`: The thermal compliance of the optimized design
 - `nelx`: The number of elements in the x-direction

engibench/problems/beams3d/model/fem_model.py

@@ -235,10 +235,8 @@ def run(  # noqa: PLR0915
             f0val = f0valm
 
             if self.eval_only:
-                vf_error = abs(np.mean(x) - volfrac)
                 return {
                     "structural_compliance": float(f0valm),
-                    "volume_fraction": vf_error,
                 }
 
             obj_values = np.array([f0valm])
@@ -323,12 +321,10 @@ def run(  # noqa: PLR0915
             opti_steps[-1].obj_values_update = np.zeros_like(opti_steps[-1].obj_values)
 
         print("3D structural optimization finished.")
-        vf_error = abs(np.mean(x) - volfrac)
 
         return {
             "design": x,
             "bcs": bcs,
             "structural_compliance": float(f0valm),
-            "volume_fraction": vf_error,
             "opti_steps": opti_steps,
         }

engibench/problems/beams3d/v0.py

@@ -14,6 +14,7 @@
 
 from engibench.constraint import bounded
 from engibench.constraint import constraint
+from engibench.constraint import Criticality
 from engibench.constraint import greater_than
 from engibench.constraint import IMPL
 from engibench.constraint import THEORY
@@ -52,6 +53,16 @@ def _fixed_elements(nelx: int, nely: int, nelz: int) -> npt.NDArray[np.int64]:
     return out
 
 
+@constraint(categories=THEORY, criticality=Criticality.Warning)
+def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None:
+    """Constraint for volume fraction of the design."""
+    actual_volfrac = design.mean()
+    tolerance = 0.01
+    assert abs(actual_volfrac - volfrac) <= tolerance, (
+        f"Volume fraction of the design {actual_volfrac:.4f} does not match target {volfrac:.4f} specified in the conditions. While the optimizer might fix it, this is likely to affect objective values as the initial design is not feasible given the constraints."
+    )
+
+
 class Beams3D(Problem[npt.NDArray]):
     """3D structural topology optimization problem."""
 
@@ -75,6 +86,7 @@ class Conditions:
         """Fractional y-position of the vertical load on the top face."""
 
     conditions = Conditions()
+    design_constraints = (volume_fraction_bound,)
     design_space = spaces.Box(low=0.0, high=1.0, shape=(NELY, NELX, NELZ), dtype=np.float32)
     dataset_id = f"IDEALLab/beams_3d_{NELX}_v0"
     container_id = None

engibench/problems/thermoelastic2d/model/fea_model.py

@@ -161,7 +161,7 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
         Args:
             bcs (dict[str, any]): A dictionary containing boundary conditions and problem parameters.
                 Expected keys include:
-                    - 'volume_fraction_target' (float): Target volume fraction.
+                    - 'volfrac' (float): Target volume fraction.
                     - 'fixed_elements' (np.ndarray): NxN binary array encoding the location of fixed elements.
                     - 'force_elements_x' (np.ndarray): NxN binary array encoding the location of loaded elements in the x direction.
                     - 'force_elements_y' (np.ndarray): NxN binary array encoding the location of loaded elements in the y direction.
@@ -174,10 +174,10 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
             Dict[str, Any]: A dictionary containing the optimization results. The dictionary includes:
                 - 'design' (np.ndarray): Final design layout.
                 - 'bcs' (Dict[str, Any]): The input boundary conditions.
-                - 'sc' (float): Structural cost component.
-                - 'tc' (float): Thermal cost component.
-                - 'vf' (float): Volume fraction error.
-            If self.eval_only is True, returns a dictionary with keys 'sc', 'tc', and 'vf' only.
+                - 'structural_compliance' (float): Structural compliance.
+                - 'thermal_compliance' (float): Thermal compliance.
+                - 'opti_steps' (list[OptiStep]): The optimization history.
+            If self.eval_only is True, returns a dictionary with keys 'structural_compliance' and 'thermal_compliance' only.
         """
         # WEIGHTING
         w1 = bcs.get("weight", 0.5)
@@ -188,7 +188,7 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
         nelx = fe_h - 1
         nely = fe_w - 1
 
-        volfrac = bcs["volume_fraction_target"]
+        volfrac = bcs["volfrac"]
         n = nely * nelx  # Total number of elements
 
         # OptiSteps records
@@ -326,16 +326,13 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
             f0val = (f0valm * w1) + (f0valt * w2)
 
             if self.eval_only is True:
-                vf_error = np.abs(np.mean(x) - volfrac)
                 return {
                     "structural_compliance": f0valm,
                     "thermal_compliance": f0valt,
-                    "volume_fraction_error": vf_error,
                 }
 
             # OptiStep Information
-            vf_error = np.abs(np.mean(x) - volfrac)
-            obj_values = np.array([f0valm, f0valt, vf_error])
+            obj_values = np.array([f0valm, f0valt])
             x_curr = x.copy()  # Design variables before the gradient update (nely, nelx)
 
             df0dx = df0dx_mat.reshape(nely * nelx, 1)
@@ -414,14 +411,12 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
                 extra_iter = True
 
         print("Optimization finished...")
-        vf_error = np.abs(np.mean(x) - volfrac)
 
         return {
             "design": x,
             "bcs": bcs,
             "structural_compliance": f0valm,
             "thermal_compliance": f0valt,
-            "volume_fraction_error": vf_error,
             "opti_steps": opti_steps,
         }
 
@@ -440,7 +435,7 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
         "fixed_elements": [lci[21], lci[32], lci[43]],
         "force_elements_y": [bri[31]],
         "heatsink_elements": [lci[31], lci[32], lci[33]],
-        "volume_fraction_target": 0.2,
+        "volfrac": 0.2,
         "rmin": 1.1,
         "weight": 1.0,  # 1.0 for pure structural, 0.0 for pure thermal
     }

engibench/problems/thermoelastic2d/v0.py

@@ -14,6 +14,7 @@
 
 from engibench.constraint import bounded
 from engibench.constraint import constraint
+from engibench.constraint import Criticality
 from engibench.constraint import IMPL
 from engibench.constraint import THEORY
 from engibench.core import ObjectiveDirection
@@ -33,6 +34,16 @@
 HEATSINK_ELEMENTS = indices_to_binary_matrix([LCI[31], LCI[32], LCI[33]], NELX + 1, NELY + 1)
 
 
+@constraint(categories=THEORY, criticality=Criticality.Warning)
+def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None:
+    """Constraint for volume fraction of the design."""
+    actual_volfrac = design.mean()
+    tolerance = 0.01
+    assert abs(actual_volfrac - volfrac) <= tolerance, (
+        f"Volume fraction of the design {actual_volfrac:.4f} does not match target {volfrac:.4f} specified in the conditions. While the optimizer might fix it, this is likely to affect objective values as the initial design is not feasible given the constraints."
+    )
+
+
 class ThermoElastic2D(Problem[npt.NDArray]):
     r"""Truss 2D integer optimization problem.
 
@@ -43,7 +54,6 @@ class ThermoElastic2D(Problem[npt.NDArray]):
     objectives: tuple[tuple[str, ObjectiveDirection], ...] = (
         ("structural_compliance", ObjectiveDirection.MINIMIZE),
         ("thermal_compliance", ObjectiveDirection.MINIMIZE),
-        ("volume_fraction_error", ObjectiveDirection.MINIMIZE),
     )
 
     @dataclass
@@ -66,7 +76,7 @@ class Conditions:
             default_factory=lambda: HEATSINK_ELEMENTS
         )
         """Binary NxN matrix specifying elements that have a heat sink"""
-        volume_fraction_target: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.3
+        volfrac: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.3
         """Target volume fraction for the volume fraction constraint"""
         rmin: Annotated[
             float, bounded(lower=1.0).category(THEORY), bounded(lower=0.0, upper=3.0).warning().category(IMPL)
@@ -76,8 +86,9 @@ class Conditions:
         """Control which objective is optimized for. 1.0 is pure structural optimization, while 0.0 is pure thermal optimization"""
 
     conditions = Conditions()
+    design_constraints = (volume_fraction_bound,)
     design_space = spaces.Box(low=0.0, high=1.0, shape=(NELX, NELY), dtype=np.float32)
-    dataset_id = "IDEALLab/thermoelastic_2d_v1"
+    dataset_id = "IDEALLab/thermoelastic_2d_v2"
     container_id = None
 
     @dataclass
@@ -130,9 +141,7 @@ def simulate_verbose(self, design: npt.NDArray, config: dict[str, Any] | None =
                     boundary_dict[key] = value
 
         results = FeaModel(plot=False, eval_only=True).run(boundary_dict, x_init=design)
-        return SimulationResult(
-            np.array([results["structural_compliance"], results["thermal_compliance"], results["volume_fraction_error"]])
-        )
+        return SimulationResult(np.array([results["structural_compliance"], results["thermal_compliance"]]))
 
     def optimize(
         self, starting_point: npt.NDArray, config: dict[str, Any] | None = None

engibench/problems/thermoelastic3d/model/fem_model.py

@@ -198,10 +198,10 @@ def run(self, bcs: dict[str, Any], x_init: np.ndarray | None = None) -> dict[str
             Dict[str, Any]: A dictionary containing the optimization results. The dictionary includes:
                 - 'design' (np.ndarray): Final design layout.
                 - 'bcs' (Dict[str, Any]): The input boundary conditions.
-                - 'sc' (float): Structural cost component.
-                - 'tc' (float): Thermal cost component.
-                - 'vf' (float): Volume fraction error.
-            If self.eval_only is True, returns a dictionary with keys 'sc', 'tc', and 'vf' only.
+                - 'structural_compliance' (float): Structural compliance.
+                - 'thermal_compliance' (float): Thermal compliance.
+                - 'opti_steps' (list[OptiStep]): The optimization history.
+            If self.eval_only is True, returns a dictionary with keys 'structural_compliance' and 'thermal_compliance' only.
         """
         # Weighting
         w1 = bcs.get("weight", 0.5)  # structural
@@ -368,14 +368,11 @@ def node_id(ix: int, iy: int, iz: int) -> int:
             f0val = (f0valm * w1) + (f0valt * w2)
 
             if self.eval_only:
-                vf_error = abs(np.mean(x) - volfrac)
                 return {
                     "structural_compliance": float(f0valm),
                     "thermal_compliance": float(f0valt),
-                    "volume_fraction": vf_error,
                 }
-            vf_error = np.abs(np.mean(x) - volfrac)
-            obj_values = np.array([f0valm, f0valt, vf_error])
+            obj_values = np.array([f0valm, f0valt])
             x_curr = x.copy()
 
             xval = x.reshape(n, 1)
@@ -473,13 +470,11 @@ def node_id(ix: int, iy: int, iz: int) -> int:
                 extra_iter = True
 
         print("3D optimization finished.")
-        vf_error = abs(np.mean(x) - volfrac)
 
         return {
             "design": x,
             "bcs": bcs,
             "structural_compliance": float(f0valm),
             "thermal_compliance": float(f0valt),
-            "volume_fraction": vf_error,
             "opti_steps": opti_steps,
         }

engibench/problems/thermoelastic3d/v0.py

@@ -12,6 +12,7 @@
 
 from engibench.constraint import bounded
 from engibench.constraint import constraint
+from engibench.constraint import Criticality
 from engibench.constraint import IMPL
 from engibench.constraint import THEORY
 from engibench.core import ObjectiveDirection
@@ -48,6 +49,16 @@ def _default_heatsink_elements(nelx: int, nely: int, nelz: int) -> npt.NDArray[n
     return out
 
 
+@constraint(categories=THEORY, criticality=Criticality.Warning)
+def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None:
+    """Constraint for volume fraction of the design."""
+    actual_volfrac = design.mean()
+    tolerance = 0.01
+    assert abs(actual_volfrac - volfrac) <= tolerance, (
+        f"Volume fraction of the design {actual_volfrac:.4f} does not match target {volfrac:.4f} specified in the conditions. While the optimizer might fix it, this is likely to affect objective values as the initial design is not feasible given the constraints."
+    )
+
+
 class ThermoElastic3D(Problem[npt.NDArray]):
     """Truss 3D integer optimization problem.
 
@@ -58,7 +69,6 @@ class ThermoElastic3D(Problem[npt.NDArray]):
     objectives: tuple[tuple[str, ObjectiveDirection], ...] = (
         ("structural_compliance", ObjectiveDirection.MINIMIZE),
         ("thermal_compliance", ObjectiveDirection.MINIMIZE),
-        ("volume_fraction", ObjectiveDirection.MINIMIZE),
     )
 
     @dataclass
@@ -87,8 +97,9 @@ class Conditions:
         weight: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.5
         """Control which objective is optimized for. 1.0 is pure structural optimization, while 0.0 is pure thermal optimization"""
 
+    design_constraints = (volume_fraction_bound,)
     design_space = spaces.Box(low=0.0, high=1.0, shape=(NELX, NELY, NELZ), dtype=np.float32)
-    dataset_id = "IDEALLab/thermoelastic_3d_v0"
+    dataset_id = "IDEALLab/thermoelastic_3d_v1"
     container_id = None
 
     @dataclass
@@ -206,9 +217,7 @@ def simulate_verbose(self, design: npt.NDArray, config: dict[str, Any] | None =
                     boundary_dict[key] = value
 
         results = FeaModel3D(eval_only=True).run(boundary_dict, x_init=design)
-        return SimulationResult(
-            np.array([results["structural_compliance"], results["thermal_compliance"], results["volume_fraction"]])
-        )
+        return SimulationResult(np.array([results["structural_compliance"], results["thermal_compliance"]]))
 
     def optimize(
         self, starting_point: npt.NDArray, config: dict[str, Any] | None = None

tests/reference/simulate/problems.thermoelastic2d.v0.ThermoElastic2D.json

@@ -1,8 +1,7 @@
 {
   "performance": [
     104399698.05000441,
-    7603.173244221379,
-    0.05000413414090871
+    7603.173244221379
   ],
   "rtol": 5e-05
 }

tests/reference/simulate/problems.thermoelastic3d.v0.ThermoElastic3D.json

@@ -1,8 +1,7 @@
 {
   "performance": [
     2330.7670793294246,
-    4050.8264802757444,
-    0.1499882920390519
+    4050.8264802757444
   ],
   "rtol": 1e-07
 }