benchmarking theta results for noiseless likelihood variance

now going to re-evaluate theta controllers for possible improvement

Author: dantzert

baseline_controllers/gamma/optimize_parameters.py

@@ -172,10 +172,78 @@ def run_swmm(constant_flows, efd_parameters=None,verbose=False):
 
     return {"data_log": env.data_log,"peak_filling_degrees":peak_filling_degrees}
 
+
+# --- Shared evaluator for objective/constraint with caching ---
+def evaluate_cf_point(params_array):
+    """
+    Compute objective and constraint once, cache into Sim_cf.computed_return_values,
+    and return a dict with 'objective' and 'constraint'.
+    """
+    key = tuple(np.array(params_array).flatten())
+    if key in Sim_cf.computed_return_values:
+        return Sim_cf.computed_return_values[key]
+
+    data = run_swmm(np.array(params_array).flatten(), None, verbose=False)
+
+    # Objective (flow exceedance)
+    flow_cost = 0.0
+    for node_id, series in data['data_log']['flow'].items():
+        if '5' in node_id or '9' in node_id:
+            continue
+        flow_exceed = [x - 3.0 for x in series]
+        flow_exceed = [x if x > 0 else 0 for x in flow_exceed]
+        flow_cost += sum(flow_exceed)
+    objective_cost = float(flow_cost)
+
+    # Constraint = flood + drainage
+    flood_cost = 0.0
+    for node_id, series in data['data_log']['flooding'].items():
+        if '5' in node_id or '9' in node_id:
+            continue
+        flood_cost += sum(series)
+    if 0.0 < flood_cost < 1.0:
+        flood_cost = 1.0
+    elif flood_cost <= 0.0:
+        flood_cost = float(max(data['peak_filling_degrees']))
+    else:
+        flood_cost = float(flood_cost)
+
+    drainage_cost = 0.0
+    for node_id, series in data['data_log']['depthN'].items():
+        if '5' in node_id or '9' in node_id:
+            continue
+        final_depth = series[-1]
+        if final_depth > 0.10:
+            drainage_cost += 10.0 * final_depth
+
+    constraint_cost = float(flood_cost + drainage_cost)
+
+    Sim_cf.computed_return_values[key] = {
+        "objective": objective_cost,
+        "constraint": constraint_cost,
+    }
+    return Sim_cf.computed_return_values[key]
+
 class Sim_cf:
     threshold = 0.99 # on the constraint function to define the feasible (safe) region
     computed_return_values = dict()
 
+    @staticmethod
+    def objective(input_data):
+        return_values = []
+        for sample in input_data:
+            res = evaluate_cf_point(sample.numpy())
+            return_values.append(res["objective"])
+        return np.array(return_values).reshape(-1, 1)
+
+    @staticmethod
+    def constraint(input_data):
+        return_values = []
+        for sample in input_data:
+            res = evaluate_cf_point(sample.numpy())
+            return_values.append(res["constraint"])
+        return np.array(return_values).reshape(-1, 1)
+    '''
     @staticmethod
     def objective(input_data):
         return_values = []
@@ -263,6 +331,7 @@ def constraint(input_data):
                 return_values.append(flood_cost)
         return_values = np.array(return_values).reshape(-1,1)
         return return_values
+    '''
 
 OBJECTIVE = "OBJECTIVE"
 CONSTRAINT = "CONSTRAINT"
@@ -406,7 +475,7 @@ def observer_efd(query_points):
         }
 
 def create_bo_model(data):
-        gpr = build_gpr(data, search_space)
+        gpr = build_gpr(data, search_space,likelihood_variance = 1e-7)
         return GaussianProcessRegression(gpr)
 
 if evaluating == "constant-flow" and mode == "optimize":
@@ -670,10 +739,17 @@ def combined_objective(params):
         function_call_count += 1
 
         params_array = np.array(params).flatten()
