Ex 1: Automatic Topological Generation & Surrogate Optimization (2D / 4-DOF)

This notebook represents the complete methodology for imprecise tolerance analysis on the 2D baseline assembly.

Building on the manual vector loops constructed previously, this notebook introduces the streamlined workflow required to scale the methodology to high-dimensional systems:

  1. Automatic Graph-Based Assembly Modeling: Automatically constructing the SystemOfConstraintsAssemblyModel from a high-level topological dictionary.

  2. Multi-Layer Perceptron (MLP) Surrogate: Replacing expensive direct geometric optimizations with a PyTorch neural network.

  3. Automated Tracking & Constraint Injection: Using otaf.optimization.OptimizationTracker and the model’s native constraint generator to manage the credal space boundaries.

  4. Min-Max Bounds Resolution: Finding a feasible warm-start (\(x_0\)) and optimizing across threshold levels (\(\alpha_k\)) to compute the P-Box.

[1]:

import numpy as np
import scipy
import openturns as ot
import matplotlib.pyplot as plt
import torch
from scipy.optimize import minimize, NonlinearConstraint, Bounds

import otaf
from otaf.example_models import model1
from gldpy import GLD

ot.Log.Show(ot.Log.NONE)
np.set_printoptions(suppress=True)

1. High-Level Topological Dictionary & Automated Loops

Instead of manually defining deviation and gap matrices, we define the assembly as a topological graph. The otaf.AssemblyDataProcessor interprets the interactions and expands the first-order loops.

(Note: For the surrogate and optimization phases, we will dynamically pull the parameters directly from ``otaf.example_models.model1``, but this dictionary illustrates how the physics of ``model1`` are constructed).

[2]:

# Nominal dimensions and base frames
X1, X2, X3 = 99.8, 100.0, 10.0
R0 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
x_, y_, z_ = R0[0], R0[1], R0[2]

# Part 1 and 2 points
P1A0, P1A1, P1A2 = np.array((0, X3/2, 0.0)), np.array((0, X3, 0.0)), np.array((0, 0, 0.0))
P1B0, P1B1, P1B2 = np.array((X1, X3/2, 0.0)), np.array((X1, X3, 0.0)), np.array((X1, 0, 0.0))
P1C0, P1C1, P1C2 = np.array((X1/2, 0, 0.0)), np.array((0, 0, 0.0)), np.array((X1, 0, 0.0))
P2A0, P2A1, P2A2 = np.array((0, X3/2, 0.0)), np.array((0, X3, 0.0)), np.array((0, 0, 0.0))
P2B0, P2B1, P2B2 = np.array((X2, X3/2, 0.0)), np.array((X2, X3, 0.0)), np.array((X2, 0, 0.0))
P2C0, P2C1, P2C2 = np.array((X2/2, 0, 0.0)), np.array((0, 0, 0.0)), np.array((X2, 0, 0.0))

2. High-Level Topological Dictionary

Instead of manually defining deviation and gap matrices, we define the assembly as a topological graph. The otaf.AssemblyDataProcessor interprets the interactions and automatically expands the first-order loops.

[3]:

system_data = {
    "PARTS" : {
        '1' : {
            "a" : {"FRAME": np.array([-1*x_, -1*y_, z_]), "POINTS": {'A0': P1A0, 'A1': P1A1, 'A2': P1A2}, "TYPE": "plane", "INTERACTIONS": ['P2a'], "CONSTRAINTS_D": ["PERFECT"], "CONSTRAINTS_G": ["FLOATING"]},
            "b" : {"FRAME": R0, "POINTS": {'B0': P1B0, 'B1': P1B1, 'B2': P1B2}, "TYPE": "plane", "INTERACTIONS": ['P2b'], "CONSTRAINTS_D": ["NONE"], "CONSTRAINTS_G": ["FLOATING"]},
            "c" : {"FRAME": np.array([-y_, x_, z_]), "POINTS": {'C0': P1C0, 'C1': P1C1, 'C2': P1C2}, "TYPE": "plane", "INTERACTIONS": ['P2c'], "CONSTRAINTS_D": ["PERFECT"], "CONSTRAINTS_G": ["SLIDING"]},
        },
        '2' : {
            "a" : {"FRAME": R0, "POINTS": {'A0': P2A0, 'A1': P2A1, 'A2': P2A2}, "TYPE": "plane", "INTERACTIONS": ['P1a'], "CONSTRAINTS_D": ["PERFECT"], "CONSTRAINTS_G": ["FLOATING"]},
            "b" : {"FRAME": np.array([-1*x_, -1*y_, z_]), "POINTS": {'B0': P2B0, 'B1': P2B1, 'B2': P2B2}, "TYPE": "plane", "INTERACTIONS": ['P1b'], "CONSTRAINTS_D": ["NONE"], "CONSTRAINTS_G": ["FLOATING"]},
            "c" : {"FRAME": np.array([y_, -1*x_, z_]), "POINTS": {'C0': P2C0, 'C1': P2C1, 'C2': P2C2}, "TYPE": "plane", "INTERACTIONS": ['P1c'], "CONSTRAINTS_D": ["PERFECT"], "CONSTRAINTS_G": ["SLIDING"]},
        }
    },
    "LOOPS": {
        "COMPATIBILITY": {
            "L0": "P1cC0 -> P2cC0 -> P2aA0 -> P1aA0",
            "L1": "P1cC0 -> P2cC0 -> P2bB0 -> P1bB0",
        },
    },
    "GLOBAL_CONSTRAINTS": "2D_NZ",
}

