From 1731af766639525cfe9b0c0f4076a8fd176edf22 Mon Sep 17 00:00:00 2001 From: nychiang Date: Tue, 8 Apr 2025 12:26:35 -0700 Subject: [PATCH 1/5] add support to use ipopt to solve constrainted problem --- .github/workflows/test_hiopbbpy.yml | 2 +- pyproject.toml | 2 +- setup.py | 4 +- .../hiopbbpy/{BODriver.py => BODriverCI.py} | 0 src/Drivers/hiopbbpy/BODriverEX.py | 93 ++++++++++++++ src/hiopbbpy/opt/__init__.py | 2 + src/hiopbbpy/opt/acquisition.py | 51 +++++++- src/hiopbbpy/opt/boalgorithm.py | 113 +++++++++++++----- src/hiopbbpy/opt/optproblem.py | 53 ++++++++ src/hiopbbpy/surrogate_modeling/gp.py | 34 +++++- src/hiopbbpy/surrogate_modeling/krg.py | 21 +++- 11 files changed, 326 insertions(+), 49 deletions(-) rename src/Drivers/hiopbbpy/{BODriver.py => BODriverCI.py} (100%) create mode 100644 src/Drivers/hiopbbpy/BODriverEX.py create mode 100644 src/hiopbbpy/opt/optproblem.py diff --git a/.github/workflows/test_hiopbbpy.yml b/.github/workflows/test_hiopbbpy.yml index e4d9289fe..2d5117bf3 100644 --- a/.github/workflows/test_hiopbbpy.yml +++ b/.github/workflows/test_hiopbbpy.yml @@ -26,7 +26,7 @@ jobs: - name: Run Tests run: | - python src/Drivers/hiopbbpy/BODriver.py 10 + python src/Drivers/hiopbbpy/BODriverCI.py 10 continue-on-error: false diff --git a/pyproject.toml b/pyproject.toml index 424783deb..955ab2bf7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,3 @@ [build-system] -requires = ["setuptools", "wheel", "numpy", "smt"] +requires = ["setuptools", "wheel", "numpy", "smt", "cyipopt"] build-backend = "setuptools.build_meta" \ No newline at end of file diff --git a/setup.py b/setup.py index 9fbaed53d..903b554f4 100644 --- a/setup.py +++ b/setup.py @@ -12,14 +12,14 @@ metadata = dict( name="hiopbbpy", - version="0.0.3", + version="0.0.4", description="HiOp black box optimization (hiopbbpy)", author="Tucker hartland et al.", author_email="hartland1@llnl.gov", license="BSD-3", packages=find_packages(where="src"), package_dir={"": "src"}, - install_requires=["smt"], + install_requires=["smt","cyipopt"], python_requires=">=3.9", zip_safe=False, url="https://github.com/LLNL/hiop", diff --git a/src/Drivers/hiopbbpy/BODriver.py b/src/Drivers/hiopbbpy/BODriverCI.py similarity index 100% rename from src/Drivers/hiopbbpy/BODriver.py rename to src/Drivers/hiopbbpy/BODriverCI.py diff --git a/src/Drivers/hiopbbpy/BODriverEX.py b/src/Drivers/hiopbbpy/BODriverEX.py new file mode 100644 index 000000000..c9841008c --- /dev/null +++ b/src/Drivers/hiopbbpy/BODriverEX.py @@ -0,0 +1,93 @@ +""" + Code description: + for a 2D example LpNormProblem + 1) randomly sample training points + 2) define a Kriging-based Gaussian-process (smt backend) + trained on said data + 3) determine the minimizer via BOAlgorithm + +Authors: Tucker Hartland + Nai-Yuan Chiang +""" + +import sys +import os +import numpy as np +import matplotlib.pyplot as plt +import warnings +warnings.filterwarnings("ignore") +from LpNormProblem import LpNormProblem +from hiopbbpy.surrogate_modeling import smtKRG +from hiopbbpy.opt import BOAlgorithm +from hiopbbpy.problems import BraninProblem + + +# Get user input for the number of repetitions from command-line arguments +if len(sys.argv) != 2 or int(sys.argv[1]) < 0: + num_repeat = 1 +else: + num_repeat = int(sys.argv[1]) + +### parameters +n_samples = 5 # number of the initial samples to train GP +theta = 1.e-2 # hyperparameter for GP kernel +nx = 2 # dimension of the problem +xlimits = np.array([[-5, 5], [-5, 5]]) # bounds on optimization variable + +prob_type_l = ["LpNorm"] # ["LpNorm", "Branin"] +acq_type_l = ["LCB"] # ["LCB", "EI"] + +def con_eq(x): + return x[0] + x[1] - 4 + +def con_jac_eq(x): + return np.array([1.0, 1.0]) + +def con_ineq(x): + return x[0] - x[1] + +def con_jac_ineq(x): + return np.array([1.0, -1.0]) + +# only 'trust-constr' method supports vector-valued constraints +user_constraint = [{'type': 'ineq', 'fun': con_ineq, 'jac': con_jac_ineq}, + {'type': 'eq', 'fun': con_eq, 'jac': con_jac_eq}] + +retval = 0 +for prob_type in prob_type_l: + print() + if prob_type == "LpNorm": + problem = LpNormProblem(nx, xlimits) + else: + problem = BraninProblem() + + for acq_type in acq_type_l: + print("Problem name: ", problem.name) + print("Acquisition type: ", acq_type) + + ### initial training set + x_train = problem.sample(n_samples) + y_train = problem.evaluate(x_train) + + ### Define the GP surrogate model + gp_model = smtKRG(theta, xlimits, nx) + gp_model.train(x_train, y_train) + + options = { + 'acquisition_type': acq_type, + 'bo_maxiter': 10, + 'opt_solver': 'IPOPT', #"SLSQP" "IPOPT" + 'solver_options': { + 'maxiter': 200 + } + } + + # Instantiate and run Bayesian Optimization + bo = BOAlgorithm(gp_model, x_train, y_train, options = options, user_constraints = user_constraint) #EI or LCB + bo.optimize(problem) + +sys.exit(retval) + + + + diff --git a/src/hiopbbpy/opt/__init__.py b/src/hiopbbpy/opt/__init__.py index 163941272..88cfcdb82 100644 --- a/src/hiopbbpy/opt/__init__.py +++ b/src/hiopbbpy/opt/__init__.py @@ -1,5 +1,6 @@ from .boalgorithm import (BOAlgorithmBase, BOAlgorithm) from .acquisition import (acquisition, LCBacquisition, EIacquisition) +from .optproblem import (IpoptProbFromScipy) __all__ = [ "BOAlgorithmBase" @@ -7,4 +8,5 @@ "acquisition" "LCBacquisition" "EIacquisition" + "IpoptProbFromScipy" ] diff --git a/src/hiopbbpy/opt/acquisition.py b/src/hiopbbpy/opt/acquisition.py index 9bc602b63..ed8ef834f 100644 --- a/src/hiopbbpy/opt/acquisition.py +++ b/src/hiopbbpy/opt/acquisition.py @@ -14,28 +14,41 @@ class acquisition(object): def __init__(self, gpsurrogate): assert isinstance(gpsurrogate, GaussianProcess) # add something here self.gpsurrogate = gpsurrogate + self.has_gradient = False # Abstract method to evaluate the acquisition function at x. def evaluate(self, x: np.ndarray) -> np.ndarray: raise NotImplementedError("Child class of acquisition should implement method evaluate") + # Abstract method to evaluate the gradient of acquisition function at x. + def eval_g(self, x: np.ndarray) -> np.ndarray: + raise NotImplementedError("Child class of acquisition should implement method evaluate") # A subclass of acquisition, implementing the Lower Confidence Bound (LCB) acquisition function. class LCBacquisition(acquisition): def __init__(self, gpsurrogate, beta=3.0): super().__init__(gpsurrogate) self.beta = beta + self.has_gradient = True # Method to evaluate the acquisition function at x. def evaluate(self, x : np.ndarray) -> np.ndarray: mu = self.gpsurrogate.mean(x) - sig = self.gpsurrogate.variance(x) - return mu - self.beta * np.sqrt(sig) + sig2 = self.gpsurrogate.variance(x) + return mu - self.beta * np.sqrt(sig2) + + def eval_g(self, x: np.ndarray) -> np.ndarray: + mu = self.gpsurrogate.mean(x) + sig2 = self.gpsurrogate.variance(x) + dsig2_dx = self.gpsurrogate.variance_gradient(x) + dmu_dx = self.gpsurrogate.mean_gradient(x) + return dmu_dx - 0.5 * self.beta * dsig2_dx / np.sqrt(sig2) # A subclass of acquisition, implementing the Expected improvement (EI) acquisition function. class EIacquisition(acquisition): def __init__(self, gpsurrogate): super().__init__(gpsurrogate) + self.has_gradient = True # Method to evaluate the acquisition function at x. def evaluate(self, x : np.ndarray) -> np.ndarray: @@ -47,12 +60,40 @@ def evaluate(self, x : np.ndarray) -> np.ndarray: retval = [] if sig.size == 1 and np.abs(sig) > 1e-12: - arg0 = (y_min - pred) / sig - retval = (y_min - pred) * norm.cdf(arg0) + sig * norm.pdf(arg0) + z = (y_min - pred) / sig + retval = (y_min - pred) * norm.cdf(z) + sig * norm.pdf(z) retval *= -1. elif sig.size == 1 and np.abs(sig) <= 1e-12: retval = 0.0 elif sig.size > 1: - NotImplementedError("TODO --- Not implemented yet!") + raise NotImplementedError("TODO --- Not implemented yet!") return retval + + def eval_g(self, x: np.ndarray) -> np.ndarray: + y_data = self.gpsurrogate.training_y + y_min = y_data[np.argmin(y_data[:, 0])] + + mean = self.gpsurrogate.mean(x) + sig2 = self.gpsurrogate.variance(x) + sig = np.sqrt(sig2) + + grad_EI = None + if sig.size == 1 and np.abs(sig) > 1e-12: + dmean_dx = self.gpsurrogate.mean_gradient(x) + dsig2_dx = self.gpsurrogate.variance_gradient(x) + dsig_dx = 0.5 * dsig2_dx / sig + + z = (y_min - mean) / sig + ncdf = norm.cdf(z) + npdf = norm.pdf(z) + EI = (y_min - mean) * ncdf + sig * npdf + + dz_dx = -dmean_dx / sig - (y_min - mean) * dsig_dx / sig**2 + grad_EI = -dmean_dx * ncdf + dsig_dx * npdf + elif sig.size == 1 and np.abs(sig) <= 1e-12: + grad_EI = 0.0 + elif sig.size > 1: + raise NotImplementedError("TODO --- Not implemented yet!") + + return grad_EI \ No newline at end of file diff --git a/src/hiopbbpy/opt/boalgorithm.py b/src/hiopbbpy/opt/boalgorithm.py index ad882e4d7..a82b390ff 100644 --- a/src/hiopbbpy/opt/boalgorithm.py +++ b/src/hiopbbpy/opt/boalgorithm.py @@ -9,9 +9,11 @@ from numpy.random import uniform from scipy.optimize import minimize from scipy.stats import qmc +from cyipopt import Problem as IpoptProb from ..surrogate_modeling.gp import GaussianProcess from .acquisition import LCBacquisition, EIacquisition from ..problems.problem import Problem +from .optproblem import IpoptProbFromScipy # A base class defining a general framework for Bayesian Optimization class BOAlgorithmBase: @@ -20,7 +22,7 @@ def __init__(self): self.xtrain = None # Training data self.ytrain = None # Training data self.prob = None # Problem structure - self.n_iter = 20 # Maximum number of optimization steps + self.bo_maxiter = 20 # Maximum number of Bayesian optimization steps self.n_start = 10 # estimating acquisition global optima by determining local optima n_start times and then determining the discrete max of that set self.q = 1 # batch size # save some internal member train @@ -62,23 +64,53 @@ def getOptimalObjective(self): # A subclass of BOAlgorithmBase implementing a full Bayesian Optimization workflow class BOAlgorithm(BOAlgorithmBase): - def __init__(self, gpsurrogate, xtrain, ytrain, acquisition_type = "LCB"): + def __init__(self, gpsurrogate, xtrain, ytrain, + user_grad = None, + user_constraints = None, + options = None): super().__init__() + assert isinstance(gpsurrogate, GaussianProcess) - assert acquisition_type in ["LCB", "EI"] + self.setTrainingData(xtrain, ytrain) - self.setAcquisitionType(acquisition_type) self.gpsurrogate = gpsurrogate - self.n_iter = 20 - self.method = "SLSQP" self.bounds = self.gpsurrogate.get_bounds() - self.constraints = () - self.options = {"maxiter": 200} - self.acqf_minimizer_callback = None + self.fun_grad = None + self.constraints = None + + if options and 'bo_maxiter' in options: + self.bo_maxiter = options['bo_maxiter'] + assert self.bo_maxiter > 0, f"Invalid bo_maxiter: {self.bo_maxiter }" + + if options and 'solver_options' in options: + self.options = options['solver_options'] + else: + self.options = {"maxiter": 200} + + if options and 'acquisition_type' in options: + acquisition_type = options['acquisition_type'] + assert acquisition_type in ["LCB", "EI"], f"Invalid acquisition_type: {acquisition_type}" + else: + acquisition_type = "LCB" + self.setAcquisitionType(acquisition_type) + + if options and 'opt_solver' in options: + opt_solver = options['opt_solver'] + assert opt_solver in ["SLSQP", "IPOPT"], f"Invalid opt_solver: {opt_solver}" + else: + opt_solver = "SLSQP" + self.set_method(opt_solver) + + if user_constraints: + self.constraints = user_constraints + + if user_grad: + self.fun_grad = user_grad + # Method to set up a callback function to minimize the acquisition function def _setup_acqf_minimizer_callback(self): - self.acqf_minimizer_callback = lambda fun, x0: pyminimize(fun, x0, self.method, self.bounds, self.constraints, self.options) + self.acqf_minimizer_callback = lambda fun, x0: pyminimize(fun, x0, self.opt_solver, self.bounds, self.constraints, self.options) # Method to train the GP model def _train_surrogate(self, x_train, y_train): @@ -95,18 +127,21 @@ def _find_best_point(self, x_train, y_train, x0 = None): else: raise NotImplementedError("No implemented acquisition_type associated to"+self.acquisition_type) - acqf_callback = lambda x: float(np.array(acqf.evaluate(np.atleast_2d(x))).flat[0]) - + acqf_obj_callback = lambda x: float(np.array(acqf.evaluate(np.atleast_2d(x))).flat[0]) + acqf_callback = {'obj': acqf_obj_callback} + if acqf.has_gradient == True: + acqf_grad_callback = lambda x: np.array(acqf.eval_g(np.atleast_2d(x))) + acqf_callback['grad'] = acqf_grad_callback + x_all = [] y_all = [] for ii in range(self.n_start): success = False # Generate random starting point if x0 is not provided if x0 is None and self.prob is not None: - x0 = self.prob.sample(1) + x0 = self.prob.sample(1)[0] else: x0 = np.array([uniform(b[0], b[1]) for b in self.bounds]) - xopt, yout, success = self.acqf_minimizer_callback(acqf_callback, x0) if success: @@ -119,7 +154,7 @@ def _find_best_point(self, x_train, y_train, x0 = None): # Set the optimization method def set_method(self, method): - self.method = method + self.opt_solver = method # Set the user options for the Bayesian optimization def set_options(self, options): @@ -127,7 +162,7 @@ def set_options(self, options): # Method to perform Bayesian optimization def optimize(self, prob:Problem): - self.problem = prob + self.prob = prob x_train = self.xtrain y_train = self.ytrain @@ -138,9 +173,9 @@ def optimize(self, prob:Problem): self.x_hist = [] self.y_hist = [] - for i in range(self.n_iter): + for i in range(self.bo_maxiter): print(f"*****************************") - print(f"Iteration {i+1}/{self.n_iter}") + print(f"Iteration {i+1}/{self.bo_maxiter}") # Get a new sample point x_new = self._find_best_point(x_train, y_train) @@ -159,29 +194,41 @@ def optimize(self, prob:Problem): print(f"Sampled point X: {x_new.flatten()}, Observation Y: {y_new.flatten()}") # Save the optimal results and all the training data - self.idx_opt = np.argmin(y_train) - self.x_opt = x_train[self.idx_opt] - self.y_opt = y_train[self.idx_opt] + self.idx_opt = np.argmin(self.y_hist) + self.x_opt = self.y_hist[self.idx_opt] + self.y_opt = self.y_hist[self.idx_opt] self.setTrainingData(x_train, y_train) print() print() - if self.idx_opt < n_init_sample: - print(f"Optimal at initial sample: {self.idx_opt+1}") - else: - print(f"Optimal at BO iteration: {self.idx_opt-n_init_sample+1} ") + print(f"Optimal at BO iteration: {self.idx_opt+1} ") + #if self.idx_opt < n_init_sample: + # print(f"Optimal at initial sample: {self.idx_opt+1}") + #else: + # print(f"Optimal at BO iteration: {self.idx_opt-n_init_sample+1} ") print(f"Optimal point: {self.x_opt.flatten()}, Optimal value: {self.y_opt}") print() - - # Find the minimum of the input objective `fun`, using the minimize function from SciPy. def pyminimize(fun, x0, method, bounds, constraints, options): - y = minimize(fun, x0, method=method, bounds=bounds, constraints=constraints, options=options) - success = y.success - if not success: - print(y.message) - xopt = y.x - yopt = y.fun + if method != "IPOPT": + if 'grad' in fun: + y = minimize(fun['obj'], x0, method=method, bounds=bounds, jac=fun['grad'], constraints=constraints, options=options) + else: + y = minimize(fun['obj'], x0, method=method, bounds=bounds, constraints=constraints, options=options) + success = y.success + if not success: + print(y.message) + xopt = y.x + yopt = y.fun + else: + ipopt_prob = IpoptProbFromScipy(fun['obj'], fun['grad'], constraints, bounds) + nlp = IpoptProb(n=ipopt_prob.nvar, m=ipopt_prob.ncon, problem_obj=ipopt_prob, lb=ipopt_prob.xl, ub=ipopt_prob.xu, cl=ipopt_prob.cl, cu=ipopt_prob.cu) + + sol = nlp.solve(x0) + success = not sol[1]['status'] # ipopt returns 0 as success + yopt = sol[1]['obj_val'] + xopt = sol[0] + return xopt, yopt, success diff --git a/src/hiopbbpy/opt/optproblem.py b/src/hiopbbpy/opt/optproblem.py new file mode 100644 index 000000000..ef251dcb5 --- /dev/null +++ b/src/hiopbbpy/opt/optproblem.py @@ -0,0 +1,53 @@ +import numpy as np + +""" +Convert a Scipy optimization problem to an Ipopt problem. + +Parameters: + objective (callable): Objective function. + gradient (callable): Gradient of the objective. + constraints_list: Scipy-styple list of dicts with 'type', 'fun', and optional 'jac'. + xbounds (list of tuple): Variable bounds [(x0_lb, x0_ub), ...] + +Returns: + Ipopt-compatible prob and bounds +""" +class IpoptProbFromScipy: + def __init__(self, objective, gradient, constraints_list, xbounds): + self.constraints_list = constraints_list + self.eval_f = objective + self.eval_g = gradient + self.xl = [b[0] for b in xbounds] + self.xu = [b[1] for b in xbounds] + self.cl = [] + self.cu = [] + self.nvar = len(xbounds) + self.ncon = len(self.constraints_list) + + for con in constraints_list: + if con['type'] == 'eq': + self.cl.append(0.0) + self.cu.append(0.0) + elif con['type'] == 'ineq': + self.cl.append(0.0) + self.cu.append(np.inf) + else: + raise ValueError(f"Unknown constraint type: {con['type']}") + + def objective(self, x): + return self.eval_f(x) + + def gradient(self, x): + return self.eval_g(x) + + def constraints(self, x): + return np.array([con['fun'](x) for con in self.constraints_list]) + + def jacobian(self, x): + jacs = [] + for con in self.constraints_list: + if 'jac' in con: + jacs.append(con['jac'](x)) + else: + raise ValueError("Jacobian not provided for constraint.") + return np.vstack(jacs) diff --git a/src/hiopbbpy/surrogate_modeling/gp.py b/src/hiopbbpy/surrogate_modeling/gp.py index 9663a75d0..e1cf1955b 100644 --- a/src/hiopbbpy/surrogate_modeling/gp.py +++ b/src/hiopbbpy/surrogate_modeling/gp.py @@ -86,4 +86,36 @@ def train(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: y : ndarray[n, 1] """ - NotImplementedError("Child class of GaussianProcess should implement method train") \ No newline at end of file + NotImplementedError("Child class of GaussianProcess should implement method train") + + # Abstract method for computing the gradient of the mean of the GP at a given input x + def mean_gradient(self, x: np.ndarray) -> np.ndarray: + """ + evaluation of the gradien of GP mean + + Parameters + --------- + x : ndarray[n, nx] + + Returns + ------- + ndarray[n, 1] + Gradient of Mean of GP at x + """ + raise NotImplementedError("Child class of GaussianProcess should implement method mean_gradient") + + # Abstract method for computing the gradient of variance of the GP at a given input x + def variance_gradient(self, x: np.ndarray) -> np.ndarray: + """ + evaluation of the gradient of GP variance + + Parameters + --------- + x: ndarray[n, nx] + + Returns + ------- + ndarray[n, n] + Gradient of variance of GP at w.r.t. x + """ + raise NotImplementedError("Child class of GaussianProcess should implement method variance_gradient") diff --git a/src/hiopbbpy/surrogate_modeling/krg.py b/src/hiopbbpy/surrogate_modeling/krg.py index 8cdb6e787..a9227b7a4 100644 --- a/src/hiopbbpy/surrogate_modeling/krg.py +++ b/src/hiopbbpy/surrogate_modeling/krg.py @@ -32,12 +32,12 @@ def __init__(self, theta, xlimits, ndim, corr="pow_exp", noise0=None, random_sta def mean(self, x): if not self.trained: - ValueError("must train kriging model before utilizing it to predict mean or variances") + raise ValueError("must train kriging model before utilizing it to predict mean or variances") return self.surrogatesmt.predict_values(x) def variance(self, x): if not self.trained: - ValueError("must train kriging model before utilizing it to predict mean or variances") + raise ValueError("must train kriging model before utilizing it to predict mean or variances") return self.surrogatesmt.predict_variances(x) def train(self, x, y): @@ -47,7 +47,16 @@ def train(self, x, y): self.surrogatesmt.train() self.trained = True - - - - + def mean_gradient(self, x: np.ndarray) -> np.ndarray: + if not self.trained: + raise ValueError("must train kriging model before utilizing it to predict gradient") + assert (np.size(x,-1) == self.ndim) + gradient = [ + self.surrogatesmt._predict_derivatives(x, kx) for kx in range(self.ndim) + ] + return np.atleast_2d(gradient).T + + def variance_gradient(self, x: np.ndarray) -> np.ndarray: + if not self.trained: + raise ValueError("must train kriging model before utilizing it to predict gradient") + return self.surrogatesmt.predict_variance_gradient(x) From 594b756e7a2965f269be7103e0660d942eac93c5 Mon Sep 17 00:00:00 2001 From: nychiang Date: Wed, 9 Apr 2025 09:20:09 -0700 Subject: [PATCH 2/5] fix CI redisgn IpoptProb --- .github/workflows/test_hiopbbpy.yml | 30 +++++++++++++++--- pyproject.toml | 2 +- setup.py | 8 ++++- src/Drivers/hiopbbpy/BODriverCI.py | 11 ++++++- src/Drivers/hiopbbpy/BODriverEX.py | 3 +- src/hiopbbpy/opt/__init__.py | 4 +-- src/hiopbbpy/opt/acquisition.py | 1 + src/hiopbbpy/opt/boalgorithm.py | 48 +++++++++++++++++------------ src/hiopbbpy/opt/optproblem.py | 27 ++++++++++++++-- 9 files changed, 102 insertions(+), 32 deletions(-) diff --git a/.github/workflows/test_hiopbbpy.yml b/.github/workflows/test_hiopbbpy.yml index 2d5117bf3..32bbf7f0c 100644 --- a/.github/workflows/test_hiopbbpy.yml +++ b/.github/workflows/test_hiopbbpy.yml @@ -14,19 +14,39 @@ jobs: uses: actions/checkout@v4 - name: Set up Python - uses: actions/setup-python@v4 # Set up Python environment + uses: actions/setup-python@v4 with: - python-version: '3.13' + python-version: "3.13" + - name: Set up Conda + uses: conda-incubator/setup-miniconda@v2 + with: + auto-activate-base: true + channels: conda-forge + + - name: Install IPOPT via conda + run: | + conda init bash + conda install -c conda-forge cyipopt pkg-config + export LD_LIBRARY_PATH=$CONDA_PREFIX/lib:$LD_LIBRARY_PATH - - name: Install HiOpBBPy and its dependencies + - name: Verify cyipopt installation run: | + source ~/.bashrc + conda activate base + python -c "import cyipopt; print('cyipopt OK')" + + - name: Install HiOpBBPy + run: | + source ~/.bashrc + conda activate base + export SKIP_CYIPOPT=1 python -m pip install --upgrade pip pip install . - name: Run Tests run: | + source ~/.bashrc + conda activate base python src/Drivers/hiopbbpy/BODriverCI.py 10 continue-on-error: false - - diff --git a/pyproject.toml b/pyproject.toml index 955ab2bf7..424783deb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,3 @@ [build-system] -requires = ["setuptools", "wheel", "numpy", "smt", "cyipopt"] +requires = ["setuptools", "wheel", "numpy", "smt"] build-backend = "setuptools.build_meta" \ No newline at end of file diff --git a/setup.py b/setup.py index 903b554f4..c211af059 100644 --- a/setup.py +++ b/setup.py @@ -6,10 +6,16 @@ ''' import sys +import os import numpy as np from setuptools import setup, find_packages +install_requires = ["smt"] + +if os.getenv("SKIP_CYIPOPT", "0") != "1": + install_requires.append("cyipopt") + metadata = dict( name="hiopbbpy", version="0.0.4", @@ -19,7 +25,7 @@ license="BSD-3", packages=find_packages(where="src"), package_dir={"": "src"}, - install_requires=["smt","cyipopt"], + install_requires=install_requires, python_requires=">=3.9", zip_safe=False, url="https://github.com/LLNL/hiop", diff --git a/src/Drivers/hiopbbpy/BODriverCI.py b/src/Drivers/hiopbbpy/BODriverCI.py index d05f53161..3f6df8ac8 100644 --- a/src/Drivers/hiopbbpy/BODriverCI.py +++ b/src/Drivers/hiopbbpy/BODriverCI.py @@ -70,6 +70,15 @@ min_obj[prob_type][acq_type] = np.inf y_opt[prob_type][acq_type] = np.zeros(num_repeat) + options = { + 'acquisition_type': acq_type, + 'bo_maxiter': 20, + 'opt_solver': 'SLSQP', #"SLSQP" "IPOPT" + 'solver_options': { + 'maxiter': 200 + } + } + print("Problem name: ", problem.name) print("Acquisition type: ", acq_type) @@ -84,7 +93,7 @@ gp_model.train(x_train, y_train) # Instantiate and run Bayesian Optimization - bo = BOAlgorithm(gp_model, x_train, y_train, acquisition_type = acq_type) #EI or LCB + bo = BOAlgorithm(gp_model, x_train, y_train, options = options) #EI or LCB bo.optimize(problem) # Retrieve optimal objec diff --git a/src/Drivers/hiopbbpy/BODriverEX.py b/src/Drivers/hiopbbpy/BODriverEX.py index c9841008c..c618a3586 100644 --- a/src/Drivers/hiopbbpy/BODriverEX.py +++ b/src/Drivers/hiopbbpy/BODriverEX.py @@ -78,7 +78,8 @@ def con_jac_ineq(x): 'bo_maxiter': 10, 'opt_solver': 'IPOPT', #"SLSQP" "IPOPT" 'solver_options': { - 'maxiter': 200 + 'max_iter': 200, + 'print_level': 1 } } diff --git a/src/hiopbbpy/opt/__init__.py b/src/hiopbbpy/opt/__init__.py index 88cfcdb82..5692e3515 100644 --- a/src/hiopbbpy/opt/__init__.py +++ b/src/hiopbbpy/opt/__init__.py @@ -1,6 +1,6 @@ from .boalgorithm import (BOAlgorithmBase, BOAlgorithm) from .acquisition import (acquisition, LCBacquisition, EIacquisition) -from .optproblem import (IpoptProbFromScipy) +from .optproblem import (IpoptProb) __all__ = [ "BOAlgorithmBase" @@ -8,5 +8,5 @@ "acquisition" "LCBacquisition" "EIacquisition" - "IpoptProbFromScipy" + "IpoptProb" ] diff --git a/src/hiopbbpy/opt/acquisition.py b/src/hiopbbpy/opt/acquisition.py index ed8ef834f..07a496e76 100644 --- a/src/hiopbbpy/opt/acquisition.py +++ b/src/hiopbbpy/opt/acquisition.py @@ -91,6 +91,7 @@ def eval_g(self, x: np.ndarray) -> np.ndarray: dz_dx = -dmean_dx / sig - (y_min - mean) * dsig_dx / sig**2 grad_EI = -dmean_dx * ncdf + dsig_dx * npdf + grad_EI *= -1. elif sig.size == 1 and np.abs(sig) <= 1e-12: grad_EI = 0.0 elif sig.size > 1: diff --git a/src/hiopbbpy/opt/boalgorithm.py b/src/hiopbbpy/opt/boalgorithm.py index a82b390ff..834cb798e 100644 --- a/src/hiopbbpy/opt/boalgorithm.py +++ b/src/hiopbbpy/opt/boalgorithm.py @@ -9,11 +9,11 @@ from numpy.random import uniform from scipy.optimize import minimize from scipy.stats import qmc -from cyipopt import Problem as IpoptProb +import warnings from ..surrogate_modeling.gp import GaussianProcess from .acquisition import LCBacquisition, EIacquisition from ..problems.problem import Problem -from .optproblem import IpoptProbFromScipy +from .optproblem import IpoptProb # A base class defining a general framework for Bayesian Optimization class BOAlgorithmBase: @@ -83,9 +83,9 @@ def __init__(self, gpsurrogate, xtrain, ytrain, assert self.bo_maxiter > 0, f"Invalid bo_maxiter: {self.bo_maxiter }" if options and 'solver_options' in options: - self.options = options['solver_options'] + self.solver_options = options['solver_options'] else: - self.options = {"maxiter": 200} + self.solver_options = {"maxiter": 200} if options and 'acquisition_type' in options: acquisition_type = options['acquisition_type'] @@ -110,7 +110,7 @@ def __init__(self, gpsurrogate, xtrain, ytrain, # Method to set up a callback function to minimize the acquisition function def _setup_acqf_minimizer_callback(self): - self.acqf_minimizer_callback = lambda fun, x0: pyminimize(fun, x0, self.opt_solver, self.bounds, self.constraints, self.options) + self.acqf_minimizer_callback = lambda fun, x0: minimizer(fun, x0, self.opt_solver, self.bounds, self.constraints, self.solver_options) # Method to train the GP model def _train_surrogate(self, x_train, y_train): @@ -143,11 +143,14 @@ def _find_best_point(self, x_train, y_train, x0 = None): else: x0 = np.array([uniform(b[0], b[1]) for b in self.bounds]) xopt, yout, success = self.acqf_minimizer_callback(acqf_callback, x0) - + if success: x_all.append(xopt) y_all.append(yout) - + + if not x_all: + raise RuntimeError("Optimization failed for all initial points — no solution found.") + best_xopt = x_all[np.argmin(np.array(y_all))] return best_xopt @@ -156,9 +159,9 @@ def _find_best_point(self, x_train, y_train, x0 = None): def set_method(self, method): self.opt_solver = method - # Set the user options for the Bayesian optimization - def set_options(self, options): - self.options = options + # Set the options for the internal optimization solver + def set_options(self, solver_options): + self.solver_options = solver_options # Method to perform Bayesian optimization def optimize(self, prob:Problem): @@ -211,24 +214,31 @@ def optimize(self, prob:Problem): print() # Find the minimum of the input objective `fun`, using the minimize function from SciPy. -def pyminimize(fun, x0, method, bounds, constraints, options): +def minimizer(fun, x0, method, bounds, constraints, solver_options): if method != "IPOPT": if 'grad' in fun: - y = minimize(fun['obj'], x0, method=method, bounds=bounds, jac=fun['grad'], constraints=constraints, options=options) + y = minimize(fun['obj'], x0, method=method, bounds=bounds, jac=fun['grad'], constraints=constraints, options=solver_options) else: - y = minimize(fun['obj'], x0, method=method, bounds=bounds, constraints=constraints, options=options) + y = minimize(fun['obj'], x0, method=method, bounds=bounds, constraints=constraints, options=solver_options) success = y.success if not success: print(y.message) xopt = y.x yopt = y.fun else: - ipopt_prob = IpoptProbFromScipy(fun['obj'], fun['grad'], constraints, bounds) - nlp = IpoptProb(n=ipopt_prob.nvar, m=ipopt_prob.ncon, problem_obj=ipopt_prob, lb=ipopt_prob.xl, ub=ipopt_prob.xu, cl=ipopt_prob.cl, cu=ipopt_prob.cu) + ipopt_prob = IpoptProb(fun['obj'], fun['grad'], constraints, bounds, solver_options) + sol, info = ipopt_prob.solve(x0) + + status = info.get('status', -1) + msg = info.get('status_msg', -1) + if status == 0: + # ipopt returns 0 as success + success = True + else: + warnings.warn(f"Ipopt failed to solve the problem. Status msg: {msg}") + success = False - sol = nlp.solve(x0) - success = not sol[1]['status'] # ipopt returns 0 as success - yopt = sol[1]['obj_val'] - xopt = sol[0] + yopt = info['obj_val'] + xopt = sol return xopt, yopt, success diff --git a/src/hiopbbpy/opt/optproblem.py b/src/hiopbbpy/opt/optproblem.py index ef251dcb5..a131fcba2 100644 --- a/src/hiopbbpy/opt/optproblem.py +++ b/src/hiopbbpy/opt/optproblem.py @@ -1,4 +1,5 @@ import numpy as np +import cyipopt """ Convert a Scipy optimization problem to an Ipopt problem. @@ -12,8 +13,8 @@ Returns: Ipopt-compatible prob and bounds """ -class IpoptProbFromScipy: - def __init__(self, objective, gradient, constraints_list, xbounds): +class IpoptProb: + def __init__(self, objective, gradient, constraints_list, xbounds, solver_options=None): self.constraints_list = constraints_list self.eval_f = objective self.eval_g = gradient @@ -23,6 +24,7 @@ def __init__(self, objective, gradient, constraints_list, xbounds): self.cu = [] self.nvar = len(xbounds) self.ncon = len(self.constraints_list) + self.ipopt_options = solver_options for con in constraints_list: if con['type'] == 'eq': @@ -34,6 +36,16 @@ def __init__(self, objective, gradient, constraints_list, xbounds): else: raise ValueError(f"Unknown constraint type: {con['type']}") + self.nlp = cyipopt.Problem( + n=self.nvar, + m=self.ncon, + problem_obj=self, + lb=self.xl, + ub=self.xu, + cl=self.cl, + cu=self.xu + ) + def objective(self, x): return self.eval_f(x) @@ -51,3 +63,14 @@ def jacobian(self, x): else: raise ValueError("Jacobian not provided for constraint.") return np.vstack(jacs) + + def solve(self, x0, solver_options=None): + ipopt_options = self.ipopt_options + if solver_options is not None: + ipopt_options = solver_options + if ipopt_options is not None: + for key, value in ipopt_options.items(): + self.nlp.add_option(key, value) + + # Solve the optimization problem + return self.nlp.solve(x0) \ No newline at end of file From a58037af7de7358502f98d80bfc0165e743d0ef0 Mon Sep 17 00:00:00 2001 From: nychiang Date: Wed, 9 Apr 2025 16:44:27 -0700 Subject: [PATCH 3/5] add authors --- src/hiopbbpy/opt/optproblem.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/hiopbbpy/opt/optproblem.py b/src/hiopbbpy/opt/optproblem.py index a131fcba2..cd0745544 100644 --- a/src/hiopbbpy/opt/optproblem.py +++ b/src/hiopbbpy/opt/optproblem.py @@ -12,6 +12,9 @@ Returns: Ipopt-compatible prob and bounds + +Authors: Tucker Hartland + Nai-Yuan Chiang """ class IpoptProb: def __init__(self, objective, gradient, constraints_list, xbounds, solver_options=None): From 38762aff192a91a8ce6c0271471cdbfd303d0593 Mon Sep 17 00:00:00 2001 From: nychiang Date: Wed, 9 Apr 2025 23:15:55 -0700 Subject: [PATCH 4/5] move constraints into Problem class --- src/Drivers/hiopbbpy/BODriverEX.py | 9 +++++---- src/Drivers/hiopbbpy/LpNormProblem.py | 4 ++-- src/hiopbbpy/opt/boalgorithm.py | 19 +++++++------------ src/hiopbbpy/opt/optproblem.py | 4 ++-- src/hiopbbpy/problems/BraninProblem.py | 13 +++---------- src/hiopbbpy/problems/problem.py | 6 +++++- 6 files changed, 24 insertions(+), 31 deletions(-) diff --git a/src/Drivers/hiopbbpy/BODriverEX.py b/src/Drivers/hiopbbpy/BODriverEX.py index c618a3586..daf3b2720 100644 --- a/src/Drivers/hiopbbpy/BODriverEX.py +++ b/src/Drivers/hiopbbpy/BODriverEX.py @@ -57,9 +57,10 @@ def con_jac_ineq(x): for prob_type in prob_type_l: print() if prob_type == "LpNorm": - problem = LpNormProblem(nx, xlimits) + problem = LpNormProblem(nx, xlimits, constraints=None) else: - problem = BraninProblem() + problem = BraninProblem(constraints=None) + problem.set_constraints(user_constraint) for acq_type in acq_type_l: print("Problem name: ", problem.name) @@ -84,8 +85,8 @@ def con_jac_ineq(x): } # Instantiate and run Bayesian Optimization - bo = BOAlgorithm(gp_model, x_train, y_train, options = options, user_constraints = user_constraint) #EI or LCB - bo.optimize(problem) + bo = BOAlgorithm(problem, gp_model, x_train, y_train, options = options) #EI or LCB + bo.optimize() sys.exit(retval) diff --git a/src/Drivers/hiopbbpy/LpNormProblem.py b/src/Drivers/hiopbbpy/LpNormProblem.py index 2ce3b2fc1..3084bd802 100644 --- a/src/Drivers/hiopbbpy/LpNormProblem.py +++ b/src/Drivers/hiopbbpy/LpNormProblem.py @@ -8,9 +8,9 @@ from hiopbbpy.problems.problem import Problem class LpNormProblem(Problem): - def __init__(self, ndim, xlimits, p=2.0): + def __init__(self, ndim, xlimits, p=2.0, constraints=None): name = "LpNormProblem" - super().__init__(ndim, xlimits, name=name) + super().__init__(ndim, xlimits, name=name, constraints=constraints) self.p = p def _evaluate(self, x): diff --git a/src/hiopbbpy/opt/boalgorithm.py b/src/hiopbbpy/opt/boalgorithm.py index 834cb798e..081572f41 100644 --- a/src/hiopbbpy/opt/boalgorithm.py +++ b/src/hiopbbpy/opt/boalgorithm.py @@ -64,19 +64,18 @@ def getOptimalObjective(self): # A subclass of BOAlgorithmBase implementing a full Bayesian Optimization workflow class BOAlgorithm(BOAlgorithmBase): - def __init__(self, gpsurrogate, xtrain, ytrain, + def __init__(self, prob:Problem, gpsurrogate:GaussianProcess, xtrain, ytrain, user_grad = None, - user_constraints = None, options = None): super().__init__() assert isinstance(gpsurrogate, GaussianProcess) self.setTrainingData(xtrain, ytrain) + self.prob = prob self.gpsurrogate = gpsurrogate self.bounds = self.gpsurrogate.get_bounds() self.fun_grad = None - self.constraints = None if options and 'bo_maxiter' in options: self.bo_maxiter = options['bo_maxiter'] @@ -101,16 +100,13 @@ def __init__(self, gpsurrogate, xtrain, ytrain, opt_solver = "SLSQP" self.set_method(opt_solver) - if user_constraints: - self.constraints = user_constraints - if user_grad: self.fun_grad = user_grad # Method to set up a callback function to minimize the acquisition function def _setup_acqf_minimizer_callback(self): - self.acqf_minimizer_callback = lambda fun, x0: minimizer(fun, x0, self.opt_solver, self.bounds, self.constraints, self.solver_options) + self.acqf_minimizer_callback = lambda fun, x0: minimizer(fun, x0, self.opt_solver, self.bounds, self.prob.constraints, self.solver_options) # Method to train the GP model def _train_surrogate(self, x_train, y_train): @@ -164,8 +160,7 @@ def set_options(self, solver_options): self.solver_options = solver_options # Method to perform Bayesian optimization - def optimize(self, prob:Problem): - self.prob = prob + def optimize(self): x_train = self.xtrain y_train = self.ytrain @@ -184,7 +179,7 @@ def optimize(self, prob:Problem): x_new = self._find_best_point(x_train, y_train) # Evaluate the new sample point - y_new = prob.evaluate(np.atleast_2d(x_new)) + y_new = self.prob.evaluate(np.atleast_2d(x_new)) # Update training set x_train = np.vstack([x_train, x_new]) @@ -229,8 +224,8 @@ def minimizer(fun, x0, method, bounds, constraints, solver_options): ipopt_prob = IpoptProb(fun['obj'], fun['grad'], constraints, bounds, solver_options) sol, info = ipopt_prob.solve(x0) - status = info.get('status', -1) - msg = info.get('status_msg', -1) + status = info.get('status', -999) + msg = info.get('status_msg', b'unknown error') if status == 0: # ipopt returns 0 as success success = True diff --git a/src/hiopbbpy/opt/optproblem.py b/src/hiopbbpy/opt/optproblem.py index cd0745544..0f4573aca 100644 --- a/src/hiopbbpy/opt/optproblem.py +++ b/src/hiopbbpy/opt/optproblem.py @@ -46,7 +46,7 @@ def __init__(self, objective, gradient, constraints_list, xbounds, solver_option lb=self.xl, ub=self.xu, cl=self.cl, - cu=self.xu + cu=self.cu ) def objective(self, x): @@ -76,4 +76,4 @@ def solve(self, x0, solver_options=None): self.nlp.add_option(key, value) # Solve the optimization problem - return self.nlp.solve(x0) \ No newline at end of file + return self.nlp.solve(x0) diff --git a/src/hiopbbpy/problems/BraninProblem.py b/src/hiopbbpy/problems/BraninProblem.py index a4485b22e..d45a7575f 100644 --- a/src/hiopbbpy/problems/BraninProblem.py +++ b/src/hiopbbpy/problems/BraninProblem.py @@ -19,11 +19,11 @@ # define the Branin problem class class BraninProblem(Problem): - def __init__(self): + def __init__(self, constraints=None): ndim = 2 xlimits = np.array([[-5.0, 10], [0.0, 15]]) name = 'Branin' - super().__init__(ndim, xlimits, name=name) + super().__init__(ndim, xlimits, name=name, constraints=constraints) def _evaluate(self, x: np.ndarray) -> np.ndarray: @@ -40,17 +40,10 @@ def _evaluate(self, x: np.ndarray) -> np.ndarray: arg1 = (x[:,1] - b * x[:,0]**2 + c * x[:,0] - r) y[:,0] = arg1**2 + s * (1 - t) * np.cos(x[:,0]) + s - ''' - # compute derivatives - dy_dx0 = 2*arg1*(-2*b*x[:,0]+c) - s*(1-t)*np.sin(x[:,0]) - dy_dx1 = 2*arg1 - - dy_dx = np.array([dy_dx0, dy_dx1]) - ''' - return y + diff --git a/src/hiopbbpy/problems/problem.py b/src/hiopbbpy/problems/problem.py index c0895ae33..1d0bae18f 100644 --- a/src/hiopbbpy/problems/problem.py +++ b/src/hiopbbpy/problems/problem.py @@ -10,7 +10,7 @@ # define the general optimization problem class class Problem: - def __init__(self, ndim, xlimits, name=None): + def __init__(self, ndim, xlimits, name=None, constraints=None): self.ndim = ndim self.xlimits = xlimits assert self.xlimits.shape[0] == ndim @@ -19,6 +19,7 @@ def __init__(self, ndim, xlimits, name=None): else: self.name = name self.sampler = qmc.LatinHypercube(ndim) + self.constraints = constraints def _evaluate(self, x: np.ndarray) -> np.ndarray: """ @@ -80,6 +81,9 @@ def sample(self, nsample: int) -> np.ndarray: return xsample + def set_constraints(self, constraints: dict): + self.constraints = constraints + From 50e138888e47387e48a4acb682378f204237c40b Mon Sep 17 00:00:00 2001 From: nychiang Date: Wed, 9 Apr 2025 23:25:52 -0700 Subject: [PATCH 5/5] update CI --- src/Drivers/hiopbbpy/BODriverCI.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Drivers/hiopbbpy/BODriverCI.py b/src/Drivers/hiopbbpy/BODriverCI.py index 3f6df8ac8..2bcc386f2 100644 --- a/src/Drivers/hiopbbpy/BODriverCI.py +++ b/src/Drivers/hiopbbpy/BODriverCI.py @@ -93,8 +93,8 @@ gp_model.train(x_train, y_train) # Instantiate and run Bayesian Optimization - bo = BOAlgorithm(gp_model, x_train, y_train, options = options) #EI or LCB - bo.optimize(problem) + bo = BOAlgorithm(problem, gp_model, x_train, y_train, options = options) + bo.optimize() # Retrieve optimal objec y_opt[prob_type][acq_type][n_repeat] = bo.getOptimalObjective()