From ffb185c9e41fadb28a06a0e93623c6059970d6c4 Mon Sep 17 00:00:00 2001 From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com> Date: Wed, 7 Aug 2024 10:58:41 +0100 Subject: [PATCH 1/6] One opt param and one constraint in solver test case 5 --- tests/integration/test_vmcon.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/integration/test_vmcon.py b/tests/integration/test_vmcon.py index c5613e008c..bed897d85e 100644 --- a/tests/integration/test_vmcon.py +++ b/tests/integration/test_vmcon.py @@ -578,6 +578,10 @@ def get_case5(): case.solver_args.x = np.array([ 5.0 ]) # Try different values, e.g. 5.0, 2.0, 1.0, 0.0... + case.solver_args.bndl = np.zeros(1) + case.solver_args.bndu = np.full(1, 5.0) + case.solver_args.ilower = np.zeros(1) + case.solver_args.iupper = np.zeros(1) # Expected values case.exp.x = np.array([3.0]) From 2dc41107fdf12187b807c39e146d4adf924e24bc Mon Sep 17 00:00:00 2001 From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com> Date: Wed, 7 Aug 2024 11:18:53 +0100 Subject: [PATCH 2/6] Implement nlopt's SLSQP solver --- process/solver.py | 194 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 194 insertions(+) diff --git a/process/solver.py b/process/solver.py index bffbfb07d7..19c793bdcf 100644 --- a/process/solver.py +++ b/process/solver.py @@ -14,6 +14,8 @@ solve, ) from scipy.optimize import fsolve +import nlopt +import time from process.evaluators import Evaluators from process.exceptions import ProcessValueError @@ -317,6 +319,196 @@ def solve(self) -> int: return self.info +class SLSQP(_Solver): + def obj_func(self, x, grad): + # Must be passed these args from nlopt requirements + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x, self.ifail) + self.constr_res = np.sqrt( + np.sum(np.square(conf[: self.meq])) + ) # only for equality constraints + + if grad.size > 0: + # Gradient required by solver; modify grad in-place + fgrd, cnorm = self.evaluators.fcnvmc2(self.n, self.m, x, self.lcnorm) + grad[...] = fgrd + + self.eval_count += 1 + status = ( + f"Evaluation {self.eval_count}, objective function = {objf:.5}, " + f"constraint residuals = {self.constr_res:.3e}" + ) + print(status) + logger.info(status) + return objf + + def constraint_eq_vec(self, result, x, grad): + # i is no. of constraint + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x, self.ifail) + self.constr_res = np.sqrt( + np.sum(np.square(conf[: self.meq])) + ) # only for equality constraints + + # Check if constraints are below tolerance yet + # conf_gt_tol = np.sum((np.abs(conf[:meq]) > CONSTR_TOL)) + # if conf_gt_tol == 0: + # logger.info("Constraints are all below tolerance") + # else: + # logger.info(f"Constraints above tolerance: {conf_gt_tol=}") + + if grad.size > 0: + # Gradient required by solver; modify grad in-place + fgrd, cnorm = self.evaluators.fcnvmc2(self.n, self.m, x, self.lcnorm) + if cnorm.ndim == 1: + # 1 constraint, 1 optimisation parameter (used in tests) + # TODO Add np.newaxis to fcnvmc2 to avoid this + grad[...] = -cnorm[:] + else: + grad[...] = np.transpose(-cnorm[: x.shape[0], : self.meq]) + + # Negative conf and cnorm: nlopt requires opposite formulation of + # VMCON's (and Process's) constraint form + result[...] = -conf[: self.meq] + + def constraint_ineq_vec(self, result, x, grad): + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x, self.ifail) + if grad.size > 0: + # Gradient required by solver; modify grad in-place + fgrd, cnorm = self.evaluators.fcnvmc2(self.n, self.m, x, self.lcnorm) + if cnorm.ndim == 1: + # 1 constraint, 1 optimisation parameter (used in tests) + # TODO Add np.newaxis to fcnvmc2 to avoid this + grad[...] = -cnorm[0] + else: + grad[...] = np.transpose(-cnorm[: x.shape[0], self.meq : self.m]) + + # Negative conf and cnorm: nlopt requires opposite formulation of + # VMCON's (and Process's) constraint form + result[...] = -conf[self.meq : self.m] + + def solve(self) -> int: + """Try running nlopt.""" + # global objf, conf, fgrd, cnorm + self.eval_count = 0 + self.constr_res = 0.0 + + self.n = self.x_0.shape[0] + self.lcnorm = 176 # (ippnvars + 1), but could just be n! + + # Solver tolerances (high) + MAIN_TOL = 1e-8 + # LOCAL_TOL = 1e-8 # for AUGLAG method only + CONSTR_TOL = 1e-10 + + # Set up optimiser object + # opt = nlopt.opt(nlopt.AUGLAG, n) + opt = nlopt.opt(nlopt.LD_SLSQP, self.n) + + # Set subsidiary optimiser + # For AUGLAG only + # local_opt = nlopt.opt(nlopt.LD_SLSQP, n) + + opt_name = opt.get_algorithm_name() + # local_opt_name = local_opt.get_algorithm_name() + # logger.info(f"{opt_name=}, {local_opt_name=}") + logger.info(f"{opt_name=}") + logger.info(f"{MAIN_TOL=:.3e}, {CONSTR_TOL=:.3e}") + # logger.info(f"{MAIN_TOL=:.3e}, {LOCAL_TOL=:.3e}, {CONSTR_TOL=:.3e}") + + # Define tolerances + opt.set_ftol_rel(MAIN_TOL) + # local_opt.set_ftol_rel(LOCAL_TOL) + + # Need to terminate in solver test case 3! + opt.set_maxeval(1000) + + # opt.set_local_optimizer(local_opt) + + # if self.maximise: + # # Maximisation required for test case 5 + # opt.set_max_objective(self.obj_func) + # else: + opt.set_min_objective(self.obj_func) + + # Check bounds are activated for all optimisation parameters (default case) + # If not, handle it + if not (np.all(self.ilower) and np.all(self.iupper)): + # Some bounds are inactive: used e.g. in solver integration tests + for i in range(self.ilower.shape[0]): + if self.ilower[i] == 0: + # Inactive lower bound: set to -inf + self.bndl[i] = -np.inf + + for i in range(self.iupper.shape[0]): + if self.iupper[i] == 0: + # Inactive upper bound: set to +inf + self.bndu[i] = np.inf + + opt.set_lower_bounds(self.bndl) + opt.set_upper_bounds(self.bndu) + + # Kludge initial normalised x into the normalised bounds range if required + for i in range(self.x_0.shape[0]): + if self.x_0[i] < self.bndl[i]: + self.x_0[i] = self.bndl[i] + elif self.x_0[i] > self.bndu[i]: + self.x_0[i] = self.bndu[i] + + # Constraints + if self.meq > 0: + eq_constr_tols = np.full(self.meq, CONSTR_TOL) + opt.add_equality_mconstraint(self.constraint_eq_vec, eq_constr_tols) + if self.meq < self.m: + ineq_constr_tols = np.full((self.m - self.meq), CONSTR_TOL) + opt.add_inequality_mconstraint(self.constraint_ineq_vec, ineq_constr_tols) + + start_time = time.time() + x_opt = opt.optimize(self.x_0) + end_time = time.time() + duration = end_time - start_time + + opt_val = opt.last_optimum_value() + return_value = opt.last_optimize_result() + + # Main opt iterations + main_opt_evals = opt.get_numevals() + + if return_value < 0: + raise RuntimeError("nlopt didn't converge") + elif return_value == nlopt.MAXEVAL_REACHED: + # For giving up in int test 3 + # TODO Reconcile MAXEVAL with above exception: int case 3 needs to pass + info = 5 + elif return_value > 0: + # TODO Might want to be aware of other return values + info = 1 + + print(f"{return_value=}") + + # Recalculate conf at optimum x + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x_opt, self.ifail) + self.constr_res = np.sqrt( + np.sum(np.square(conf[: self.meq])) + ) # only for equality constraints + + logger.info(f"Main opt evaluations = {main_opt_evals}") + logger.info(f"{opt_val=:.3e}, {self.constr_res=:.3e}") + logger.info(f"{duration=:.1f}\n") + logger.info(f"{conf[:self.meq]=}") + + print(f"{self.constr_res=:.3e}") + + # Check how many conf elements are above tolerance + conf_gt_tol = np.sum((np.abs(conf[: self.meq]) > CONSTR_TOL)) + logger.info(f"Constraints above tolerance: {conf_gt_tol}") + + # Store required results on object + self.objf = objf + self.conf = conf + self.x = x_opt + + return info + + def get_solver(solver_name: str = "vmcon") -> _Solver: """Return a solver instance. @@ -333,6 +525,8 @@ def get_solver(solver_name: str = "vmcon") -> _Solver: solver = VmconBounded() elif solver_name == "fsolve": solver = FSolve() + elif solver_name == "slsqp": + solver = SLSQP() else: try: solver = load_external_solver(solver_name) From 2cf86ca52f98863e74625577712ad255d339d6bd Mon Sep 17 00:00:00 2001 From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com> Date: Wed, 14 Aug 2024 09:50:58 +0100 Subject: [PATCH 3/6] Add scipy's SLSQP --- process/solver.py | 135 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 134 insertions(+), 1 deletion(-) diff --git a/process/solver.py b/process/solver.py index 19c793bdcf..34bb30327d 100644 --- a/process/solver.py +++ b/process/solver.py @@ -16,6 +16,7 @@ from scipy.optimize import fsolve import nlopt import time +from scipy import optimize from process.evaluators import Evaluators from process.exceptions import ProcessValueError @@ -387,7 +388,7 @@ def constraint_ineq_vec(self, result, x, grad): def solve(self) -> int: """Try running nlopt.""" - # global objf, conf, fgrd, cnorm + print("Using SLSQP") self.eval_count = 0 self.constr_res = 0.0 @@ -509,6 +510,136 @@ def solve(self) -> int: return info +class Scipy_SLSQP(_Solver): + """Minimise using scipy's SLSQP.""" + + print("Running scipy's SLSQP") + SOLVER_TOL = 1e-6 + EQ_CONSTRAINT_TOL = 1e-8 + + def obj_func(self, x): + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x, self.ifail) + # constr_res = np.sqrt( + # np.sum(np.square(conf[:meq])) + # ) # only for equality constraints + # logger.debug(f"Constraint residuals: {constr_res:.3e}") + + # print( + # f"Evaluation {eval_count}, objective function = {objf:.5}, " + # f"constraint residuals = {constr_res:.3e}" + # ) + return objf + + def constraint_eq_vec(self, x): + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x, self.ifail) + # constr_res = np.sqrt( + # np.sum(np.square(conf[:meq])) + # ) # only to equality constraints + + return conf[: self.meq] + + def constraint_ineq_vec(self, x): + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, x, self.ifail) + conf_gt_tol = np.sum((conf[self.meq :] < 0.0)) + logger.info(f"{conf_gt_tol} inequality constraints above 0.0") + return conf[self.meq : self.m] + + def convergence_progress(self, x_current): + conf = self.constraint_ineq_vec(x_current) + ineq_rms = np.sqrt(np.mean(np.square(conf[conf < 0.0]))) + print(f"{ineq_rms = :.3e}") + + def solve(self): + self.n = self.x_0.shape[0] + + # Check bounds are activated for all optimisation parameters (default case) + # If not, handle it + if not (np.all(self.ilower) and np.all(self.iupper)): + # Some bounds are inactive: used e.g. in solver integration tests + for i in range(self.ilower.shape[0]): + if self.ilower[i] == 0: + # Inactive lower bound: set to -inf + self.bndl[i] = -np.inf + + for i in range(self.iupper.shape[0]): + if self.iupper[i] == 0: + # Inactive upper bound: set to +inf + self.bndu[i] = np.inf + + # Kludge initial normalised x into the normalised bounds range if required + for i in range(self.n): + if self.x_0[i] < self.bndl[i]: + self.x_0[i] = self.bndl[i] + elif self.x_0[i] > self.bndu[i]: + self.x_0[i] = self.bndu[i] + + bounds = optimize.Bounds(lb=self.bndl, ub=self.bndu) + + constraints = [] + if self.meq > 0: + eq_constraints = optimize.NonlinearConstraint( + self.constraint_eq_vec, -self.EQ_CONSTRAINT_TOL, self.EQ_CONSTRAINT_TOL + ) + constraints.append(eq_constraints) + + if self.meq < self.m: + ineq_constraints = optimize.NonlinearConstraint( + self.constraint_ineq_vec, + 0.0, + np.inf, + ) + constraints.append(ineq_constraints) + + start_time = time.time() + + result = optimize.minimize( + self.obj_func, + self.x_0, + method="SLSQP", + jac=None, + bounds=bounds, + constraints=constraints, + tol=self.SOLVER_TOL, + callback=self.convergence_progress, + options={"disp": True, "eps": 1e-3}, + ) + end_time = time.time() + duration = end_time - start_time + + # Log stuff + logger.info(f"{self.SOLVER_TOL=}, {self.EQ_CONSTRAINT_TOL=}") + logger.info(f"Iterations: {result.nit}") + logger.info(f"Evaluations: {result.nfev}") + logger.info(f"Duration: {duration:.3}") + + # Recalculate constraints at optimium x + objf, conf = self.evaluators.fcnvmc1(self.n, self.m, result.x, self.ifail) + + # Check if constraints are all below tolerance + conf_gt_tol = np.sum((conf[self.meq :] < 0.0)) + logger.info(f"{conf_gt_tol} inequality constraints violated") + logger.info(f"Constraint residuals: {conf}") + + # constr_res = np.sqrt( + # np.sum(np.square(conf[: self.meq])) + # ) # only for equality constraints + # logger.info(f"Constraint residuals: {constr_res:.3e}") + # logger.info(f"{conf=}") + + if result.success: + info = 1 + else: + # Want to write error code to MFILE + # raise RuntimeError("scipy failed to converge") + info = 2 + + self.objf = result.fun + self.conf = conf + self.x = result.x + + return info + + def get_solver(solver_name: str = "vmcon") -> _Solver: """Return a solver instance. @@ -527,6 +658,8 @@ def get_solver(solver_name: str = "vmcon") -> _Solver: solver = FSolve() elif solver_name == "slsqp": solver = SLSQP() + elif solver_name == "scipy_slsqp": + solver = Scipy_SLSQP() else: try: solver = load_external_solver(solver_name) From 2b3f4a7982ca5f7ee69ed1e0480ccd9bfc71f76d Mon Sep 17 00:00:00 2001 From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com> Date: Thu, 22 Aug 2024 10:35:16 +0100 Subject: [PATCH 4/6] Use custom finite difference step for SLSQP --- process/solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/solver.py b/process/solver.py index 34bb30327d..5c9ce1ed69 100644 --- a/process/solver.py +++ b/process/solver.py @@ -601,7 +601,7 @@ def solve(self): constraints=constraints, tol=self.SOLVER_TOL, callback=self.convergence_progress, - options={"disp": True, "eps": 1e-3}, + options={"disp": True, "eps": numerics.epsfcn}, ) end_time = time.time() duration = end_time - start_time From 8b9c47d786d0500a5624f4db863b94e292d22cfd Mon Sep 17 00:00:00 2001 From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com> Date: Thu, 22 Aug 2024 10:35:44 +0100 Subject: [PATCH 5/6] Add constraint residual reporting to SLSQP --- process/solver.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/process/solver.py b/process/solver.py index 5c9ce1ed69..637dfbf569 100644 --- a/process/solver.py +++ b/process/solver.py @@ -545,9 +545,17 @@ def constraint_ineq_vec(self, x): return conf[self.meq : self.m] def convergence_progress(self, x_current): - conf = self.constraint_ineq_vec(x_current) - ineq_rms = np.sqrt(np.mean(np.square(conf[conf < 0.0]))) - print(f"{ineq_rms = :.3e}") + eqs = self.constraint_eq_vec(x_current) + ineqs = self.constraint_ineq_vec(x_current) + cons = np.concatenate((eqs, ineqs)) + print("\nIteration results:") + ineqs_rms = np.sqrt(np.mean(np.square(ineqs[ineqs < 0.0]))) + print(f"{ineqs_rms = :.3e}") + + # Print constraints sorted by value (most negative (most violated) first) + sorted_con_indexes = cons.argsort() + for i in sorted_con_indexes: + print(f"Constraint {numerics.icc[i]} = {cons[i]:.3e}") def solve(self): self.n = self.x_0.shape[0] From 75261167c36b239daa091452c7c62c0d65148a75 Mon Sep 17 00:00:00 2001 From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com> Date: Fri, 30 May 2025 15:02:51 +0100 Subject: [PATCH 6/6] Add nlopt to setup.py --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index 45e59451a6..1c3f520e56 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,7 @@ "matplotlib>=2.1.1", "seaborn>=0.12.2", "tabulate", + "nlopt", ], "extras_require": { "test": ["pytest>=5.4.1", "requests>=2.30", "testbook>=0.4"],