[4]:

# Process the topology
SDA = otaf.AssemblyDataProcessor(system_data)
SDA.generate_expanded_loops()

# Extract constraints automatically and embed the slack variable (s)
CLH = otaf.CompatibilityLoopHandling(SDA)
ILH = otaf.InterfaceLoopHandling(SDA, CLH, circle_resolution=20)
SOCAM = otaf.SystemOfConstraintsAssemblyModel(
    CLH.get_compatibility_expressions(),
    ILH.get_interface_expressions()
)
SOCAM.embedOptimizationVariable()

print(f"Total DOFs: {len(SOCAM.deviation_symbols)}")
print(f"Variables: {SOCAM.deviation_symbols}")

Total DOFs: 4
Variables: [u_d_4, gamma_d_4, u_d_5, gamma_d_5]

2. Surrogate Training (MLP)

To bypass the expensive direct optimizations, we fit a PyTorch surrogate. We draw the distribution parameters directly from the model1 definition, taking advantage of the multiply_composed_distribution_standard_with_constants utility to push the sampling into the failure space for robust training.

[5]:

# Fetch distributions and parameters natively from model1
jointDist, symbols, max_std_vect, mu_vect = model1.getDistributionParams(tol=0.31, capa=1.0)
dim = model1.dim

# Artificial dispersion multiplier to ensure limit-state representation
mult = 1.35
expanded_dist = otaf.distribution.multiply_composed_distribution_standard_with_constants(jointDist, [mult]*dim)

# Generate Sample
sample_size = 100000
np.random.seed(420)
TRAIN_SAMPLE = np.array(expanded_dist.getSample(sample_size), dtype="float32")

# Compute exact gap optimizations for training
TRAIN_RESULTS = otaf.uncertainty.compute_gap_optimizations_on_sample_batch(
    SOCAM, TRAIN_SAMPLE, bounds=None, n_cpu=-2, progress_bar=True, batch_size=500, dtype="float32"
)

failure_ratio = np.where(TRAIN_RESULTS[:, -1] < 0, 1, 0).sum() / sample_size
print(f"Ratio of failed simulations in sample: {failure_ratio:.5f}")

Ratio of failed simulations in sample: 0.14665

[ ]:

# Neural Network Initialization and Training
neural_model = otaf.surrogate.NeuralRegressorNetwork(
    input_dim=dim, output_dim=1, X=TRAIN_SAMPLE, y=TRAIN_RESULTS[:, -1],
    clamping=True, finish_critertion_epoch=5, loss_finish=1e-6,
    metric_finish=0.9998, max_epochs=500, batch_size=30000,
    compile_model=False, train_size=0.6, display_progress_disable=False
)

neural_model.model = torch.nn.Sequential(
    *otaf.surrogate.get_custom_mlp_layers([dim, 100, 70, 30, 1], activation_class=torch.nn.GELU)
)
neural_model.optimizer = torch.optim.AdamW(neural_model.parameters(), lr=0.003, weight_decay=1e-4)
otaf.surrogate.initialize_model_weights(neural_model)
neural_model.scheduler = torch.optim.lr_scheduler.ExponentialLR(neural_model.optimizer, 1.0001)
neural_model.loss_fn = torch.nn.MSELoss()

neural_model.train_model()

3. Defining the Evaluation Closure and Optimization Logic

We define a closure get_model_evaluator to map the normalized inputs (\(x \in [0,1]\)) directly to the un-normalized model evaluation.

We also wrap the optimization execution using the otaf.optimization.OptimizationTracker to reliably store all configuration evaluations, successful GLD parameters, and failure probabilities for later visualization.

[ ]:

def get_model_evaluator(sample, mu_vect, surrogate_model):
    """Closure mapping scaled input x to unstandardized surrogate space."""
    def evaluate(x):
        x_transformed = (sample - mu_vect) * x + mu_vect
        prediction = surrogate_model.evaluate_model_non_standard_space(x_transformed)
        return np.squeeze(prediction.numpy())
    return evaluate

