libMesh
nlopt_optimization_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #include "libmesh/libmesh_common.h"
21 
22 #if defined(LIBMESH_HAVE_NLOPT) && !defined(LIBMESH_USE_COMPLEX_NUMBERS)
23 
24 
25 // C++ includes
26 
27 // Local Includes
28 #include "libmesh/dof_map.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/nlopt_optimization_solver.h"
31 #include "libmesh/sparse_matrix.h"
32 
33 namespace libMesh
34 {
35 
36 double __libmesh_nlopt_objective(unsigned n,
37  const double * x,
38  double * gradient,
39  void * data)
40 {
41  LOG_SCOPE("objective()", "NloptOptimizationSolver");
42 
43  // ctx should be a pointer to the solver (it was passed in as void *)
45  static_cast<NloptOptimizationSolver<Number> *> (data);
46 
47  OptimizationSystem & sys = solver->system();
48 
49  // We'll use current_local_solution below, so let's ensure that it's consistent
50  // with the vector x that was passed in.
51  for (unsigned int i=sys.solution->first_local_index();
52  i<sys.solution->last_local_index(); i++)
53  sys.solution->set(i, x[i]);
54 
55  // Make sure the solution vector is parallel-consistent
56  sys.solution->close();
57 
58  // Impose constraints on X
59  sys.get_dof_map().enforce_constraints_exactly(sys);
60 
61  // Update sys.current_local_solution based on X
62  sys.update();
63 
64  Real objective;
65  if (solver->objective_object != libmesh_nullptr)
66  {
67  objective =
68  solver->objective_object->objective(
69  *(sys.current_local_solution), sys);
70  }
71  else
72  {
73  libmesh_error_msg("Objective function not defined in __libmesh_nlopt_objective");
74  }
75 
76  // If the gradient has been requested, fill it in
77  if (gradient)
78  {
79  if (solver->gradient_object != libmesh_nullptr)
80  {
81  solver->gradient_object->gradient(
82  *(sys.current_local_solution), *(sys.rhs), sys);
83 
84  // we've filled up sys.rhs with the gradient data, now copy it
85  // to the nlopt data structure
86  libmesh_assert(sys.rhs->size() == n);
87 
88  std::vector<double> grad;
89  sys.rhs->localize_to_one(grad);
90  for (unsigned i=0; i<n; ++i)
91  gradient[i] = grad[i];
92  }
93  else
94  libmesh_error_msg("Gradient function not defined in __libmesh_nlopt_objective");
95  }
96 
97  // Increment the iteration count.
98  solver->get_iteration_count()++;
99 
100  // Possibly print the current value of the objective function
101  if (solver->verbose)
102  libMesh::out << objective << std::endl;
103 
104  return objective;
105 }
106 
107 
108 
110  double * result,
111  unsigned n,
112  const double * x,
113  double * gradient,
114  void * data)
115 {
116  LOG_SCOPE("equality_constraints()", "NloptOptimizationSolver");
117 
118  libmesh_assert(data);
119 
120  // data should be a pointer to the solver (it was passed in as void *)
122  static_cast<NloptOptimizationSolver<Number> *> (data);
123 
124  OptimizationSystem & sys = solver->system();
125 
126  // We'll use current_local_solution below, so let's ensure that it's consistent
127  // with the vector x that was passed in.
128  if (sys.solution->size() != n)
129  libmesh_error_msg("Error: Input vector x has different length than sys.solution!");
130 
131  for (unsigned int i=sys.solution->first_local_index(); i<sys.solution->last_local_index(); i++)
132  sys.solution->set(i, x[i]);
133  sys.solution->close();
134 
135  // Impose constraints on the solution vector
136  sys.get_dof_map().enforce_constraints_exactly(sys);
137 
138  // Update sys.current_local_solution based on the solution vector
139  sys.update();
140 
141  // Call the user's equality constraints function if there is one.
143  if (eco)
144  {
145  eco->equality_constraints(*sys.current_local_solution,
146  *sys.C_eq,
147  sys);
148 
149  sys.C_eq->close();
150 
151  // Copy the values out of eq_constraints into 'result'.
152  // TODO: Even better would be if we could use 'result' directly
153  // as the storage of eq_constraints. Perhaps a serial-only
154  // NumericVector variant which supports this option?
155  for (unsigned i=0; i<m; ++i)
156  result[i] = (*sys.C_eq)(i);
157 
158  // If gradient != NULL, then the Jacobian matrix of the equality
159  // constraints has been requested. The incoming 'gradient'
160  // array is of length m*n and d(c_i)/d(x_j) = gradient[n*i+j].
161  if (gradient)
162  {
165 
166  if (eco_jac)
167  {
168  eco_jac->equality_constraints_jacobian(*sys.current_local_solution,
169  *sys.C_eq_jac,
170  sys);
171 
172  sys.C_eq_jac->close();
173 
174  // copy the Jacobian data to the gradient array
175  for (numeric_index_type i=0; i<m; i++)
176  {
177  std::set<numeric_index_type>::iterator it = sys.eq_constraint_jac_sparsity[i].begin();
178  std::set<numeric_index_type>::iterator it_end = sys.eq_constraint_jac_sparsity[i].end();
179  for ( ; it != it_end; ++it)
180  {
181  numeric_index_type dof_index = *it;
182  gradient[n*i+dof_index] = (*sys.C_eq_jac)(i,dof_index);
183  }
184  }
185  }
186  else
187  libmesh_error_msg("Jacobian function not defined in __libmesh_nlopt_equality_constraints");
188  }
189 
190  }
191  else
192  libmesh_error_msg("Constraints function not defined in __libmesh_nlopt_equality_constraints");
193 }
194 
195 
197  double * result,
198  unsigned n,
199  const double * x,
200  double * gradient,
201  void * data)
202 {
203  LOG_SCOPE("inequality_constraints()", "NloptOptimizationSolver");
204 
205  libmesh_assert(data);
206 
207  // data should be a pointer to the solver (it was passed in as void *)
209  static_cast<NloptOptimizationSolver<Number> *> (data);
210 
211  OptimizationSystem & sys = solver->system();
212 
213  // We'll use current_local_solution below, so let's ensure that it's consistent
214  // with the vector x that was passed in.
215  if (sys.solution->size() != n)
216  libmesh_error_msg("Error: Input vector x has different length than sys.solution!");
217 
218  for (unsigned int i=sys.solution->first_local_index(); i<sys.solution->last_local_index(); i++)
219  sys.solution->set(i, x[i]);
220  sys.solution->close();
221 
222  // Impose constraints on the solution vector
223  sys.get_dof_map().enforce_constraints_exactly(sys);
224 
225  // Update sys.current_local_solution based on the solution vector
226  sys.update();
227 
228  // Call the user's inequality constraints function if there is one.
230  if (ineco)
231  {
232  ineco->inequality_constraints(*sys.current_local_solution,
233  *sys.C_ineq,
234  sys);
235 
236  sys.C_ineq->close();
237 
238  // Copy the values out of ineq_constraints into 'result'.
239  // TODO: Even better would be if we could use 'result' directly
240  // as the storage of ineq_constraints. Perhaps a serial-only
241  // NumericVector variant which supports this option?
242  for (unsigned i=0; i<m; ++i)
243  result[i] = (*sys.C_ineq)(i);
244 
245  // If gradient != NULL, then the Jacobian matrix of the equality
246  // constraints has been requested. The incoming 'gradient'
247  // array is of length m*n and d(c_i)/d(x_j) = gradient[n*i+j].
248  if (gradient)
249  {
252 
253  if (ineco_jac)
254  {
255  ineco_jac->inequality_constraints_jacobian(*sys.current_local_solution,
256  *sys.C_ineq_jac,
257  sys);
258 
259  sys.C_ineq_jac->close();
260 
261  // copy the Jacobian data to the gradient array
262  for (numeric_index_type i=0; i<m; i++)
263  {
264  std::set<numeric_index_type>::iterator it = sys.ineq_constraint_jac_sparsity[i].begin();
265  std::set<numeric_index_type>::iterator it_end = sys.ineq_constraint_jac_sparsity[i].end();
266  for ( ; it != it_end; ++it)
267  {
268  numeric_index_type dof_index = *it;
269  gradient[n*i+dof_index] = (*sys.C_ineq_jac)(i,dof_index);
270  }
271  }
272  }
273  else
274  libmesh_error_msg("Jacobian function not defined in __libmesh_nlopt_inequality_constraints");
275  }
276 
277  }
278  else
279  libmesh_error_msg("Constraints function not defined in __libmesh_nlopt_inequality_constraints");
280 }
281 
282 //---------------------------------------------------------------------
283 
284 
285 
286 //---------------------------------------------------------------------
287 // NloptOptimizationSolver<> methods
288 template <typename T>
290  :
291  OptimizationSolver<T>(system_in),
292  _opt(libmesh_nullptr),
293  _result(NLOPT_SUCCESS),
294  _iteration_count(0),
295  _constraints_tolerance(1.e-8)
296 {
297 }
298 
299 
300 
301 template <typename T>
303 {
304  this->clear ();
305 }
306 
307 
308 
309 template <typename T>
311 {
312  if (this->initialized())
313  {
314  this->_is_initialized = false;
315 
316  nlopt_destroy(_opt);
317  }
318 }
319 
320 
321 
322 template <typename T>
324 {
325  // Initialize the data structures if not done so already.
326  if (!this->initialized())
327  {
328  this->_is_initialized = true;
329 
330  // By default, use the LD_SLSQP solver
331  std::string nlopt_algorithm_name = "LD_SLSQP";
332 
333  if (libMesh::on_command_line("--nlopt-algorithm"))
334  nlopt_algorithm_name = libMesh::command_line_next ("--nlopt-algorithm",
335  nlopt_algorithm_name);
336 
337  // Convert string to an nlopt algorithm type
338  std::map<std::string, nlopt_algorithm>::iterator it =
339  _nlopt_algorithms.find(nlopt_algorithm_name);
340 
341  if (it == _nlopt_algorithms.end())
342  libmesh_error_msg("Invalid nlopt algorithm requested on command line: " \
343  << nlopt_algorithm_name);
344 
345  _opt = nlopt_create(it->second, this->system().solution->size());
346  }
347 }
348 
349 
350 
351 template <typename T>
353 {
354  LOG_SCOPE("solve()", "NloptOptimizationSolver");
355 
356  this->init ();
357 
358  unsigned int nlopt_size = this->system().solution->size();
359 
360  // We have to have an objective function
362 
363  // Set routine for objective and (optionally) gradient evaluation
364  {
365  nlopt_result ierr =
366  nlopt_set_min_objective(_opt,
368  this);
369  if (ierr < 0)
370  libmesh_error_msg("NLopt failed to set min objective: " << ierr);
371  }
372 
374  {
375  // Need to actually compute the bounds vectors first
377 
378  std::vector<Real> nlopt_lb(nlopt_size);
379  std::vector<Real> nlopt_ub(nlopt_size);
380  for (unsigned int i=0; i<nlopt_size; i++)
381  {
382  nlopt_lb[i] = this->system().get_vector("lower_bounds")(i);
383  nlopt_ub[i] = this->system().get_vector("upper_bounds")(i);
384  }
385 
386  nlopt_set_lower_bounds(_opt, &nlopt_lb[0]);
387  nlopt_set_upper_bounds(_opt, &nlopt_ub[0]);
388  }
389 
390  // If we have an equality constraints object, tell NLopt about it.
391  if (this->equality_constraints_object)
392  {
393  // NLopt requires a vector to specify the tolerance for each constraint.
394  // NLopt makes a copy of this vector internally, so it's safe for us to
395  // let it go out of scope.
396  std::vector<double> equality_constraints_tolerances(this->system().C_eq->size(),
398 
399  // It would be nice to call the C interface directly, at least it should return an error
400  // code we could parse... unfortunately, there does not seem to be a way to extract
401  // the underlying nlopt_opt object from the nlopt::opt class!
402  nlopt_result ierr =
403  nlopt_add_equality_mconstraint(_opt,
404  equality_constraints_tolerances.size(),
406  this,
407  &equality_constraints_tolerances[0]);
408 
409  if (ierr < 0)
410  libmesh_error_msg("NLopt failed to add equality constraint: " << ierr);
411  }
412 
413  // If we have an inequality constraints object, tell NLopt about it.
415  {
416  // NLopt requires a vector to specify the tolerance for each constraint
417  std::vector<double> inequality_constraints_tolerances(this->system().C_ineq->size(),
419 
420  nlopt_add_inequality_mconstraint(_opt,
421  inequality_constraints_tolerances.size(),
423  this,
424  &inequality_constraints_tolerances[0]);
425  }
426 
427  // Set a relative tolerance on the optimization parameters
428  nlopt_set_ftol_rel(_opt, this->objective_function_relative_tolerance);
429 
430  // Set the maximum number of allowed objective function evaluations
431  nlopt_set_maxeval(_opt, this->max_objective_function_evaluations);
432 
433  // Reset internal iteration counter
434  this->_iteration_count = 0;
435 
436  // Perform the optimization
437  std::vector<Real> x(nlopt_size);
438  Real min_val = 0.;
439  _result = nlopt_optimize(_opt, &x[0], &min_val);
440 
441  if (_result < 0)
442  libMesh::out << "NLopt failed!" << std::endl;
443  else
444  libMesh::out << "NLopt obtained optimal value: "
445  << min_val
446  << " in "
447  << this->get_iteration_count()
448  << " iterations."
449  << std::endl;
450 }
451 
452 
453 template <typename T>
455 {
456  libMesh::out << "NLopt optimization solver convergence/divergence reason: "
457  << this->get_converged_reason() << std::endl;
458 }
459 
460 
461 
462 template <typename T>
464 {
465  return static_cast<int>(_result);
466 }
467 
468 
469 //------------------------------------------------------------------
470 // Explicit instantiations
471 template class NloptOptimizationSolver<Number>;
472 
473 } // namespace libMesh
474 
475 
476 
477 #endif // #if defined(LIBMESH_HAVE_NLOPT) && !defined(LIBMESH_USE_COMPLEX_NUMBERS)
Abstract base class to be used to calculate the inequality constraints.
OptimizationSystem::ComputeObjective * objective_object
Object that computes the objective function f(X) at the input iterate X.
double _constraints_tolerance
NLopt requires us to specify a tolerance for the constraints.
virtual void gradient(const NumericVector< Number > &X, NumericVector< Number > &grad_f, sys_type &S)=0
This function will be called to compute the gradient of the objective function, and must be implement...
nlopt_opt _opt
Optimization solver context.
This class provides an interface to the NLopt optimization solvers.
nlopt_result _result
Store the result (i.e.
friend void __libmesh_nlopt_equality_constraints(unsigned m, double *result, unsigned n, const double *x, double *gradient, void *data)
ImplicitSystem & sys
Real objective_function_relative_tolerance
Required change in objective function which signals convergence.
friend void __libmesh_nlopt_inequality_constraints(unsigned m, double *result, unsigned n, const double *x, double *gradient, void *data)
T command_line_next(const std::string &, T)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:964
virtual void inequality_constraints_jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &C_ineq_jac, sys_type &S)=0
This function will be called to evaluate the Jacobian of C_ineq(X).
const class libmesh_nullptr_t libmesh_nullptr
OptimizationSystem::ComputeGradient * gradient_object
Object that computes the gradient grad_f(X) of the objective function at the input iterate X...
The libMesh namespace provides an interface to certain functionality in the library.
unsigned _iteration_count
Stores the current iteration index (incremented at each call of __libmesh_nlopt_objective).
bool _is_initialized
Flag indicating if the data structures have been initialized.
virtual void print_converged_reason() libmesh_override
Prints a useful message about why the latest optimization solve con(di)verged.
const sys_type & system() const
libmesh_assert(j)
PetscDiffSolver & solver
OptimizationSystem::ComputeInequalityConstraints * inequality_constraints_object
Object that computes the inequality constraints vector C_ineq(X).
friend double __libmesh_nlopt_objective(unsigned n, const double *x, double *gradient, void *data)
PetscErrorCode Vec x
OptimizationSystem::ComputeEqualityConstraintsJacobian * equality_constraints_jacobian_object
Object that computes the Jacobian of C_eq(X).
Abstract base class to be used to calculate the equality constraints.
dof_id_type numeric_index_type
Definition: id_types.h:92
This System subclass enables us to assemble an objective function, gradient, Hessian and bounds for o...
NloptOptimizationSolver(sys_type &system)
Constructor.
virtual void init() libmesh_override
Initialize data structures if not done so already.
unsigned int max_objective_function_evaluations
Maximum number of objective function evaluations allowed.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
virtual void equality_constraints(const NumericVector< Number > &X, NumericVector< Number > &C_eq, sys_type &S)=0
This function will be called to evaluate the equality constraints vector C_eq(X). ...
virtual void solve() libmesh_override
Call the NLopt solver.
void __libmesh_nlopt_inequality_constraints(unsigned m, double *result, unsigned n, const double *x, double *gradient, void *data)
virtual void lower_and_upper_bounds(sys_type &S)=0
This function should update the following two vectors: this->get_vector("lower_bounds"), this->get_vector("upper_bounds").
virtual int get_converged_reason() libmesh_override
Abstract base class to be used to calculate the Jacobian of the equality constraints.
virtual Number objective(const NumericVector< Number > &X, sys_type &S)=0
This function will be called to compute the objective function to be minimized, and must be implement...
virtual void equality_constraints_jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &C_eq_jac, sys_type &S)=0
This function will be called to evaluate the Jacobian of C_eq(X).
OptimizationSystem::ComputeLowerAndUpperBounds * lower_and_upper_bounds_object
Object that computes the lower and upper bounds vectors.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool verbose
Control how much is output from the OptimizationSolver as it&#39;s running.
This base class can be inherited from to provide interfaces to optimization solvers from different pa...
PetscErrorCode ierr
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
OStreamProxy out
double __libmesh_nlopt_objective(unsigned n, const double *x, double *gradient, void *data)
IterBase * data
Ideally this private member data should have protected access.
virtual void clear() libmesh_override
Release all memory and clear data structures.
virtual void inequality_constraints(const NumericVector< Number > &X, NumericVector< Number > &C_ineq, sys_type &S)=0
This function will be called to evaluate the equality constraints vector C_ineq(X).
OptimizationSystem::ComputeInequalityConstraintsJacobian * inequality_constraints_jacobian_object
Object that computes the Jacobian of C_ineq(X).
static std::map< std::string, nlopt_algorithm > _nlopt_algorithms
OptimizationSystem::ComputeEqualityConstraints * equality_constraints_object
Object that computes the equality constraints vector C_eq(X).
Abstract base class to be used to calculate the Jacobian of the inequality constraints.
void __libmesh_nlopt_equality_constraints(unsigned m, double *result, unsigned n, const double *x, double *gradient, void *data)