redisgn IpoptProb

Author: nychiang

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
            }
       }
 

src/hiopbbpy/opt/__init__.py

@@ -1,12 +1,12 @@
 from .boalgorithm import (BOAlgorithmBase, BOAlgorithm)
 from .acquisition import (acquisition, LCBacquisition, EIacquisition)
-from .optproblem import (IpoptProbFromScipy)
+from .optproblem import (IpoptProb)
 
 __all__ = [
         "BOAlgorithmBase"
         "BOAlgorithm"
         "acquisition"
         "LCBacquisition"
         "EIacquisition"
-        "IpoptProbFromScipy"
+        "IpoptProb"
         ]

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)
         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

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)