@otaf.optimization.scaling(scale_factor=1.0)
def optimization_function(
        x, failure_slack=0.0, gld=None, model=None,
        experiment_key=None, tracker=None, logprob=False, minimize=True):

    multiplier = 1 if minimize else -1
    slack = model(x)
    gld_params = gld.fit_LMM(slack, disp_fit=False, disp_optimizer=False)

    fp_slack = np.where(slack < failure_slack, 1, 0).mean()
    if np.any(np.isnan(gld_params)):
        fp_out = fp_slack
    else:
        fp_out = gld.CDF_num(failure_slack, gld_params, xtol=1e-6)

    tracker.update_objective_data(
        exp_key=experiment_key, x=x, fp_gld=fp_out,
        fp_slack=fp_slack, gld_params=gld_params, failure_slack=failure_slack
    )

    if logprob:
        return multiplier * np.log(1e-16 + fp_out)
    return multiplier * fp_out

4. Constraint Factory and Feasible Warm Start (\(x_0\))

Instead of manually constructing tolerance equations, we pull the encapsulated credal set constraint function natively from model1. We use SLSQP to push a generic uniform array to the closest valid point on the boundary, serving as a reliable warm start for COBYQA.

[ ]:

def optimize_scaling_vector(constraint_fn, n_vars, x_warm=None):
    """Finds the independent scaling factors x closest to 1.0 satisfying the credal constraints."""
    if x_warm is None:
        x_warm = np.full(n_vars, 0.5)

    result = scipy.optimize.minimize(
        fun=lambda x: np.sum((x - 1.0) ** 2),
        jac=lambda x: 2.0 * (x - 1.0),
        x0=x_warm, method="SLSQP",
        bounds=Bounds(lb=1e-7, ub=1.0, keep_feasible=True),
        constraints={"type": "ineq", "fun": lambda x: -np.array(constraint_fn(x))},
        options={"ftol": 1e-10, "maxiter": 500, "disp": False},
    )
    return result.x

def create_directional_constraint(tr, expKey, is_minimization):
    lb = 0.0 if is_minimization else -np.inf
    ub = np.inf if is_minimization else 0.0
    return NonlinearConstraint(
        fun=model1.getScaledCredalSetConstraintsFunction(max_std_vect, tr, expKey),
        lb=lb, ub=ub, keep_feasible=True
    )

5. Min-Max Global Optimization Loop

With the tracker, closure, and constraints configured, we iterate over the alpha levels (\(\alpha_k\)).

[ ]:

def pf_min_max_optimizer(failure_slack, tracker, experiment_key, model_eval_fn, x0, dim):
    gld = GLD('VSL')
    normalized_bounds = Bounds(1e-9, 1.0, keep_feasible=True)

    print(f"\nOptimizing bounds for Slack = {failure_slack}...")

    res_maxi = scipy.optimize.minimize(
        optimization_function, x0,
        args=(failure_slack, gld, model_eval_fn, experiment_key, tracker, True, False),
        method="COBYQA", bounds=normalized_bounds,
        constraints=create_directional_constraint(tracker, experiment_key, False),
        options={"f_target": -np.inf, "maxiter": 1000, "feasibility_tol": 1e-5, "disp": False}
    )

    res_mini = scipy.optimize.minimize(
        optimization_function, x0,
        args=(failure_slack, gld, model_eval_fn, experiment_key, tracker, True, True),
        method="COBYQA", bounds=normalized_bounds,
        constraints=create_directional_constraint(tracker, experiment_key, True),
        options={"f_target": -np.inf, "maxiter": 1000, "feasibility_tol": 1e-5, "disp": False}
    )
    return res_mini, res_maxi

[ ]:

# Execution
sample_gld = np.array(jointDist.getSample(100000))
evaluator = get_model_evaluator(sample_gld, mu_vect, neural_model)
tracker = otaf.optimization.OptimizationTracker(bounds=Bounds(1e-7, 1.0), constraint_tolerance=1e-5, precision_decimals=8)

# Obtain feasible starting point
print("Finding feasible start (x0)...")
base_constraint = model1.getScaledCredalSetConstraintsFunction(max_std_vect)
x0 = optimize_scaling_vector(base_constraint, dim)

[ ]:

# Sweep alpha levels
alpha_levels = [0.0, 0.025, 0.05, 0.075]
for slack in alpha_levels:
    pf_min_max_optimizer(
        failure_slack=slack,
        tracker=tracker,
        experiment_key=f"model1_slack_{slack}",
        model_eval_fn=evaluator,
        x0=x0,
        dim=dim
    )

[ ]:

# Save the centralized tracker dictionary to CSV for analysis
df = tracker.to_dataframe()
df.to_csv("OptimizationResults_model1_4_dof.csv")
print("\nOptimization Complete. Data saved to OptimizationResults_model1_4_dof.csv")

[ ]: