Constrained Optimization Example

This script demonstrates the optimization of a constrained problem using Multi-Objective Bayesian Optimization.

It leverages the QPOTS framework to optimize the Disc Brake problem, a multi-objective problem with constraints.

Overview

  • Uses BoTorch and PyTorch for Gaussian Process (GP) modeling.

  • Implements Pareto Optimal Thompson Sampling (QPOTS) for optimization.

  • Evaluates hypervolume (HV) to measure performance.

  • Saves results (candidates, hypervolume, and timing) for post-analysis.

Script Details

constrained_optimization.py
  1"""
  2This file demonstrates an optimization of a constrained problem
  3"""
  4
  5import warnings
  6import time
  7import os
  8import numpy as np
  9
 10warnings.filterwarnings('ignore')
 11
 12from qpots.acquisition import Acquisition
 13from qpots.config import DEFAULT_DEVICE, DEFAULT_DTYPE
 14from qpots.config import DEFAULT_DEVICE, DEFAULT_DTYPE
 15from qpots.model_object import ModelObject
 16from qpots.utils.utils import expected_hypervolume
 17from qpots.utils.utils import posterior_mean_fill
 18from qpots.function import Function
 19
 20import torch
 21from botorch.utils.transforms import unnormalize, normalize
 22
 23test_function_tag="weldedbeam"
 24if test_function_tag =="discbrake":
 25    dim_pass=4
 26    nobj_pass=2
 27    ntrain_pass=10*dim_pass
 28    ncons_pass=4
 29    ref_pass=[5.8, 4.0]
 30    var_pass=[-0.5,-0.5]
 31elif test_function_tag =="weldedbeam":
 32    dim_pass=4
 33    nobj_pass=2
 34    ntrain_pass=10*dim_pass
 35    ncons_pass=4
 36    ref_pass=[40, 10]
 37    var_pass=[-0.5,-0.5]
 38
 39device = DEFAULT_DEVICE
 40args = dict(
 41        {
 42            "ntrain": ntrain_pass,
 43            "iters": 50,
 44            "reps": 20,
 45            "q": 2,
 46            "wd": "..",
 47            "ref_point": -1*torch.tensor(ref_pass, device=device, dtype=DEFAULT_DTYPE),
 48            "dim": dim_pass,
 49            "nobj": nobj_pass,
 50            "ncons": ncons_pass,
 51            "nystrom": 0,
 52            "nychoice": "pareto",
 53            "ngen": 10,
 54            "mt": 1,
 55            "partial_info": 0,
 56            "variance_threshold": None, #torch.tensor([9.790e-5,9.790e-5,9.790e-5,9.790e-5,9.790e-5,9.790e-5])
 57        }
 58    )
 59
 60tf = Function(test_function_tag, dim=args["dim"], nobj=args["nobj"])
 61f = tf.evaluate
 62bounds = tf.get_bounds()
 63cons = tf.get_cons()
 64
 65os.makedirs(args["wd"], exist_ok=True)
 66torch.manual_seed(1023) #1023
 67
 68train_x = torch.rand([args["ntrain"], args["dim"]], device=device, dtype=DEFAULT_DTYPE)
 69train_y = f(unnormalize(train_x, bounds))
 70train_y = torch.column_stack([train_y, cons(unnormalize(train_x, bounds))]) # Stack constraints on top of objectives
 71full_y=train_y
 72
 73print(train_y.shape, train_x.shape) # This should be n_train x (nobj + ncons) tensor
 74
 75gps = ModelObject(train_x, train_y, bounds, args["nobj"], args["ncons"], args["ntrain"], device=device)
 76if args["mt"]==1:
 77    gps.fit_multitask_gp()
 78else:
 79    gps.fit_gp()
 80
 81acq = Acquisition(tf, gps, cons=cons, device=device, q=args["q"])
 82
 83hvs, times = [], []
 84for i in range(args["iters"]):
 85
 86    t1 = time.time() # tracking time
 87    if args["partial_info"]==1:
 88        newx,new_task_id = acq.qpots(bounds, i, **args)
 89    else:
 90        newx = acq.qpots(bounds, i, **args)
 91        
 92    t2 = time.time()
 93    times.append(t2 - t1)
 94    
 95
 96    if args["partial_info"]==1:
 97        #Getting full cons and y
 98        full_newy = f(unnormalize(newx.reshape(-1, args["dim"]), bounds))
 99        full_newconsy = cons(unnormalize(newx.reshape(-1, args["dim"]), bounds))
100
101        #Attaching constraints 
102        full_newy = torch.column_stack([full_newy.reshape(newx.shape[0], args["nobj"]),
103                                full_newconsy.reshape(newx.shape[0], args["ncons"])])
104        
105        newy = torch.full_like(full_newy, float('nan'))
106
107        for j in range(newx.shape[0]):
108            cols = new_task_id[j]
109            valid_mask = ~torch.isnan(cols)           
110            cols = cols[valid_mask].long()             
111            newy[j, cols] = full_newy[j, cols]
112
113        #newy = torch.column_stack([newy.reshape(newx.shape[0], args["nobj"]), newconsy.reshape(newx.shape[0], args["ncons"])])
114        
115        
116        print("Partial_info newy:\n",newy)
117    else:
118        newy = f(unnormalize(newx.reshape(-1, args["dim"]), bounds))
119        newconsy = cons(unnormalize(newx.reshape(-1, args["dim"]), bounds))
120        newy = torch.column_stack([newy.reshape(args["q"], args["nobj"]),
121                                newconsy.reshape(args["q"], args["ncons"])])
122    
123    hv, _ = expected_hypervolume(gps, ref_point=args['ref_point'])
124    hvs.append(hv)
125
126    print(f"Iteration: {i}, New candidate: {newx}, Time: {t2 - t1}, HV: {hv}")   
127
128    train_x = torch.row_stack([train_x, newx.view(-1, args["dim"])])
129    train_y = torch.row_stack([train_y, newy])
130    
131    if args["partial_info"]==1:
132        full_y=torch.row_stack([full_y, full_newy])
133
134    gps = ModelObject(train_x, train_y, bounds, args["nobj"], args["ncons"], args["ntrain"], device=device)
135    if args["mt"]==1:
136        if args["variance_threshold"] is None:
137            if args["partial_info"] == 1:
138                tag="rand"
139            else:
140                tag="joint"
141        else:
142            tag="var_thresh"
143        gps.fit_multitask_gp()
144        np.save(f"{args['wd']}/cons_"+tag+"_train_x.npy", train_x)
145        np.save(f"{args['wd']}/cons_"+tag+"_train_y.npy", train_y)
146        hvs_tensor = torch.stack(hvs)
147        np.save(f"{args['wd']}/cons_"+tag+"_hv.npy", hvs_tensor.detach().cpu().numpy())
148        np.save(f"{args['wd']}/cons_"+tag+"_times.npy", times)
149        if args["partial_info"]==1:
150            np.save(f"{args['wd']}/cons_"+tag+"_full_y.npy", full_y) #Full y is without the NaNs, using for pareto sorting later
151    else:
152        gps.fit_gp()
153        np.save(f"{args['wd']}/cons_Model_list_train_x.npy", train_x)
154        np.save(f"{args['wd']}/cons_Model_list_train_y.npy", train_y)
155        np.save(f"{args['wd']}/cons_Model_list_hv.npy", hvs)
156        np.save(f"{args['wd']}/cons_Model_list_times.npy", times)
157
158
159if args["partial_info"]==1:
160    train_y_filled=posterior_mean_fill(gps)
161    np.save(f"{args['wd']}/cons_"+tag+"_train_y_filled.npy", train_y_filled.detach().cpu().numpy())
162

Example Output

Iteration: 0, New candidate: tensor([...]), Time: 0.23s, HV: 107.278
Iteration: 1, New candidate: tensor([...]), Time: 0.19s, HV: 123.899
...
Iteration: 199, New candidate: tensor([...]), Time: 0.30s, HV: 151.326

Usage

Run the script using:

python examples/constrained_example.py

Ensure that dependencies such as BoTorch, PyTorch, and PyMoo are installed.