-        
+        res = evaluate_cf_point(params_array)
+        # If the constraint is violated, return a large penalty
+        if res['constraint'] > Sim_cf.threshold:
+            return 1e6*(res['constraint'] - Sim_cf.threshold) # if all initial samples have the same value some methods will error out.
+        return float(res['objective'])
+
+        '''
         # Use existing simulation infrastructure to evaluate
         data = run_swmm(params_array, None, verbose=False)
         return sum(data['data_log']['performance_measure'])
+        '''
 
     # For tracking performance across methods
     results = {
@@ -684,9 +760,9 @@ def combined_objective(params):
     }
 
     # Generate common initial points for all methods
-    num_initial_points = 50  # Number of initial points
-    num_steps = 200  
-    initial_seed = 7  # Use fixed seed for reproducibility
+    num_initial_points = 5#25  # Number of initial points
+    num_steps = 5#25  
+    initial_seed = 42  # Use fixed seed for reproducibility
     tf.random.set_seed(initial_seed)
     np.random.seed(initial_seed)
 
@@ -724,46 +800,46 @@ def combined_objective(params):
     datasets = opt_result.try_get_final_datasets()
     obj_data = datasets[OBJECTIVE]
     con_data = datasets[CONSTRAINT]
+    models = opt_result.try_get_final_models()  # capture models (not used for 9D plotting)
     bouc_fcalls = len(obj_data.query_points)
-        # Process BOUC results and calculate pystorms costs
+
+    # Process BOUC results and compute combined costs for each point (no re-simulation)
+    bouc_combined_costs = []
     for i in range(len(obj_data.query_points)):
         point = obj_data.query_points[i].numpy()
-        
-        # Calculate pystorms cost for this point
-        data = run_swmm(point, None, verbose=False)
-        pystorms_cost = sum(data['data_log']['performance_measure'])
-        bouc_pystorms_costs.append(pystorms_cost)
-        
-        # Record original BOUC objective value
+
+        combined_cost = combined_objective(point)
+        bouc_combined_costs.append(combined_cost)
+
+        # Record original BOUC objective/constraint
         obj_value = float(obj_data.observations[i][0])
         con_value = float(con_data.observations[i][0])
-        
+
         # Track best feasible point according to BOUC's objective
         if con_value <= Sim_cf.threshold and obj_value < best_feasible_obj:
             best_feasible_obj = obj_value
             best_feasible_point = point
-            best_feasible_idx = i  # Store the index here for later reference
-            
+            best_feasible_idx = i
+
         bouc_costs.append(obj_value)
         bouc_fcalls_arr.append(i + 1)
 
     # Store results for comparison
     results["BOUC"]["costs"] = bouc_costs
-    results["BOUC"]["pystorms_costs"] = bouc_pystorms_costs
+    results["BOUC"]["combined_costs"] = bouc_combined_costs
     results["BOUC"]["fcalls"] = bouc_fcalls_arr
-    
-    # Use the already calculated pystorms cost for the best point
+
+    # Use the combined cost for the best point
     if best_feasible_point is not None:
         results["BOUC"]["best_params"] = best_feasible_point
-        results["BOUC"]["best_cost"] = bouc_pystorms_costs[best_feasible_idx]
+        results["BOUC"]["best_cost"] = bouc_combined_costs[best_feasible_idx]
     else:
-        # If no feasible point, use the best objective point
-        idx = np.argmin(obj_data.observations)
+        idx = int(np.argmin(bouc_combined_costs))
         results["BOUC"]["best_params"] = obj_data.query_points[idx].numpy()
-        results["BOUC"]["best_cost"] = bouc_pystorms_costs[idx]
+        results["BOUC"]["best_cost"] = bouc_combined_costs[idx]
 
     print("BOUC function calls:", bouc_fcalls_arr)
-    print("BOUC pystorms costs:", bouc_pystorms_costs)
+    print("BOUC combined costs:", bouc_combined_costs)
 
     # 2. Vanilla Bayesian Optimization using Trieste
     print("\nRunning vanilla BO optimization...")
