Files
Atomizer/projects/hydrotech-beam/studies/01_doe_landscape/run_doe.py
Antoine 0e459028fe Fix: FEM part lookup (exclude _i.prt), hole_count unit (Constant not mm), add file logging
- solve_simulation.py: FEM finder now excludes idealized parts, falls back to loading .fem
- solve_simulation.py: hole_count written as [Constant] not [MilliMeter] in .exp
- run_doe.py: dual logging to console + results/doe_run.log
2026-02-11 14:17:43 +00:00

682 lines
24 KiB
Python

#!/usr/bin/env python3
"""Phase 1 LHS DoE — Hydrotech Beam Optimization.
Main entry point for the Design of Experiments landscape mapping study.
Generates 50 LHS sample points + 1 baseline (Trial 0), runs pre-flight
geometric checks, evaluates via NX (or stub), and manages the Optuna study.
Usage:
# Dry run with stub solver (development/testing):
python run_doe.py --backend stub --study-name hydrotech_doe_dev
# Real run with NXOpen (on Windows/dalidou):
python run_doe.py --backend nxopen --model-dir /path/to/nx/models
# Resume interrupted study:
python run_doe.py --backend stub --study-name hydrotech_doe_dev --resume
References:
OPTIMIZATION_STRATEGY.md — Full strategy document
CONTEXT.md — Confirmed expression names and values
"""
from __future__ import annotations
import argparse
import csv
import json
import logging
import os
import sys
import time
from datetime import datetime, timezone
from pathlib import Path
from typing import Any
import numpy as np
import optuna
from geometric_checks import (
DesignPoint,
FeasibilityResult,
check_feasibility,
)
from nx_interface import TrialInput, TrialResult, create_solver
from sampling import DV_DEFINITIONS, generate_lhs_samples, points_to_dicts
# ---------------------------------------------------------------------------
# Configuration
# ---------------------------------------------------------------------------
DEFAULT_STUDY_NAME = "hydrotech_beam_doe_phase1"
DEFAULT_N_SAMPLES = 50
DEFAULT_SEED = 42
DEFAULT_RESULTS_DIR = "results"
DEFAULT_DB_PATH = "results/optuna_study.db"
# Constraint limits (hard limits — OPTIMIZATION_STRATEGY.md §1.2)
DISPLACEMENT_LIMIT = 10.0 # mm
STRESS_LIMIT = 130.0 # MPa
# Infeasible placeholder values (for geometric pre-check failures)
INFEASIBLE_MASS = 99999.0 # kg
INFEASIBLE_DISPLACEMENT = 99999.0 # mm
INFEASIBLE_STRESS = 99999.0 # MPa
logger = logging.getLogger(__name__)
# ---------------------------------------------------------------------------
# Optuna constraint callback
# ---------------------------------------------------------------------------
def constraints_func(trial: optuna.trial.FrozenTrial) -> list[float]:
"""Compute constraint violations for Deb's feasibility rules.
Returns list of constraint values where:
value ≤ 0 → feasible (constraint satisfied)
value > 0 → infeasible (constraint violated)
This implements Deb's feasibility rules (Deb 2000):
1. Feasible vs feasible → lower objective wins
2. Feasible vs infeasible → feasible wins
3. Infeasible vs infeasible → lower total violation wins
References:
OPTIMIZATION_STRATEGY.md §3.2
Args:
trial: Optuna frozen trial with user attributes set.
Returns:
[displacement_violation, stress_violation]
"""
disp = trial.user_attrs.get("tip_displacement", INFEASIBLE_DISPLACEMENT)
stress = trial.user_attrs.get("max_von_mises", INFEASIBLE_STRESS)
return [
disp - DISPLACEMENT_LIMIT, # ≤ 0 means displacement ≤ 10 mm
stress - STRESS_LIMIT, # ≤ 0 means stress ≤ 130 MPa
]
# ---------------------------------------------------------------------------
# Trial evaluation
# ---------------------------------------------------------------------------
def evaluate_trial(
trial: optuna.Trial,
solver: Any,
) -> float:
"""Evaluate a single trial: geometric check → NX solve → extract.
Args:
trial: Optuna trial (with parameters already suggested/enqueued).
solver: NX solver instance (stub or real).
Returns:
Objective value (mass in kg). Returns INFEASIBLE_MASS for
geometrically infeasible or failed evaluations.
"""
# Extract design variables from trial
dv1 = trial.suggest_float(
"beam_half_core_thickness", 10.0, 40.0,
)
dv2 = trial.suggest_float(
"beam_face_thickness", 10.0, 40.0,
)
dv3 = trial.suggest_float(
"holes_diameter", 150.0, 450.0,
)
dv4 = trial.suggest_int(
"hole_count", 5, 15,
)
trial_num = trial.number
logger.info(
"Trial %d: DV1=%.2f, DV2=%.2f, DV3=%.1f, DV4=%d",
trial_num, dv1, dv2, dv3, dv4,
)
# Store DVs in user attributes for logging
trial.set_user_attr("beam_half_core_thickness", dv1)
trial.set_user_attr("beam_face_thickness", dv2)
trial.set_user_attr("holes_diameter", dv3)
trial.set_user_attr("hole_count", dv4)
# Pre-flight geometric check
point = DesignPoint(
beam_half_core_thickness=dv1,
beam_face_thickness=dv2,
holes_diameter=dv3,
hole_count=dv4,
)
geo_result: FeasibilityResult = check_feasibility(point)
trial.set_user_attr("geo_feasible", geo_result.feasible)
trial.set_user_attr("ligament", geo_result.ligament)
trial.set_user_attr("web_clearance", geo_result.web_clearance)
if not geo_result.feasible:
logger.warning(
"Trial %d: GEOMETRICALLY INFEASIBLE — %s",
trial_num, geo_result.reason,
)
trial.set_user_attr("status", "geo_infeasible")
trial.set_user_attr("geo_reason", geo_result.reason)
trial.set_user_attr("tip_displacement", INFEASIBLE_DISPLACEMENT)
trial.set_user_attr("max_von_mises", INFEASIBLE_STRESS)
trial.set_user_attr("mass", INFEASIBLE_MASS)
return INFEASIBLE_MASS
# NX evaluation
trial_input = TrialInput(
beam_half_core_thickness=dv1,
beam_face_thickness=dv2,
holes_diameter=dv3,
hole_count=dv4,
)
t_start = time.monotonic()
nx_result: TrialResult = solver.solve(trial_input)
t_elapsed = time.monotonic() - t_start
trial.set_user_attr("solve_time_s", round(t_elapsed, 2))
if not nx_result.success:
logger.error(
"Trial %d: NX SOLVE FAILED — %s (%.1fs)",
trial_num, nx_result.error_message, t_elapsed,
)
trial.set_user_attr("status", "solve_failed")
trial.set_user_attr("error_message", nx_result.error_message)
trial.set_user_attr("tip_displacement", INFEASIBLE_DISPLACEMENT)
trial.set_user_attr("max_von_mises", INFEASIBLE_STRESS)
trial.set_user_attr("mass", INFEASIBLE_MASS)
return INFEASIBLE_MASS
# Record successful results
trial.set_user_attr("status", "solved")
trial.set_user_attr("mass", nx_result.mass)
trial.set_user_attr("tip_displacement", nx_result.tip_displacement)
trial.set_user_attr("max_von_mises", nx_result.max_von_mises)
# Check constraint feasibility
disp_ok = nx_result.tip_displacement <= DISPLACEMENT_LIMIT
stress_ok = nx_result.max_von_mises <= STRESS_LIMIT
trial.set_user_attr("displacement_feasible", disp_ok)
trial.set_user_attr("stress_feasible", stress_ok)
trial.set_user_attr("fully_feasible", disp_ok and stress_ok)
logger.info(
"Trial %d: mass=%.2f kg, disp=%.2f mm%s, stress=%.1f MPa%s (%.1fs)",
trial_num,
nx_result.mass,
nx_result.tip_displacement,
"" if disp_ok else "",
nx_result.max_von_mises,
"" if stress_ok else "",
t_elapsed,
)
return nx_result.mass
# ---------------------------------------------------------------------------
# Results export
# ---------------------------------------------------------------------------
def export_csv(study: optuna.Study, output_path: str) -> None:
"""Export all trial results to CSV.
Columns: trial_number, status, beam_half_core_thickness, beam_face_thickness,
holes_diameter, hole_count, mass, tip_displacement, max_von_mises,
geo_feasible, displacement_feasible, stress_feasible, fully_feasible,
ligament, web_clearance, solve_time_s
Args:
study: Completed Optuna study.
output_path: Path to output CSV file.
"""
fieldnames = [
"trial_number",
"status",
"beam_half_core_thickness",
"beam_face_thickness",
"holes_diameter",
"hole_count",
"mass_kg",
"tip_displacement_mm",
"max_von_mises_MPa",
"geo_feasible",
"displacement_feasible",
"stress_feasible",
"fully_feasible",
"ligament_mm",
"web_clearance_mm",
"solve_time_s",
]
with open(output_path, "w", newline="") as f:
writer = csv.DictWriter(f, fieldnames=fieldnames)
writer.writeheader()
for trial in study.trials:
ua = trial.user_attrs
writer.writerow({
"trial_number": trial.number,
"status": ua.get("status", "unknown"),
"beam_half_core_thickness": ua.get("beam_half_core_thickness", ""),
"beam_face_thickness": ua.get("beam_face_thickness", ""),
"holes_diameter": ua.get("holes_diameter", ""),
"hole_count": ua.get("hole_count", ""),
"mass_kg": ua.get("mass", ""),
"tip_displacement_mm": ua.get("tip_displacement", ""),
"max_von_mises_MPa": ua.get("max_von_mises", ""),
"geo_feasible": ua.get("geo_feasible", ""),
"displacement_feasible": ua.get("displacement_feasible", ""),
"stress_feasible": ua.get("stress_feasible", ""),
"fully_feasible": ua.get("fully_feasible", ""),
"ligament_mm": ua.get("ligament", ""),
"web_clearance_mm": ua.get("web_clearance", ""),
"solve_time_s": ua.get("solve_time_s", ""),
})
logger.info("CSV results exported to %s (%d trials)", output_path, len(study.trials))
def export_summary(study: optuna.Study, output_path: str) -> None:
"""Export study summary as JSON.
Includes metadata, statistics, best feasible design, and
constraint satisfaction rates.
Args:
study: Completed Optuna study.
output_path: Path to output JSON file.
"""
trials = study.trials
n_total = len(trials)
# Count by status
n_solved = sum(1 for t in trials if t.user_attrs.get("status") == "solved")
n_geo_infeasible = sum(
1 for t in trials if t.user_attrs.get("status") == "geo_infeasible"
)
n_failed = sum(
1 for t in trials if t.user_attrs.get("status") == "solve_failed"
)
n_feasible = sum(
1 for t in trials if t.user_attrs.get("fully_feasible", False)
)
# Best feasible trial
best_feasible = None
best_mass = float("inf")
for t in trials:
if t.user_attrs.get("fully_feasible", False):
mass = t.user_attrs.get("mass", float("inf"))
if mass < best_mass:
best_mass = mass
best_feasible = t
best_info = None
if best_feasible is not None:
ua = best_feasible.user_attrs
best_info = {
"trial_number": best_feasible.number,
"mass_kg": ua.get("mass"),
"tip_displacement_mm": ua.get("tip_displacement"),
"max_von_mises_MPa": ua.get("max_von_mises"),
"design_variables": {
"beam_half_core_thickness": ua.get("beam_half_core_thickness"),
"beam_face_thickness": ua.get("beam_face_thickness"),
"holes_diameter": ua.get("holes_diameter"),
"hole_count": ua.get("hole_count"),
},
}
summary = {
"study_name": study.study_name,
"phase": "Phase 1 — LHS DoE",
"project": "Hydrotech Beam Structural Optimization",
"timestamp": datetime.now(timezone.utc).isoformat(),
"configuration": {
"n_lhs_samples": DEFAULT_N_SAMPLES,
"seed": DEFAULT_SEED,
"baseline_included": True,
"algorithm": "Latin Hypercube Sampling (scipy.stats.qmc)",
"constraint_handling": "Deb's feasibility rules",
"displacement_limit_mm": DISPLACEMENT_LIMIT,
"stress_limit_MPa": STRESS_LIMIT,
},
"design_variables": [
{
"name": dv.name,
"nx_expression": dv.nx_expression,
"lower": dv.lower,
"upper": dv.upper,
"baseline": dv.baseline,
"type": "integer" if dv.is_integer else "continuous",
}
for dv in DV_DEFINITIONS
],
"results": {
"total_trials": n_total,
"solved": n_solved,
"geo_infeasible": n_geo_infeasible,
"solve_failed": n_failed,
"fully_feasible": n_feasible,
"solve_success_rate": round(n_solved / max(n_total, 1) * 100, 1),
"feasibility_rate": round(n_feasible / max(n_solved, 1) * 100, 1),
},
"best_feasible": best_info,
"phase1_gate_check": {
"min_feasible_5": n_feasible >= 5,
"solve_success_80pct": (n_solved / max(n_total, 1)) >= 0.80,
"gate_passed": n_feasible >= 5 and (n_solved / max(n_total, 1)) >= 0.80,
},
}
with open(output_path, "w") as f:
json.dump(summary, f, indent=2)
logger.info("Summary exported to %s", output_path)
# ---------------------------------------------------------------------------
# Main study runner
# ---------------------------------------------------------------------------
def run_study(args: argparse.Namespace) -> None:
"""Execute the Phase 1 LHS DoE study.
Steps:
1. Generate LHS sample points + baseline
2. Create/load Optuna study with SQLite storage
3. Enqueue all trials (baseline first, then LHS)
4. Run optimization (all trials evaluated via objective fn)
5. Export results to CSV and JSON
Args:
args: Parsed command-line arguments.
"""
results_dir = Path(args.results_dir)
results_dir.mkdir(parents=True, exist_ok=True)
# -----------------------------------------------------------------------
# 1. Generate sample points
# -----------------------------------------------------------------------
logger.info("=" * 70)
logger.info("HYDROTECH BEAM — Phase 1 LHS DoE")
logger.info("=" * 70)
points = generate_lhs_samples(
n_samples=args.n_samples,
seed=args.seed,
include_baseline=True,
)
n_trials = len(points)
logger.info("Generated %d trial points (1 baseline + %d LHS)", n_trials, n_trials - 1)
# -----------------------------------------------------------------------
# 2. Create Optuna study
# -----------------------------------------------------------------------
db_path = results_dir / "optuna_study.db"
storage = f"sqlite:///{db_path}"
if args.clean and db_path.exists():
logger.info("--clean flag: deleting existing DB at %s", db_path)
db_path.unlink()
if args.resume:
logger.info("Resuming existing study: %s", args.study_name)
study = optuna.load_study(
study_name=args.study_name,
storage=storage,
)
logger.info("Loaded study with %d existing trials", len(study.trials))
else:
study = optuna.create_study(
study_name=args.study_name,
storage=storage,
direction="minimize", # minimize mass
load_if_exists=True, # safe re-run: reuse if exists
sampler=optuna.samplers.TPESampler(seed=args.seed),
)
logger.info("Created new study: %s (storage: %s)", args.study_name, db_path)
# Enqueue all trial points
trial_dicts = points_to_dicts(points)
for i, td in enumerate(trial_dicts):
study.enqueue_trial(td)
if i == 0:
logger.info("Enqueued Trial 0 (baseline): %s", td)
logger.info("Enqueued %d trials (1 baseline + %d LHS)", n_trials, n_trials - 1)
# -----------------------------------------------------------------------
# 3. Create solver
# -----------------------------------------------------------------------
solver = create_solver(
backend=args.backend,
model_dir=args.model_dir,
)
# -----------------------------------------------------------------------
# 4. Run all trials
# -----------------------------------------------------------------------
logger.info("-" * 70)
logger.info("Starting trial evaluations...")
logger.info("-" * 70)
t_study_start = time.monotonic()
# Suppress Optuna's verbose logging during trials
optuna.logging.set_verbosity(optuna.logging.WARNING)
study.optimize(
lambda trial: evaluate_trial(trial, solver),
n_trials=n_trials,
callbacks=[_progress_callback],
)
t_study_elapsed = time.monotonic() - t_study_start
# -----------------------------------------------------------------------
# 5. Export results
# -----------------------------------------------------------------------
logger.info("-" * 70)
logger.info("Study complete. Exporting results...")
logger.info("-" * 70)
csv_path = str(results_dir / "doe_results.csv")
json_path = str(results_dir / "doe_summary.json")
export_csv(study, csv_path)
export_summary(study, json_path)
# -----------------------------------------------------------------------
# 6. Print summary
# -----------------------------------------------------------------------
_print_summary(study, t_study_elapsed)
# Cleanup
solver.close()
def _progress_callback(study: optuna.Study, trial: optuna.trial.FrozenTrial) -> None:
"""Log progress after each trial."""
n_complete = len(study.trials)
status = trial.user_attrs.get("status", "unknown")
mass = trial.user_attrs.get("mass", "N/A")
feasible = trial.user_attrs.get("fully_feasible", False)
logger.info(
" [%d/%d] Trial %d: status=%s, mass=%s, feasible=%s",
n_complete,
DEFAULT_N_SAMPLES + 1,
trial.number,
status,
f"{mass:.2f} kg" if isinstance(mass, (int, float)) else mass,
feasible,
)
def _print_summary(study: optuna.Study, elapsed: float) -> None:
"""Print a human-readable summary to stdout."""
trials = study.trials
n_total = len(trials)
n_solved = sum(1 for t in trials if t.user_attrs.get("status") == "solved")
n_geo_inf = sum(1 for t in trials if t.user_attrs.get("status") == "geo_infeasible")
n_failed = sum(1 for t in trials if t.user_attrs.get("status") == "solve_failed")
n_feasible = sum(1 for t in trials if t.user_attrs.get("fully_feasible", False))
print("\n" + "=" * 70)
print("PHASE 1 DoE — RESULTS SUMMARY")
print("=" * 70)
print(f" Total trials: {n_total}")
print(f" Solved: {n_solved}")
print(f" Geometrically infeasible: {n_geo_inf}")
print(f" Solve failures: {n_failed}")
print(f" Fully feasible: {n_feasible}")
print(f" Solve success rate: {n_solved/max(n_total,1)*100:.1f}%")
print(f" Feasibility rate: {n_feasible/max(n_solved,1)*100:.1f}%")
print(f" Total time: {elapsed:.1f}s ({elapsed/60:.1f} min)")
# Phase 1 → Phase 2 gate check
print("\n GATE CHECK (Phase 1 → Phase 2):")
gate_feasible = n_feasible >= 5
gate_solve = (n_solved / max(n_total, 1)) >= 0.80
print(f" ≥5 feasible points: {'✓ PASS' if gate_feasible else '✗ FAIL'} ({n_feasible})")
print(f" ≥80% solve success: {'✓ PASS' if gate_solve else '✗ FAIL'} ({n_solved/max(n_total,1)*100:.0f}%)")
print(f" GATE: {'✓ PASSED' if gate_feasible and gate_solve else '✗ BLOCKED'}")
# Best feasible
best_mass = float("inf")
best_trial = None
for t in trials:
if t.user_attrs.get("fully_feasible", False):
m = t.user_attrs.get("mass", float("inf"))
if m < best_mass:
best_mass = m
best_trial = t
if best_trial is not None:
ua = best_trial.user_attrs
print(f"\n BEST FEASIBLE DESIGN (Trial {best_trial.number}):")
print(f" Mass: {ua['mass']:.2f} kg")
print(f" Displacement: {ua['tip_displacement']:.2f} mm (limit: {DISPLACEMENT_LIMIT} mm)")
print(f" VM Stress: {ua['max_von_mises']:.1f} MPa (limit: {STRESS_LIMIT} MPa)")
print(f" Core thickness: {ua['beam_half_core_thickness']:.3f} mm")
print(f" Face thickness: {ua['beam_face_thickness']:.3f} mm")
print(f" Hole diameter: {ua['holes_diameter']:.1f} mm")
print(f" Hole count: {ua['hole_count']}")
else:
print("\n ⚠️ NO FEASIBLE DESIGN FOUND — see OPTIMIZATION_STRATEGY.md §7.1")
print("=" * 70)
# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------
def parse_args() -> argparse.Namespace:
"""Parse command-line arguments."""
parser = argparse.ArgumentParser(
description="Hydrotech Beam — Phase 1 LHS DoE Study",
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog=(
"Examples:\n"
" # Development (stub solver):\n"
" python run_doe.py --backend stub\n\n"
" # Real NX evaluation:\n"
" python run_doe.py --backend nxopen --model-dir /path/to/models\n\n"
" # Resume interrupted study:\n"
" python run_doe.py --backend stub --resume\n"
),
)
parser.add_argument(
"--backend",
choices=["stub", "nxopen"],
default="stub",
help="NX solver backend: 'stub' for testing, 'nxopen' for real (default: stub)",
)
parser.add_argument(
"--model-dir",
default="",
help="Path to NX model directory (required for --backend nxopen)",
)
parser.add_argument(
"--study-name",
default=DEFAULT_STUDY_NAME,
help=f"Optuna study name (default: {DEFAULT_STUDY_NAME})",
)
parser.add_argument(
"--n-samples",
type=int,
default=DEFAULT_N_SAMPLES,
help=f"Number of LHS sample points (default: {DEFAULT_N_SAMPLES})",
)
parser.add_argument(
"--seed",
type=int,
default=DEFAULT_SEED,
help=f"Random seed for LHS (default: {DEFAULT_SEED})",
)
parser.add_argument(
"--results-dir",
default=DEFAULT_RESULTS_DIR,
help=f"Output directory for results (default: {DEFAULT_RESULTS_DIR})",
)
parser.add_argument(
"--resume",
action="store_true",
help="Resume an existing study instead of creating a new one",
)
parser.add_argument(
"--clean",
action="store_true",
help="Delete existing results DB before starting (fresh run)",
)
parser.add_argument(
"--verbose", "-v",
action="store_true",
help="Enable verbose (DEBUG) logging",
)
return parser.parse_args()
def main() -> None:
"""Entry point."""
args = parse_args()
# Configure logging — console + file
log_level = logging.DEBUG if args.verbose else logging.INFO
log_format = "%(asctime)s [%(levelname)-7s] %(name)s: %(message)s"
log_datefmt = "%Y-%m-%d %H:%M:%S"
# Ensure results dir exists for log file
results_dir = Path(args.results_dir)
results_dir.mkdir(parents=True, exist_ok=True)
log_file = results_dir / "doe_run.log"
logging.basicConfig(
level=log_level,
format=log_format,
datefmt=log_datefmt,
handlers=[
logging.StreamHandler(), # console
logging.FileHandler(log_file, mode="a", encoding="utf-8"), # file
],
)
logger.info("Log file: %s", log_file.resolve())
# Run
try:
run_study(args)
except KeyboardInterrupt:
logger.warning("Study interrupted by user. Progress saved to Optuna DB.")
logger.info("Resume with: python run_doe.py --resume --study-name %s", args.study_name)
sys.exit(1)
except Exception:
logger.exception("Study failed with unexpected error")
sys.exit(2)
if __name__ == "__main__":
main()