@@ -791,7 +867,7 @@ def observer_vanilla_bo(query_points):
         return Dataset(query_points, Sim_vanilla_bo.objective(query_points))
 
     def vanilla_bo_create_model(data):
-        gpr = build_gpr(data, search_space)
+        gpr = build_gpr(data, search_space, likelihood_variance = 1e-7)
         return GaussianProcessRegression(gpr)
 
     initial_data_bo = observer_vanilla_bo(initial_points)
@@ -940,8 +1016,8 @@ def de_callback(x, convergence):
 
     # Compare results
     print("\nOptimization Results Comparison (using consistent performance measure):")
-    print("Starting costs (from pystorms performance measure):")
-    print(f"  BOUC: {bouc_pystorms_costs[0]:.4f}, BO: {bo_costs[0]:.4f}, DA: {da_costs[0]:.4f}, DE: {de_costs[0]:.4f}")
+    print("Starting costs (combined objective):")
+    print(f"  BOUC: {bouc_combined_costs[0]:.4f}, BO: {bo_costs[0]:.4f}, DA: {da_costs[0]:.4f}, DE: {de_costs[0]:.4f}")
     print("-" * 80)
     print(f"{'Method':<10} | {'Best Cost':<15} | {'Function Calls':<15} | {'Time (s)':<15}")
     print("-" * 80)
@@ -1021,7 +1097,7 @@ def rolling_min(costs, feasible=None):
         return best
 
     # For BOUC, use pystorms cost directly
-    results["BOUC"]["best_cost_so_far"] = rolling_min(results["BOUC"]["pystorms_costs"])
+    results["BOUC"]["best_cost_so_far"] = rolling_min(results["BOUC"]["combined_costs"])
 
     # For other methods, all are assumed feasible
     for method in ["BO", "DA", "DE"]:

baseline_controllers/gamma/v2/optimal_constant_flows_summary.csv

@@ -1,5 +1,5 @@
 method,best_cost,optimization_time_sec,num_function_calls,param_1,param_2,param_3,param_4,param_5,param_6,param_7,param_8,param_9
-BOUC,54604.68450075984,3228.6177713871,250,5.436529783332854,5.2185728521246135,4.7906485474226,4.763275302682986,4.502900992326628,4.938501015614886,4.48096928667281,3.2572279236826436,2.5745598571846338
-BO,54604.68450075984,2494.4057030677795,250,5.436529783332854,5.2185728521246135,4.7906485474226,4.763275302682986,4.502900992326628,4.938501015614886,4.48096928667281,3.2572279236826436,2.5745598571846338
-DA,6684.060085752899,3958.8013985157013,450,3.1841683238744736,3.71714245271869,3.4829832687973976,4.632918010640424,4.454999765657703,2.623481869697571,4.822746708989143,2.5657217912375927,2.5587263798806816
-DE,11068.532715935731,1497.0523014068604,200,3.7389973237442766,4.0230993228590535,3.80864666266616,4.2640294186878425,4.775061835437011,3.8027835724900605,4.63957584796589,4.013554629503549,5.99345575483958
+BOUC,99290.67299849597,34477.51378130913,50,3.6077411389789056,5.3807122663141955,3.442493276477518,5.7232907117410665,5.337139884533713,4.190866514161556,2.729240969602435,4.5890491048784785,2.7649846856003766
+BO,98377.50329707502,535.5013947486877,50,3.931620357112107,5.10573759247095,3.6278875241137896,4.84199169640375,4.3314561496114035,2.5151946820151045,4.246216602957377,2.6639378348946607,5.570273166335934
+DA,68970.11886400828,897.9846568107605,75,3.262722931802273,4.638648509979248,3.1899866349995136,5.378272503614426,4.7226347625255585,3.721466556424275,5.884222760796547,2.5597217231988907,5.8646092265844345
+DE,87445.70084261333,954.190426826477,100,3.2419585373433937,3.9608576461699503,3.395929034032598,5.83247917716061,3.1556220277169844,2.5558112183094126,2.8870635771375257,4.553924883217957,5.1252473575821496