Scippy

SCIP

Solving Constraint Integer Programs

nlpi_ipopt.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file nlpi_ipopt.cpp
17  * @ingroup DEFPLUGINS_NLPI
18  * @brief Ipopt NLP interface
19  * @author Stefan Vigerske
20  * @author Benjamin Müller
21  *
22  * @todo if too few degrees of freedom, solve a slack-minimization problem instead?
23  * @todo automatically switch to Hessian approximation if Hessian is dense or slow? (only do so if opttol/solvertol is large?)
24  *
25  * This file can only be compiled if Ipopt is available.
26  * Otherwise, to resolve public functions, use nlpi_ipopt_dummy.c.
27  * Since the dummy code is C instead of C++, it has been moved into a separate file.
28  */
29 
30 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
31 
32 #include "scip/nlpi_ipopt.h"
33 
34 #include "scip/nlpioracle.h"
35 #include "scip/exprinterpret.h"
36 #include "scip/scip_nlpi.h"
37 #include "scip/scip_randnumgen.h"
38 #include "scip/scip_mem.h"
39 #include "scip/scip_message.h"
40 #include "scip/scip_general.h"
41 #include "scip/scip_numerics.h"
42 #include "scip/scip_param.h"
43 #include "scip/scip_solve.h"
44 #include "scip/scip_copy.h"
45 #include "scip/pub_misc.h"
46 #include "scip/pub_paramset.h"
47 
48 #include <new> /* for std::bad_alloc */
49 #include <sstream>
50 #include <cstring>
51 
52 /* turn off some lint warnings for file */
53 /*lint --e{1540,750,3701}*/
54 
55 #include "IpoptConfig.h"
56 
57 #if defined(__GNUC__) && IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
58 #pragma GCC diagnostic ignored "-Wshadow"
59 #endif
60 #include "IpIpoptApplication.hpp"
61 #include "IpIpoptCalculatedQuantities.hpp"
62 #include "IpSolveStatistics.hpp"
63 #include "IpJournalist.hpp"
64 #include "IpIpoptData.hpp"
65 #include "IpTNLPAdapter.hpp"
66 #include "IpOrigIpoptNLP.hpp"
67 #include "IpLapack.hpp"
68 #if defined(__GNUC__) && IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
69 #pragma GCC diagnostic warning "-Wshadow"
70 #endif
71 
72 #if IPOPT_VERSION_MAJOR < 3 || (IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 12)
73 #error "The Ipopt interface requires at least 3.12.0"
74 #endif
75 
76 /* MUMPS that can be used by Ipopt is not threadsafe
77  * If we want SCIP to be threadsafe (SCIP_THREADSAFE), have std::mutex (C++11 or higher), and use Ipopt before 3.14,
78  * then we protect the call to Ipopt by a mutex if MUMPS is used as linear solver.
79  * Thus, we allow only one Ipopt run at a time.
80  * Ipopt 3.14 has this build-in to its MUMPS interface, so we won't have to take care of this.
81  */
82 #if defined(SCIP_THREADSAFE) && __cplusplus >= 201103L && IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
83 #define PROTECT_SOLVE_BY_MUTEX
84 #include <mutex>
85 static std::mutex solve_mutex; /*lint !e1756*/
86 #endif
87 
88 using namespace Ipopt;
89 
90 #define NLPI_NAME "ipopt" /**< short concise name of solver */
91 #define NLPI_DESC "Ipopt interface" /**< description of solver */
92 #define NLPI_PRIORITY 1000 /**< priority */
93 
94 #define MAXPERTURB 0.01 /**< maximal perturbation of bounds in starting point heuristic */
95 #define FEASTOLFACTOR 0.9 /**< factor for user-given feasibility tolerance to get feasibility tolerance that is actually passed to Ipopt */
96 
97 #define DEFAULT_RANDSEED 71 /**< initial random seed */
98 
99 // enable this to collect statistics on number of iterations and problem characteristics in csv-form in log
100 // note that this overwrites given iterlimit
101 // see https://git.zib.de/integer/scip/-/snippets/1213 for some script that evaluates the collected data
102 // #define COLLECT_SOLVESTATS
103 
104 /* Convergence check (see ScipNLP::intermediate_callback)
105  *
106  * If the fastfail option is set to aggressive, then we stop Ipopt if the reduction in
107  * primal infeasibility is not sufficient for a consecutive number of iterations.
108  * With the parameters as given below, we require Ipopt to
109  * - not increase the primal infeasibility after 5 iterations
110  * - reduce the primal infeasibility by at least 50% within 10 iterations
111  * - reduce the primal infeasibility by at least 90% within 30 iterations
112  * The targets are updated once they are reached and the limit on allowed iterations to reach the new target is reset.
113  *
114  * In certain situations, it is allowed to exceed an iteration limit:
115  * - If we are in the first 10 (convcheck_startiter) iterations.
116  * - If we are within 10 (convcheck_startiter) iterations after the restoration phase ended.
117  * The reason for this is that during feasibility restoration phase Ipopt aims completely on
118  * reducing constraint violation, completely forgetting the objective function.
119  * When returning from feasibility restoration and considering the original objective again,
120  * it is unlikely that Ipopt will continue to decrease primal infeasibility, since it may now target on
121  * more on optimality again. Thus, we do not check convergence for a number of iterations.
122  * - If the target on dual infeasibility reduction has been achieved, we are below twice the iteration limit, and
123  * we are not in restoration mode.
124  * The reason for this is that if Ipopt makes good progress towards optimality,
125  * we want to allow some more iterations where primal infeasibility is not reduced.
126  * However, in restoration mode, dual infeasibility does not correspond to the original problem and
127  * the complete aim is to restore primal infeasibility.
128  */
129 static const int convcheck_nchecks = 3; /**< number of convergence checks */
130 static const int convcheck_startiter = 10; /**< iteration where to start convergence checking */
131 static const int convcheck_maxiter[convcheck_nchecks] = { 5, 15, 30 }; /**< maximal number of iterations to achieve each convergence check */
132 static const SCIP_Real convcheck_minred[convcheck_nchecks] = { 1.0, 0.5, 0.1 }; /**< minimal required infeasibility reduction in each convergence check */
133 
134 /// integer parameters of Ipopt to make available via SCIP parameters
135 static const char* ipopt_int_params[] =
136  { "print_level" }; // print_level must be first
137 
138 /// string parameters of Ipopt to make available via SCIP parameters
139 static const char* ipopt_string_params[] =
140  { "linear_solver",
141  "hsllib",
142  "pardisolib",
143  "linear_system_scaling",
144  "nlp_scaling_method",
145  "mu_strategy",
146  "hessian_approximation"
147  };
148 
149 class ScipNLP;
150 
151 struct SCIP_NlpiData
152 {
153 public:
154  char* optfile; /**< Ipopt options file to read */
155  int print_level; /**< print_level set via nlpi/ipopt/print_level option */
156  SCIP_Real warm_start_push; /**< value to use for Ipopt's warm_start_bound_push/frac options */
157 
158  /** constructor */
159  explicit SCIP_NlpiData()
160  : optfile(NULL), print_level(-1), warm_start_push(1e-9)
161  { }
162 };
163 
164 struct SCIP_NlpiProblem
165 {
166 public:
167  SCIP_NLPIORACLE* oracle; /**< Oracle-helper to store and evaluate NLP */
168  SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
169 
170  SmartPtr<IpoptApplication> ipopt; /**< Ipopt application */
171  SmartPtr<ScipNLP> nlp; /**< NLP in Ipopt form */
172 
173  bool firstrun; /**< whether the next NLP solve will be the first one */
174  bool samestructure;/**< whether the NLP solved next will still have the same (Ipopt-internal) structure (same number of variables, constraints, bounds, and nonzero pattern) */
175 
176  SCIP_NLPSOLSTAT solstat; /**< status of current solution (if any) */
177  SCIP_NLPTERMSTAT termstat; /**< termination status of last solve (if any) */
178  bool solprimalvalid;/**< whether primal solution values are available (solprimals has meaningful values) */
179  bool solprimalgiven;/**< whether primal solution values were set by caller */
180  bool soldualvalid; /**< whether dual solution values are available (soldual* have meaningful values) */
181  bool soldualgiven; /**< whether dual solution values were set by caller */
182  SCIP_Real* solprimals; /**< primal solution values, if available */
183  SCIP_Real* soldualcons; /**< dual solution values of constraints, if available */
184  SCIP_Real* soldualvarlb; /**< dual solution values of variable lower bounds, if available */
185  SCIP_Real* soldualvarub; /**< dual solution values of variable upper bounds, if available */
186  SCIP_Real solobjval; /**< objective function value in solution from last run */
187  SCIP_Real solconsviol; /**< constraint violation of primal solution, if available */
188  SCIP_Real solboundviol; /**< variable bound violation of primal solution, if available */
189  int lastniter; /**< number of iterations in last run */
190  SCIP_Real lasttime; /**< time spend in last run */
191 
192  /** constructor */
194  : oracle(NULL), randnumgen(NULL),
195  firstrun(true), samestructure(true),
197  solprimalvalid(false), solprimalgiven(false), soldualvalid(false), soldualgiven(false),
198  solprimals(NULL), soldualcons(NULL), soldualvarlb(NULL), soldualvarub(NULL),
199  solobjval(SCIP_INVALID), solconsviol(SCIP_INVALID), solboundviol(SCIP_INVALID),
200  lastniter(-1), lasttime(-1.0)
201  { }
202 };
203 
204 /** TNLP implementation for SCIPs NLP */
205 class ScipNLP : public TNLP
206 {
207 private:
208  SCIP_NLPIPROBLEM* nlpiproblem; /**< NLPI problem data */
209  SCIP* scip; /**< SCIP data structure */
210  SCIP_NLPPARAM param; /**< NLP solve parameters */
211 
212  SCIP_Real conv_prtarget[convcheck_nchecks]; /**< target primal infeasibility for each convergence check */
213  SCIP_Real conv_dutarget[convcheck_nchecks]; /**< target dual infeasibility for each convergence check */
214  int conv_iterlim[convcheck_nchecks]; /**< iteration number where target primal infeasibility should to be achieved */
215  int conv_lastrestoiter; /**< last iteration number in restoration mode, or -1 if none */
216 
217  unsigned int current_x; /**< unique number that identifies current iterate (x): incremented when Ipopt calls with new_x=true */
218  unsigned int last_f_eval_x; /**< the number of the iterate for which the objective was last evaluated (eval_f) */
219  unsigned int last_g_eval_x; /**< the number of the iterate for which the constraints were last evaluated (eval_g) */
220 
221 public:
222  bool approxhessian; /**< do we tell Ipopt to approximate the hessian? (may also be false if user set to approx. hessian via option file) */
223 
224  // cppcheck-suppress uninitMemberVar
225  /** constructor */
226  ScipNLP(
227  SCIP_NLPIPROBLEM* nlpiproblem_ = NULL,/**< NLPI problem data */
228  SCIP* scip_ = NULL /**< SCIP data structure */
229  )
230  : nlpiproblem(nlpiproblem_), scip(scip_),
231  conv_lastrestoiter(-1),
232  current_x(1), last_f_eval_x(0), last_g_eval_x(0),
233  approxhessian(false)
234  {
235  assert(scip != NULL);
236  }
237 
238  /** destructor */
239  ~ScipNLP()
240  { /*lint --e{1540}*/
241  }
242 
243  /** initialize for new solve */
244  void initializeSolve(
245  SCIP_NLPIPROBLEM* nlpiproblem_, /**< NLPI problem */
246  const SCIP_NLPPARAM& nlpparam /**< NLP solve parameters */
247  )
248  {
249  assert(nlpiproblem_ != NULL);
250  nlpiproblem = nlpiproblem_;
251  param = nlpparam;
252 
253  // since we are about to start a new solve, use this opportunity to reset the counts on x
254  current_x = 1;
255  last_f_eval_x = 0;
256  last_g_eval_x = 0;
257  }
258 
259  /** Method to return some info about the nlp */
260  bool get_nlp_info(
261  Index& n, /**< place to store number of variables */
262  Index& m, /**< place to store number of constraints */
263  Index& nnz_jac_g, /**< place to store number of nonzeros in jacobian */
264  Index& nnz_h_lag, /**< place to store number of nonzeros in hessian */
265  IndexStyleEnum& index_style /**< place to store used index style (0-based or 1-based) */
266  );
267 
268  /** Method to return the bounds for my problem */
269  bool get_bounds_info(
270  Index n, /**< number of variables */
271  Number* x_l, /**< buffer to store lower bounds on variables */
272  Number* x_u, /**< buffer to store upper bounds on variables */
273  Index m, /**< number of constraints */
274  Number* g_l, /**< buffer to store lower bounds on constraints */
275  Number* g_u /**< buffer to store lower bounds on constraints */
276  );
277 
278  /** Method to return the starting point for the algorithm */
279  bool get_starting_point(
280  Index n, /**< number of variables */
281  bool init_x, /**< whether initial values for primal values are requested */
282  Number* x, /**< buffer to store initial primal values */
283  bool init_z, /**< whether initial values for dual values of variable bounds are requested */
284  Number* z_L, /**< buffer to store dual values for variable lower bounds */
285  Number* z_U, /**< buffer to store dual values for variable upper bounds */
286  Index m, /**< number of constraints */
287  bool init_lambda, /**< whether initial values for dual values of constraints are required */
288  Number* lambda /**< buffer to store dual values of constraints */
289  );
290 
291  /** Method to return the number of nonlinear variables. */
292  Index get_number_of_nonlinear_variables();
293 
294  /** Method to return the indices of the nonlinear variables */
295  bool get_list_of_nonlinear_variables(
296  Index num_nonlin_vars, /**< number of nonlinear variables */
297  Index* pos_nonlin_vars /**< array to fill with indices of nonlinear variables */
298  );
299 
300  /** Method to return metadata about variables and constraints */
301  bool get_var_con_metadata(
302  Index n, /**< number of variables */
303  StringMetaDataMapType& var_string_md, /**< variable meta data of string type */
304  IntegerMetaDataMapType& var_integer_md,/**< variable meta data of integer type */
305  NumericMetaDataMapType& var_numeric_md,/**< variable meta data of numeric type */
306  Index m, /**< number of constraints */
307  StringMetaDataMapType& con_string_md, /**< constraint meta data of string type */
308  IntegerMetaDataMapType& con_integer_md,/**< constraint meta data of integer type */
309  NumericMetaDataMapType& con_numeric_md /**< constraint meta data of numeric type */
310  );
311 
312  /** Method to return the objective value */
313  bool eval_f(
314  Index n, /**< number of variables */
315  const Number* x, /**< point to evaluate */
316  bool new_x, /**< whether some function evaluation method has been called for this point before */
317  Number& obj_value /**< place to store objective function value */
318  );
319 
320  /** Method to return the gradient of the objective */
321  bool eval_grad_f(
322  Index n, /**< number of variables */
323  const Number* x, /**< point to evaluate */
324  bool new_x, /**< whether some function evaluation method has been called for this point before */
325  Number* grad_f /**< buffer to store objective gradient */
326  );
327 
328  /** Method to return the constraint residuals */
329  bool eval_g(
330  Index n, /**< number of variables */
331  const Number* x, /**< point to evaluate */
332  bool new_x, /**< whether some function evaluation method has been called for this point before */
333  Index m, /**< number of constraints */
334  Number* g /**< buffer to store constraint function values */
335  );
336 
337  /** Method to return:
338  * 1) The structure of the jacobian (if "values" is NULL)
339  * 2) The values of the jacobian (if "values" is not NULL)
340  */
341  bool eval_jac_g(
342  Index n, /**< number of variables */
343  const Number* x, /**< point to evaluate */
344  bool new_x, /**< whether some function evaluation method has been called for this point before */
345  Index m, /**< number of constraints */
346  Index nele_jac, /**< number of nonzero entries in jacobian */
347  Index* iRow, /**< buffer to store row indices of nonzero jacobian entries, or NULL if values
348  * are requested */
349  Index* jCol, /**< buffer to store column indices of nonzero jacobian entries, or NULL if values
350  * are requested */
351  Number* values /**< buffer to store values of nonzero jacobian entries, or NULL if structure is
352  * requested */
353  );
354 
355  /** Method to return:
356  * 1) The structure of the hessian of the lagrangian (if "values" is NULL)
357  * 2) The values of the hessian of the lagrangian (if "values" is not NULL)
358  */
359  bool eval_h(
360  Index n, /**< number of variables */
361  const Number* x, /**< point to evaluate */
362  bool new_x, /**< whether some function evaluation method has been called for this point before */
363  Number obj_factor, /**< weight for objective function */
364  Index m, /**< number of constraints */
365  const Number* lambda, /**< weights for constraint functions */
366  bool new_lambda, /**< whether the hessian has been evaluated for these values of lambda before */
367  Index nele_hess, /**< number of nonzero entries in hessian */
368  Index* iRow, /**< buffer to store row indices of nonzero hessian entries, or NULL if values
369  * are requested */
370  Index* jCol, /**< buffer to store column indices of nonzero hessian entries, or NULL if values
371  * are requested */
372  Number* values /**< buffer to store values of nonzero hessian entries, or NULL if structure is requested */
373  );
374 
375  /** Method called by the solver at each iteration.
376  *
377  * Checks whether Ctrl-C was hit.
378  */
379  bool intermediate_callback(
380  AlgorithmMode mode, /**< current mode of algorithm */
381  Index iter, /**< current iteration number */
382  Number obj_value, /**< current objective value */
383  Number inf_pr, /**< current primal infeasibility */
384  Number inf_du, /**< current dual infeasibility */
385  Number mu, /**< current barrier parameter */
386  Number d_norm, /**< current gradient norm */
387  Number regularization_size,/**< current size of regularization */
388  Number alpha_du, /**< current dual alpha */
389  Number alpha_pr, /**< current primal alpha */
390  Index ls_trials, /**< current number of linesearch trials */
391  const IpoptData* ip_data, /**< pointer to Ipopt Data */
392  IpoptCalculatedQuantities* ip_cq /**< pointer to current calculated quantities */
393  );
394 
395  /** This method is called when the algorithm is complete so the TNLP can store/write the solution. */
396  void finalize_solution(
397  SolverReturn status, /**< solve and solution status */
398  Index n, /**< number of variables */
399  const Number* x, /**< primal solution values */
400  const Number* z_L, /**< dual values of variable lower bounds */
401  const Number* z_U, /**< dual values of variable upper bounds */
402  Index m, /**< number of constraints */
403  const Number* g, /**< values of constraints */
404  const Number* lambda, /**< dual values of constraints */
405  Number obj_value, /**< objective function value */
406  const IpoptData* data, /**< pointer to Ipopt Data */
407  IpoptCalculatedQuantities* cq /**< pointer to calculated quantities */
408  );
409 };
410 
411 /** A particular Ipopt::Journal implementation that uses the SCIP message routines for output. */
412 class ScipJournal : public Ipopt::Journal {
413 private:
414  SCIP* scip; /**< SCIP data structure */
415 
416 public:
417  ScipJournal(
418  const char* name, /**< name of journal */
419  Ipopt::EJournalLevel default_level, /**< default verbosity level */
420  SCIP* scip_ /**< SCIP data structure */
421  )
422  : Ipopt::Journal(name, default_level),
423  scip(scip_)
424  { }
425 
426  ~ScipJournal() { }
427 
428 protected:
429  /*lint -e{715}*/
430  void PrintImpl(
431  Ipopt::EJournalCategory category, /**< category of message */
432  Ipopt::EJournalLevel level, /**< verbosity level of message */
433  const char* str /**< message to print */
434  )
435  { /*lint --e{715} */
436  if( level == J_ERROR )
437  {
438  SCIPerrorMessage("%s", str);
439  }
440  else
441  {
442  SCIPinfoMessage(scip, NULL, "%s", str);
443  }
444  }
445 
446  /*lint -e{715}*/
447  void PrintfImpl(
448  Ipopt::EJournalCategory category, /**< category of message */
449  Ipopt::EJournalLevel level, /**< verbosity level of message */
450  const char* pformat, /**< message printing format */
451  va_list ap /**< arguments of message */
452  )
453  { /*lint --e{715} */
454  if( level == J_ERROR )
455  {
456  SCIPmessagePrintErrorHeader(__FILE__, __LINE__);
457  SCIPmessageVPrintError(pformat, ap);
458  }
459  else
460  {
461  SCIPmessageVPrintInfo(SCIPgetMessagehdlr(scip), pformat, ap);
462  }
463  }
464 
465  void FlushBufferImpl() { }
466 };
467 
468 /** sets status codes to mark that last NLP solve is no longer valid (usually because the NLP changed) */
469 static
471  SCIP_NLPIPROBLEM* problem /**< data structure of problem */
472  )
473 {
474  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
475  problem->termstat = SCIP_NLPTERMSTAT_OTHER;
476  problem->solobjval = SCIP_INVALID;
477  problem->solconsviol = SCIP_INVALID;
478  problem->solboundviol = SCIP_INVALID;
479  problem->lastniter = -1;
480  problem->lasttime = -1.0;
481 }
482 
483 /** sets solution values to be invalid and calls invalidateSolved() */
484 static
486  SCIP_NLPIPROBLEM* problem /**< data structure of problem */
487  )
488 {
489  assert(problem != NULL);
490 
491  problem->solprimalvalid = false;
492  problem->solprimalgiven = false;
493  problem->soldualvalid = false;
494  problem->soldualgiven = false;
495 
496  invalidateSolved(problem);
497 }
498 
499 /** makes sure a starting point (initial guess) is available */
500 static
502  SCIP* scip, /**< SCIP data structure */
503  SCIP_NLPIPROBLEM* problem, /**< data structure of problem */
504  SCIP_Bool& warmstart /**< whether a warmstart has been requested */
505  )
506 {
507  SCIP_Real lb, ub;
508  int n;
509 
510  assert(problem != NULL);
511 
512  // disable warmstart if no primal or dual solution values are available
513  if( warmstart && (!problem->solprimalvalid || !problem->soldualvalid ))
514  {
515  SCIPdebugMsg(scip, "Disable warmstart as no primal or dual solution available.\n");
516  warmstart = false;
517  }
518 
519  // continue below with making up a random primal starting point if
520  // the user did not set a starting point and warmstart is disabled (so the last solution shouldn't be used)
521  // (if warmstart, then due to the checks above we must now have valid primal and dual solution values)
522  if( problem->solprimalgiven || warmstart )
523  {
524  // so we must have a primal solution to start from
525  // if warmstart, then we also need to have a dual solution to start from
526  // if warmstart and primal solution is given by user, then also dual solution should have been given by user
527  assert(problem->solprimalvalid);
528  assert(problem->solprimals != NULL);
529  assert(!warmstart || !problem->solprimalgiven || problem->soldualgiven);
530  assert(!warmstart || problem->soldualcons != NULL);
531  assert(!warmstart || problem->soldualvarlb != NULL);
532  assert(!warmstart || problem->soldualvarub != NULL);
533  SCIPdebugMsg(scip, "Starting solution for %sstart available from %s.\n",
534  warmstart ? "warm" : "cold",
535  problem->solprimalgiven ? "user" : "previous solve");
536  return SCIP_OKAY;
537  }
538 
539  SCIPdebugMsg(scip, "Starting solution for coldstart not available. Making up something by projecting 0 onto variable bounds and adding a random perturbation.\n");
540 
541  n = SCIPnlpiOracleGetNVars(problem->oracle);
542 
543  if( problem->randnumgen == NULL )
544  {
546  }
547 
548  if( problem->solprimals == NULL )
549  {
550  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->solprimals, n) );
551  }
552 
553  for( int i = 0; i < n; ++i )
554  {
555  lb = SCIPnlpiOracleGetVarLbs(problem->oracle)[i];
556  ub = SCIPnlpiOracleGetVarUbs(problem->oracle)[i];
557  if( lb > 0.0 )
558  problem->solprimals[i] = SCIPrandomGetReal(problem->randnumgen, lb, lb + MAXPERTURB*MIN(1.0, ub-lb));
559  else if( ub < 0.0 )
560  problem->solprimals[i] = SCIPrandomGetReal(problem->randnumgen, ub - MAXPERTURB*MIN(1.0, ub-lb), ub);
561  else
562  problem->solprimals[i] = SCIPrandomGetReal(problem->randnumgen, MAX(lb, -MAXPERTURB*MIN(1.0, ub-lb)), MIN(ub, MAXPERTURB*MIN(1.0, ub-lb)));
563  }
564  problem->solprimalvalid = true;
565 
566  return SCIP_OKAY;
567 }
568 
569 /// pass NLP solve parameters to Ipopt
570 static
572  SCIP* scip, /**< SCIP data structure */
573  SCIP_NLPIDATA* nlpidata, /**< NLPI data */
574  SCIP_NLPIPROBLEM* nlpiproblem, /**< NLP */
575  const SCIP_NLPPARAM param /**< solve parameters */
576  )
577 {
578  assert(scip != NULL);
579  assert(nlpiproblem != NULL);
580 
581  if( nlpidata->print_level < 0 ) // if nlpi/ipopt/print_level param has not been set
582  {
583  switch( param.verblevel )
584  {
585  case 0:
586  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("print_level", J_ERROR);
587  break;
588  case 1:
589  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("print_level", J_SUMMARY);
590  break;
591  case 2:
592  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("print_level", J_ITERSUMMARY);
593  break;
594  case 3:
595  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("print_level", J_DETAILED);
596  break;
597  default:
598  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("print_level", MIN(J_ITERSUMMARY + (param.verblevel-1), J_ALL));
599  break;
600  }
601  }
602 
603 #ifdef SCIP_DISABLED_CODE
604  if( param.iterlimit < 0 )
605  {
606  if( nlpidata->autoiterlim > 0 )
607  {
608  param.iterlimit = nlpidata->autoiterlim;
609  }
610  else
611  {
612  int nvars = 0; // number of variables, including slacks
613  int nnlvars = 0; // number of nonlinear variables
614  int nvarbnd = 0; // number of finite lower and upper bounds, including slacks
615  int nlincons = 0; // number of linear constraints
616  int nnlcons = 0; // number of nonlinear constraints
617  int jacnnz = 0; // number of nonzeros in Jacobian
618 
619  int n = SCIPnlpiOracleGetNVars(nlpiproblem->oracle);
620  int m = SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle);
621  const SCIP_Real* lbs = SCIPnlpiOracleGetVarLbs(nlpiproblem->oracle);
622  const SCIP_Real* ubs = SCIPnlpiOracleGetVarUbs(nlpiproblem->oracle);
623 
624  for( int i = 0; i < n; ++i )
625  {
626  // skip fixed vars
627  if( SCIPisEQ(scip, lbs[i], ubs[i]) )
628  continue;
629 
630  ++nvars;
631 
632  if( SCIPnlpiOracleIsVarNonlinear(scip, nlpiproblem->oracle, i) )
633  ++nnlvars;
634 
635  // every variable lower bound contributes a ln(x-lb) term in the barrier problem
636  if( !SCIPisInfinity(scip, -lbs[i]) )
637  ++nvarbnd;
638 
639  // every variable upper bound contributes a ln(ub-x) term in the barrier problem
640  if( !SCIPisInfinity(scip, ubs[i]) )
641  ++nvarbnd;
642  }
643 
644  for( int i = 0; i < m; ++i )
645  {
646  if( SCIPnlpiOracleIsConstraintNonlinear(nlpiproblem->oracle, i) )
647  ++nnlcons;
648  else
649  ++nlincons;
650 
651  SCIP_Real lhs = SCIPnlpiOracleGetConstraintLhs(nlpiproblem->oracle, i);
652  SCIP_Real rhs = SCIPnlpiOracleGetConstraintRhs(nlpiproblem->oracle, i);
653 
654  // every inequality contributes a slack variable
655  if( !SCIPisEQ(scip, lhs, rhs) )
656  ++nvars;
657 
658  // every constraint lhs contributes a ln(slack-lhs) term in the barrier problem
659  if( !SCIPisInfinity(scip, -lhs) )
660  ++nvarbnd;
661  // every constraint rhs contributes a ln(rhs-slack) term in the barrier problem
662  if( !SCIPisInfinity(scip, rhs) )
663  ++nvarbnd;
664  }
665 
666  const int* offset;
667  SCIP_CALL( SCIPnlpiOracleGetJacobianSparsity(scip, nlpiproblem->oracle, &offset, NULL) );
668  jacnnz = offset[m];
669 
670  /* fitting data from NLP runs gave the following coefficients (see also !2634):
671  * intercept +40.2756726751
672  * jacnnz -0.0021259769
673  * jacnnz_sqrt +2.0121042012
674  * nlincons -0.0374801925
675  * nlincons_sqrt +2.9562232443
676  * nnlcons -0.0133039200
677  * nnlcons_sqrt -0.0412118434
678  * nnlvars -0.0702890379
679  * nnlvars_sqrt +7.0920920430
680  * nvarbnd +0.0183592749
681  * nvarbnd_sqrt -4.7218258847
682  * nvars +0.0112944627
683  * nvars_sqrt -0.8365873360
684  */
685  param.iterlimit = SCIPfloor(scip, 40.2756726751
686  -0.0021259769 * jacnnz +2.0121042012 * sqrt(jacnnz)
687  -0.0374801925 * nlincons +2.9562232443 * sqrt(nlincons)
688  -0.0133039200 * nnlcons -0.0412118434 * sqrt(nnlcons)
689  -0.0702890379 * nnlvars +7.0920920430 * sqrt(nnlvars)
690  +0.0183592749 * nvarbnd -4.7218258847 * sqrt(nvarbnd)
691  +0.0112944627 * nvars -0.8365873360 * sqrt(nvars));
692  SCIPdebugMsg(scip, "Iteration limit guess: %d\n", param.iterlimit);
693  /* but with all the negative coefficients, let's also ensure some minimal number of iterations */
694  if( param.iterlimit < 50 )
695  param.iterlimit = 50;
696  }
697  if( nlpidata->print_level >= J_SUMMARY || param.verblevel > 0 )
698  SCIPinfoMessage(scip, NULL, "Chosen iteration limit to be %d\n", param.iterlimit);
699  }
700 #endif
701  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("max_iter", param.iterlimit);
702 
703  (void) nlpiproblem->ipopt->Options()->SetNumericValue("constr_viol_tol", FEASTOLFACTOR * param.feastol);
704  (void) nlpiproblem->ipopt->Options()->SetNumericValue("acceptable_constr_viol_tol", FEASTOLFACTOR * param.feastol);
705 
706  /* set optimality tolerance parameters in Ipopt
707  *
708  * Sets dual_inf_tol and compl_inf_tol to opttol and tol to solvertol.
709  * We leave acceptable_dual_inf_tol and acceptable_compl_inf_tol untouched for now, which means that if Ipopt has convergence problems, then
710  * it can stop with a solution that is still feasible, but essentially without a proof of local optimality.
711  * Note, that in this case we report only feasibility and not optimality of the solution (see ScipNLP::finalize_solution).
712  */
713  (void) nlpiproblem->ipopt->Options()->SetNumericValue("dual_inf_tol", param.opttol);
714  (void) nlpiproblem->ipopt->Options()->SetNumericValue("compl_inf_tol", param.opttol);
715  if( param.solvertol > 0.0 )
716  (void) nlpiproblem->ipopt->Options()->SetNumericValue("tol", param.solvertol);
717  else
718 #if IPOPT_VERSION_MAJOR > 3 || IPOPT_VERSION_MINOR > 14 || (IPOPT_VERSION_MINOR == 14 && IPOPT_VERSION_RELEASE >= 2)
719  (void) nlpiproblem->ipopt->Options()->UnsetValue("tol");
720 #else
721  (void) nlpiproblem->ipopt->Options()->SetNumericValue("tol", 1e-8); // 1e-8 is the default
722 #endif
723 
724  /* Ipopt doesn't like a setting of exactly 0 for the max_*_time, so increase as little as possible in that case */
725 #if IPOPT_VERSION_MAJOR > 3 || IPOPT_VERSION_MINOR >= 14
726  (void) nlpiproblem->ipopt->Options()->SetNumericValue("max_wall_time", MAX(param.timelimit, DBL_MIN));
727 #else
728  (void) nlpiproblem->ipopt->Options()->SetNumericValue("max_cpu_time", MAX(param.timelimit, DBL_MIN));
729 #endif
730 
731  // disable acceptable-point heuristic iff fastfail is completely off
732  // by default (fastfail=conservative), it seems useful to have Ipopt stop when it obviously doesn't make progress (like one of the NLPs in the bendersqp ctest)
733  if( param.fastfail == SCIP_NLPPARAM_FASTFAIL_OFF )
734  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("acceptable_iter", 0);
735  else
736 #if IPOPT_VERSION_MAJOR > 3 || IPOPT_VERSION_MINOR > 14 || (IPOPT_VERSION_MINOR == 14 && IPOPT_VERSION_RELEASE >= 2)
737  (void) nlpiproblem->ipopt->Options()->UnsetValue("acceptable_iter");
738 #else
739  (void) nlpiproblem->ipopt->Options()->SetIntegerValue("acceptable_iter", 15); // 15 is the default
740 #endif
741 
742  (void) nlpiproblem->ipopt->Options()->SetStringValue("expect_infeasible_problem", param.expectinfeas ? "yes" : "no");
743 
744  if( !nlpiproblem->ipopt->Options()->SetStringValue("warm_start_init_point", param.warmstart ? "yes" : "no") && !param.warmstart )
745  {
746  // if we cannot disable warmstarts in Ipopt, then we have a big problem
747  SCIPerrorMessage("Failed to set Ipopt warm_start_init_point option to no.");
748  return SCIP_ERROR;
749  }
750 
751  return SCIP_OKAY;
752 }
753 
754 #ifdef COLLECT_SOLVESTATS
755 /// writes out some solve status, number of iterations, time, and problem properties
756 static
757 void collectStatistic(
758  SCIP* scip,
759  ApplicationReturnStatus status,
760  SCIP_NLPIPROBLEM* problem,
761  SmartPtr<SolveStatistics> stats
762  )
763 {
764  int nvars = 0; // number of variables, including slacks
765  int nslacks = 0; // number of slack variables
766  int nnlvars = 0; // number of nonlinear variables
767  int nvarlb = 0; // number of finite lower bounds, including slacks
768  int nvarub = 0; // number of finite upper bounds, including slacks
769  int nlincons = 0; // number of linear constraints
770  int nnlcons = 0; // number of nonlinear constraints
771  int objnl = 0; // number of nonlinear objectives
772  int jacnnz = 0; // number of nonzeros in Jacobian
773  int hesnnz = 0; // number of nonzeros in Hessian of Lagrangian
774  int linsys11nz = 0; // number of nonzeros in linear system (11) solved by Ipopt
775  int linsys13nz = 0; // number of nonzeros in linear system (13) solved by Ipopt
776 
777  int n = SCIPnlpiOracleGetNVars(problem->oracle);
778  int m = SCIPnlpiOracleGetNConstraints(problem->oracle);
779 
780  for( int i = 0; i < n; ++i )
781  {
782  SCIP_Real lb, ub;
783 
784  lb = SCIPnlpiOracleGetVarLbs(problem->oracle)[i];
785  ub = SCIPnlpiOracleGetVarUbs(problem->oracle)[i];
786 
787  // skip fixed vars
788  if( SCIPisEQ(scip, lb, ub) )
789  continue;
790 
791  ++nvars;
792 
793  if( SCIPnlpiOracleIsVarNonlinear(scip, problem->oracle, i) )
794  ++nnlvars;
795 
796  // every variable lower bound contributes a ln(x-lb) term in the barrier problem
797  if( !SCIPisInfinity(scip, -lb) )
798  ++nvarlb;
799 
800  // every variable upper bound contributes a ln(ub-x) term in the barrier problem
801  if( !SCIPisInfinity(scip, ub) )
802  ++nvarub;
803  }
804 
805  for( int i = 0; i < m; ++i )
806  {
807  SCIP_Real lhs, rhs;
808 
810  ++nnlcons;
811  else
812  ++nlincons;
813 
814  lhs = SCIPnlpiOracleGetConstraintLhs(problem->oracle, i);
815  rhs = SCIPnlpiOracleGetConstraintRhs(problem->oracle, i);
816 
817  // every inequality contributes a slack variable
818  if( !SCIPisEQ(scip, lhs, rhs) )
819  {
820  ++nvars;
821  ++nslacks;
822  }
823 
824  // every constraint lhs contributes a ln(slack-lhs) term in the barrier problem
825  if( !SCIPisInfinity(scip, -lhs) )
826  ++nvarlb;
827  // every constraint rhs contributes a ln(rhs-slack) term in the barrier problem
828  if( !SCIPisInfinity(scip, rhs) )
829  ++nvarub;
830  }
831 
832  objnl = SCIPnlpiOracleIsConstraintNonlinear(problem->oracle, -1);
833 
834  const int* offset;
835  const int* col;
836  SCIP_CALL_ABORT( SCIPnlpiOracleGetJacobianSparsity(scip, problem->oracle, &offset, NULL) );
837  jacnnz = offset[m];
838 
839  SCIP_CALL_ABORT( SCIPnlpiOracleGetHessianLagSparsity(scip, problem->oracle, &offset, &col) );
840  hesnnz = offset[n];
841 
842  // number of nonzeros of matrix in linear system of barrier problem ((11) in Ipopt paper):
843  // off-diagonal elements of Hessian of Lagrangian (W_k without diagonal)
844  // + number of variables (diagonal of W_k + Sigma_k)
845  // + 2*(elements in Jacobian + number of slacks) (A_k)
846  for( int i = 0; i < n; ++i )
847  for( int j = offset[i]; j < offset[i+1]; ++j )
848  if( col[j] != i )
849  linsys11nz += 2; // off-diagonal element of Lagrangian, once for each triangle of matrix
850  linsys11nz += nvars; // number of variables
851  linsys11nz += 2 * (jacnnz + nslacks); // because each slack var contributes one entry to the Jacobian
852 
853  // number of nonzeros of matrix in perturbed linear system of barrier problem ((13) in Ipopt paper):
854  // linsys11nz + number of constraints
855  linsys13nz = linsys11nz + m;
856 
857  SCIP_Real linsys11density = (SCIP_Real)linsys11nz / SQR((SCIP_Real)(nvars+m));
858  SCIP_Real linsys13density = (SCIP_Real)linsys13nz / SQR((SCIP_Real)(nvars+m));
859 
860  bool expectinfeas;
861  problem->ipopt->Options()->GetBoolValue("expect_infeasible_problem", expectinfeas, "");
862 
863  static bool firstwrite = true;
864  if( firstwrite )
865  {
866  printf("IPOPTSTAT status,iter,time,nvars,nnlvars,nvarlb,nvarub,nlincons,nnlcons,objnl,jacnnz,hesnnz,linsys11nz,linsys13nz,linsys11density,linsys13density,expectinfeas\n");
867  firstwrite = false;
868  }
869 
870  printf("IPOPTSTAT %d,%d,%g,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%f,%f,%d\n",
871  status, stats->IterationCount(), stats->TotalWallclockTime(),
872  nvars, nnlvars, nvarlb, nvarub, nlincons, nnlcons, objnl, jacnnz, hesnnz, linsys11nz, linsys13nz, linsys11density, linsys13density, expectinfeas);
873 }
874 #endif
875 
876 /** copy method of NLP interface (called when SCIP copies plugins) */
877 static
878 SCIP_DECL_NLPICOPY(nlpiCopyIpopt)
879 {
881 
882  return SCIP_OKAY;
883 }
884 
885 /** destructor of NLP interface to free nlpi data */
886 static
887 SCIP_DECL_NLPIFREE(nlpiFreeIpopt)
888 {
889  assert(nlpidata != NULL);
890 
891  delete *nlpidata;
892  *nlpidata = NULL;
893 
894  return SCIP_OKAY;
895 }
896 
897 /** gets pointer for NLP solver to do dirty stuff */
898 static
899 SCIP_DECL_NLPIGETSOLVERPOINTER(nlpiGetSolverPointerIpopt)
900 {
901  if( problem == NULL )
902  return NULL;
903 
904  return (void*)GetRawPtr(problem->ipopt);
905 }
906 
907 /** creates a problem instance */
908 static
909 SCIP_DECL_NLPICREATEPROBLEM(nlpiCreateProblemIpopt)
910 {
911  SCIP_NLPIDATA* data;
912 
913  assert(nlpi != NULL);
914  assert(problem != NULL);
915 
916  data = SCIPnlpiGetData(nlpi);
917  assert(data != NULL);
918 
919  SCIP_ALLOC( *problem = new SCIP_NLPIPROBLEM ); /*lint !e774*/
920 
921  SCIP_CALL( SCIPnlpiOracleCreate(scip, &(*problem)->oracle) );
922  SCIP_CALL( SCIPnlpiOracleSetProblemName(scip, (*problem)->oracle, name) );
923 
924  try
925  {
926  /* initialize IPOPT without default journal */
927  (*problem)->ipopt = new IpoptApplication(false);
928 
929  /* plugin our journal to get output through SCIP message handler */
930  SmartPtr<Journal> jrnl = new ScipJournal("console", J_ITERSUMMARY, scip);
931  jrnl->SetPrintLevel(J_DBG, J_NONE);
932  if( !(*problem)->ipopt->Jnlst()->AddJournal(jrnl) )
933  {
934  SCIPerrorMessage("Failed to register ScipJournal for IPOPT output.");
935  }
936 
937  /* initialize Ipopt/SCIP NLP interface */
938  (*problem)->nlp = new ScipNLP(*problem, scip);
939  }
940  catch( const std::bad_alloc& )
941  {
942  SCIPerrorMessage("Not enough memory to initialize Ipopt.\n");
943  return SCIP_NOMEMORY;
944  }
945 
946  for( size_t i = 0; i < sizeof(ipopt_string_params) / sizeof(const char*); ++i )
947  {
948  SCIP_PARAM* param;
950  char* paramval;
951 
952  strcpy(paramname, "nlpi/" NLPI_NAME "/");
953  strcat(paramname, ipopt_string_params[i]);
954  param = SCIPgetParam(scip, paramname);
955 
956  // skip parameters that we didn't add to SCIP because they didn't exist in this build of Ipopt
957  if( param == NULL )
958  continue;
959 
960  // if value wasn't left at the default, then pass to Ipopt and forbid overwriting
961  paramval = SCIPparamGetString(param);
962  assert(paramval != NULL);
963  if( *paramval != '\0' )
964  (void) (*problem)->ipopt->Options()->SetStringValue(ipopt_string_params[i], paramval, false);
965  }
966 
967  for( size_t i = 0; i < sizeof(ipopt_int_params) / sizeof(const char*); ++i )
968  {
969  SCIP_PARAM* param;
971  int paramval;
972 
973  strcpy(paramname, "nlpi/" NLPI_NAME "/");
974  strcat(paramname, ipopt_int_params[i]);
975  param = SCIPgetParam(scip, paramname);
976 
977  // skip parameters that we didn't add to SCIP because they didn't exist in this build of Ipopt
978  if( param == NULL )
979  continue;
980 
981  // if value wasn't left at the default, then pass to Ipopt and forbid overwriting
982  paramval = SCIPparamGetInt(param);
983  if( paramval != SCIPparamGetIntDefault(param) )
984  (void) (*problem)->ipopt->Options()->SetIntegerValue(ipopt_int_params[i], paramval, false);
985  }
986 
987 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
988  /* Turn off bound relaxation for older Ipopt, as solutions may be out of bounds by more than constr_viol_tol.
989  * For Ipopt 3.14, bounds are relaxed by at most constr_viol_tol, so can leave bound_relax_factor at its default.
990  */
991  (void) (*problem)->ipopt->Options()->SetNumericValue("bound_relax_factor", 0.0);
992 #endif
993 
994  /* modify Ipopt's default settings to what we believe is appropriate */
995  /* (*problem)->ipopt->Options()->SetStringValue("print_timing_statistics", "yes"); */
996 #ifdef SCIP_DEBUG
997  (void) (*problem)->ipopt->Options()->SetStringValue("print_user_options", "yes");
998 #endif
999  (void) (*problem)->ipopt->Options()->SetStringValue("sb", "yes");
1000  (void) (*problem)->ipopt->Options()->SetStringValueIfUnset("mu_strategy", "adaptive");
1001  (void) (*problem)->ipopt->Options()->SetIntegerValue("max_iter", INT_MAX);
1002  (void) (*problem)->ipopt->Options()->SetNumericValue("nlp_lower_bound_inf", -SCIPinfinity(scip), false);
1003  (void) (*problem)->ipopt->Options()->SetNumericValue("nlp_upper_bound_inf", SCIPinfinity(scip), false);
1004  (void) (*problem)->ipopt->Options()->SetNumericValue("diverging_iterates_tol", SCIPinfinity(scip), false);
1005 
1006  /* when warmstarting, then reduce how much Ipopt modified the starting point */
1007  (void) (*problem)->ipopt->Options()->SetNumericValue("warm_start_bound_push", data->warm_start_push);
1008  (void) (*problem)->ipopt->Options()->SetNumericValue("warm_start_bound_frac", data->warm_start_push);
1009  (void) (*problem)->ipopt->Options()->SetNumericValue("warm_start_slack_bound_push", data->warm_start_push);
1010  (void) (*problem)->ipopt->Options()->SetNumericValue("warm_start_slack_bound_frac", data->warm_start_push);
1011  (void) (*problem)->ipopt->Options()->SetNumericValue("warm_start_mult_bound_push", data->warm_start_push);
1012 
1013  /* apply user's given options file */
1014  assert(data->optfile != NULL);
1015  if( (*problem)->ipopt->Initialize(data->optfile) != Solve_Succeeded )
1016  {
1017  SCIPerrorMessage("Error during initialization of Ipopt using optionfile \"%s\"\n", data->optfile);
1018  return SCIP_ERROR;
1019  }
1020 
1021  return SCIP_OKAY;
1022 }
1023 
1024 /** free a problem instance */
1025 static
1026 SCIP_DECL_NLPIFREEPROBLEM(nlpiFreeProblemIpopt)
1027 {
1028  int n;
1029  int m;
1030 
1031  assert(nlpi != NULL);
1032  assert(problem != NULL);
1033  assert(*problem != NULL);
1034  assert((*problem)->oracle != NULL);
1035 
1036  n = SCIPnlpiOracleGetNVars((*problem)->oracle);
1037  m = SCIPnlpiOracleGetNConstraints((*problem)->oracle);
1038 
1039  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->solprimals, n);
1040  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->soldualcons, m);
1041  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->soldualvarlb, n);
1042  SCIPfreeBlockMemoryArrayNull(scip, &(*problem)->soldualvarub, n);
1043 
1044  SCIP_CALL( SCIPnlpiOracleFree(scip, &(*problem)->oracle) );
1045 
1046  if( (*problem)->randnumgen != NULL )
1047  {
1048  SCIPfreeRandom(scip, &(*problem)->randnumgen);
1049  }
1050 
1051  delete *problem;
1052  *problem = NULL;
1053 
1054  return SCIP_OKAY;
1055 }
1056 
1057 /** gets pointer to solver-internal problem instance to do dirty stuff */
1058 static
1059 SCIP_DECL_NLPIGETPROBLEMPOINTER(nlpiGetProblemPointerIpopt)
1060 {
1061  assert(nlpi != NULL);
1062  assert(problem != NULL);
1063 
1064  return GetRawPtr(problem->nlp);
1065 }
1066 
1067 /** add variables */
1068 static
1069 SCIP_DECL_NLPIADDVARS(nlpiAddVarsIpopt)
1070 {
1071  int oldnvars;
1072 
1073  assert(nlpi != NULL);
1074  assert(problem != NULL);
1075  assert(problem->oracle != NULL);
1076 
1077  oldnvars = SCIPnlpiOracleGetNVars(problem->oracle);
1078 
1079  SCIPfreeBlockMemoryArrayNull(scip, &problem->solprimals, oldnvars);
1080  SCIPfreeBlockMemoryArrayNull(scip, &problem->soldualvarlb, oldnvars);
1081  SCIPfreeBlockMemoryArrayNull(scip, &problem->soldualvarub, oldnvars);
1082  invalidateSolution(problem);
1083 
1084  SCIP_CALL( SCIPnlpiOracleAddVars(scip, problem->oracle, nvars, lbs, ubs, varnames) );
1085 
1086  problem->samestructure = false;
1087 
1088  return SCIP_OKAY;
1089 }
1090 
1091 /** add constraints */
1092 static
1093 SCIP_DECL_NLPIADDCONSTRAINTS(nlpiAddConstraintsIpopt)
1094 {
1095  int oldncons;
1096 
1097  assert(nlpi != NULL);
1098  assert(problem != NULL);
1099  assert(problem->oracle != NULL);
1100 
1101  oldncons = SCIPnlpiOracleGetNConstraints(problem->oracle);
1102 
1103  SCIPfreeBlockMemoryArrayNull(scip, &problem->soldualcons, oldncons);
1104  problem->soldualvalid = false;
1105  problem->soldualgiven = false;
1106 
1107  SCIP_CALL( SCIPnlpiOracleAddConstraints(scip, problem->oracle, nconss, lhss, rhss, nlininds, lininds, linvals, exprs, names) );
1108 
1109  problem->samestructure = false;
1110 
1111  return SCIP_OKAY;
1112 }
1113 
1114 /** sets or overwrites objective, a minimization problem is expected */
1115 static
1116 SCIP_DECL_NLPISETOBJECTIVE(nlpiSetObjectiveIpopt)
1117 {
1118  assert(nlpi != NULL);
1119  assert(problem != NULL);
1120  assert(problem->oracle != NULL);
1121 
1122  /* We pass the objective gradient in dense form to Ipopt, so if the sparsity of that gradient changes, we do not change the structure of the problem inside Ipopt.
1123  * However, if the sparsity of the Hessian matrix of the objective changes, then the sparsity pattern of the Hessian of the Lagrangian may change.
1124  * Thus, set samestructure=false if the objective was and/or becomes nonlinear, but leave samestructure untouched if it was and stays linear.
1125  */
1126  if( expr != NULL || SCIPnlpiOracleIsConstraintNonlinear(problem->oracle, -1) )
1127  problem->samestructure = false;
1128 
1129  SCIP_CALL( SCIPnlpiOracleSetObjective(scip, problem->oracle, constant, nlins, lininds, linvals, expr) );
1130 
1131  /* keep solution as valid, but reset solve status and objective value */
1132  invalidateSolved(problem);
1133 
1134  return SCIP_OKAY;
1135 }
1136 
1137 /** change variable bounds */
1138 static
1139 SCIP_DECL_NLPICHGVARBOUNDS(nlpiChgVarBoundsIpopt)
1140 {
1141  assert(nlpi != NULL);
1142  assert(problem != NULL);
1143  assert(problem->oracle != NULL);
1144 
1145  /* Check whether the structure of the Ipopt internal NLP changes, if problem->samestructure at the moment.
1146  * We need to check whether variables become fixed or unfixed and whether bounds are added or removed.
1147  *
1148  * Project primal solution onto new bounds if currently valid.
1149  */
1150  if( problem->samestructure || problem->solprimalvalid )
1151  {
1152  for( int i = 0; i < nvars; ++i )
1153  {
1154  SCIP_Real oldlb;
1155  SCIP_Real oldub;
1156  oldlb = SCIPnlpiOracleGetVarLbs(problem->oracle)[indices[i]];
1157  oldub = SCIPnlpiOracleGetVarUbs(problem->oracle)[indices[i]];
1158 
1159  if( (oldlb == oldub) != (lbs[i] == ubs[i]) ) /*lint !e777*/
1160  problem->samestructure = false;
1161  else if( SCIPisInfinity(scip, -oldlb) != SCIPisInfinity(scip, -lbs[i]) )
1162  problem->samestructure = false;
1163  else if( SCIPisInfinity(scip, oldub) != SCIPisInfinity(scip, ubs[i]) )
1164  problem->samestructure = false;
1165 
1166  if( problem->solprimalvalid )
1167  {
1168  assert(problem->solprimals != NULL);
1169  problem->solprimals[i] = MIN(MAX(problem->solprimals[indices[i]], lbs[i]), ubs[i]);
1170  }
1171  }
1172  }
1173 
1174  SCIP_CALL( SCIPnlpiOracleChgVarBounds(scip, problem->oracle, nvars, indices, lbs, ubs) );
1175 
1176  invalidateSolved(problem);
1177 
1178  return SCIP_OKAY;
1179 }
1180 
1181 /** change constraint bounds */
1182 static
1183 SCIP_DECL_NLPICHGCONSSIDES(nlpiChgConsSidesIpopt)
1184 {
1185  assert(nlpi != NULL);
1186  assert(problem != NULL);
1187  assert(problem->oracle != NULL);
1188 
1189  /* Check whether the structure of the Ipopt internal NLP changes, if problem->samestructure at the moment.
1190  * We need to check whether constraints change from equality to inequality and whether sides are added or removed.
1191  */
1192  for( int i = 0; i < nconss && problem->samestructure; ++i )
1193  {
1194  SCIP_Real oldlhs;
1195  SCIP_Real oldrhs;
1196  oldlhs = SCIPnlpiOracleGetConstraintLhs(problem->oracle, indices[i]);
1197  oldrhs = SCIPnlpiOracleGetConstraintRhs(problem->oracle, indices[i]);
1198 
1199  if( (oldlhs == oldrhs) != (lhss[i] == rhss[i]) ) /*lint !e777*/
1200  problem->samestructure = false;
1201  else if( SCIPisInfinity(scip, -oldlhs) != SCIPisInfinity(scip, -lhss[i]) )
1202  problem->samestructure = false;
1203  else if( SCIPisInfinity(scip, oldrhs) != SCIPisInfinity(scip, rhss[i]) )
1204  problem->samestructure = false;
1205  }
1206 
1207  SCIP_CALL( SCIPnlpiOracleChgConsSides(scip, problem->oracle, nconss, indices, lhss, rhss) );
1208 
1209  invalidateSolved(problem);
1210 
1211  return SCIP_OKAY;
1212 }
1213 
1214 /** delete a set of variables */
1215 static
1216 SCIP_DECL_NLPIDELVARSET(nlpiDelVarSetIpopt)
1217 {
1218  int nvars;
1219 
1220  assert(nlpi != NULL);
1221  assert(problem != NULL);
1222  assert(problem->oracle != NULL);
1223  assert(SCIPnlpiOracleGetNVars(problem->oracle) == dstatssize);
1224 
1225  SCIP_CALL( SCIPnlpiOracleDelVarSet(scip, problem->oracle, dstats) );
1226 
1227  nvars = SCIPnlpiOracleGetNVars(problem->oracle);
1228 
1229  if( problem->solprimalvalid || problem->soldualvalid )
1230  {
1231  // update existing solution, if valid
1232  assert(!problem->solprimalvalid || problem->solprimals != NULL);
1233  assert(!problem->soldualvalid || problem->soldualvarlb != NULL);
1234  assert(!problem->soldualvalid || problem->soldualvarub != NULL);
1235 
1236  int i;
1237  for( i = 0; i < dstatssize; ++i )
1238  {
1239  if( dstats[i] != -1 )
1240  {
1241  assert(dstats[i] >= 0);
1242  assert(dstats[i] < nvars);
1243  if( problem->solprimals != NULL )
1244  problem->solprimals[dstats[i]] = problem->solprimals[i];
1245  if( problem->soldualvarlb != NULL )
1246  {
1247  assert(problem->soldualvarub != NULL);
1248  problem->soldualvarlb[dstats[i]] = problem->soldualvarlb[i];
1249  problem->soldualvarub[dstats[i]] = problem->soldualvarub[i];
1250  }
1251  }
1252  }
1253  }
1254 
1255  /* resize solution point arrays */
1256  if( problem->solprimals != NULL )
1257  {
1258  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->solprimals, dstatssize, nvars) );
1259  }
1260  if( problem->soldualvarlb != NULL )
1261  {
1262  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->soldualvarlb, dstatssize, nvars) );
1263  }
1264  if( problem->soldualvarub != NULL )
1265  {
1266  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->soldualvarub, dstatssize, nvars) );
1267  }
1268 
1269  problem->samestructure = false;
1270 
1271  invalidateSolved(problem);
1272 
1273  return SCIP_OKAY;
1274 }
1275 
1276 /** delete a set of constraints */
1277 static
1278 SCIP_DECL_NLPIDELCONSSET(nlpiDelConstraintSetIpopt)
1279 {
1280  int ncons;
1281 
1282  assert(nlpi != NULL);
1283  assert(problem != NULL);
1284  assert(problem->oracle != NULL);
1285  assert(SCIPnlpiOracleGetNConstraints(problem->oracle) == dstatssize);
1286 
1287  SCIP_CALL( SCIPnlpiOracleDelConsSet(scip, problem->oracle, dstats) );
1288 
1289  ncons = SCIPnlpiOracleGetNConstraints(problem->oracle);
1290 
1291  if( problem->soldualvalid )
1292  {
1293  // update existing dual solution
1294  assert(problem->soldualcons != NULL);
1295 
1296  int i;
1297  for( i = 0; i < dstatssize; ++i )
1298  {
1299  if( dstats[i] != -1 )
1300  {
1301  assert(dstats[i] >= 0);
1302  assert(dstats[i] < ncons);
1303  problem->soldualcons[dstats[i]] = problem->soldualcons[i];
1304  }
1305  }
1306  }
1307 
1308  /* resize dual solution point array */
1309  if( problem->soldualcons != NULL )
1310  {
1311  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &problem->soldualcons, dstatssize, ncons) );
1312  }
1313 
1314  problem->samestructure = false;
1315 
1316  invalidateSolved(problem);
1317 
1318  return SCIP_OKAY;
1319 }
1320 
1321 /** change one linear coefficient in a constraint or objective */
1322 static
1323 SCIP_DECL_NLPICHGLINEARCOEFS(nlpiChgLinearCoefsIpopt)
1324 {
1325  assert(nlpi != NULL);
1326  assert(problem != NULL);
1327  assert(problem->oracle != NULL);
1328 
1329  SCIP_CALL( SCIPnlpiOracleChgLinearCoefs(scip, problem->oracle, idx, nvals, varidxs, vals) );
1330 
1331  problem->samestructure = false; // nonzero patterns may have changed; TODO SCIPnlpiOracleChgLinearCoefs() should let us know
1332  invalidateSolved(problem);
1333 
1334  return SCIP_OKAY;
1335 }
1336 
1337 /** replaces the expression tree of a constraint or objective */
1338 static
1339 SCIP_DECL_NLPICHGEXPR(nlpiChgExprIpopt)
1340 {
1341  assert(nlpi != NULL);
1342  assert(problem != NULL);
1343  assert(problem->oracle != NULL);
1344 
1345  SCIP_CALL( SCIPnlpiOracleChgExpr(scip, problem->oracle, idxcons, expr) );
1346 
1347  problem->samestructure = false; // nonzero patterns may have changed
1348  invalidateSolved(problem);
1349 
1350  return SCIP_OKAY;
1351 }
1352 
1353 /** change the constant offset in the objective */
1354 static
1355 SCIP_DECL_NLPICHGOBJCONSTANT(nlpiChgObjConstantIpopt)
1356 {
1357  SCIP_Real oldconstant;
1358 
1359  assert(nlpi != NULL);
1360  assert(problem != NULL);
1361  assert(problem->oracle != NULL);
1362 
1363  oldconstant = SCIPnlpiOracleGetObjectiveConstant(problem->oracle);
1364 
1365  SCIP_CALL( SCIPnlpiOracleChgObjConstant(scip, problem->oracle, objconstant) );
1366 
1367  if( problem->solobjval != SCIP_INVALID )
1368  problem->solobjval += objconstant - oldconstant;
1369 
1370  return SCIP_OKAY;
1371 }
1372 
1373 /** sets initial guess for primal variables */
1374 static
1375 SCIP_DECL_NLPISETINITIALGUESS(nlpiSetInitialGuessIpopt)
1376 {
1377  int nvars;
1378 
1379  assert(nlpi != NULL);
1380  assert(problem != NULL);
1381  assert(problem->oracle != NULL);
1382 
1383  nvars = SCIPnlpiOracleGetNVars(problem->oracle);
1384 
1385  if( primalvalues != NULL )
1386  {
1387  // copy primal solution
1388  SCIPdebugMsg(scip, "set initial guess primal values to user-given\n");
1389  if( problem->solprimals == NULL )
1390  {
1391  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->solprimals, nvars) );
1392  }
1393  BMScopyMemoryArray(problem->solprimals, primalvalues, nvars);
1394  problem->solprimalvalid = true;
1395  problem->solprimalgiven = true;
1396  }
1397  else
1398  {
1399  // invalid current primal solution (if any)
1400  if( problem->solprimalvalid )
1401  {
1402  SCIPdebugMsg(scip, "invalidate initial guess primal values on user-request\n");
1403  }
1404  problem->solprimalvalid = false;
1405  problem->solprimalgiven = false;
1406  }
1407 
1408  if( consdualvalues != NULL && varlbdualvalues != NULL && varubdualvalues != NULL )
1409  {
1410  // copy dual solution, if completely given
1411  SCIPdebugMsg(scip, "set initial guess dual values to user-given\n");
1412  int ncons = SCIPnlpiOracleGetNConstraints(problem->oracle);
1413  if( problem->soldualcons == NULL )
1414  {
1415  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->soldualcons, ncons) );
1416  }
1417  BMScopyMemoryArray(problem->soldualcons, consdualvalues, ncons);
1418 
1419  assert((problem->soldualvarlb == NULL) == (problem->soldualvarub == NULL));
1420  if( problem->soldualvarlb == NULL )
1421  {
1422  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->soldualvarlb, nvars) );
1423  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &problem->soldualvarub, nvars) );
1424  }
1425  BMScopyMemoryArray(problem->soldualvarlb, varlbdualvalues, nvars);
1426  BMScopyMemoryArray(problem->soldualvarub, varubdualvalues, nvars);
1427  problem->soldualvalid = true;
1428  problem->soldualgiven = true;
1429  }
1430  else
1431  {
1432  // invalid current dual solution (if any)
1433  if( problem->soldualvalid )
1434  {
1435  SCIPdebugMsg(scip, "invalidate initial guess dual values\n");
1436  }
1437  problem->soldualvalid = false;
1438  problem->soldualgiven = false;
1439  }
1440 
1441  return SCIP_OKAY;
1442 }
1443 
1444 /** tries to solve NLP */
1445 static
1446 SCIP_DECL_NLPISOLVE(nlpiSolveIpopt)
1447 {
1448  SCIP_NLPIDATA* nlpidata;
1449  ApplicationReturnStatus status;
1450 
1451  assert(nlpi != NULL);
1452  assert(problem != NULL);
1453  assert(problem->oracle != NULL);
1454 
1455  assert(IsValid(problem->ipopt));
1456  assert(IsValid(problem->nlp));
1457 
1458  nlpidata = SCIPnlpiGetData(nlpi);
1459  assert(nlpidata != NULL);
1460 
1461  // print parameters if either nlpi/ipopt/print_level has been set high enough or solve called with verblevel>0
1462  if( nlpidata->print_level >= J_SUMMARY || param.verblevel > 0 )
1463  {
1464  SCIPinfoMessage(scip, NULL, "Ipopt solve for problem %s at subSCIP depth %d", SCIPnlpiOracleGetProblemName(problem->oracle), SCIPgetSubscipDepth(scip));
1465  SCIPinfoMessage(scip, NULL, " with parameters " SCIP_NLPPARAM_PRINT(param));
1466  }
1467 
1468  SCIP_CALL( SCIPnlpiOracleResetEvalTime(scip, problem->oracle) );
1469 
1470  if( param.timelimit == 0.0 )
1471  {
1472  /* there is nothing we can do if we are not given any time */
1473  problem->lastniter = 0;
1474  problem->lasttime = 0.0;
1476  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
1477 
1478  return SCIP_OKAY;
1479  }
1480 
1481 #ifdef COLLECT_SOLVESTATS
1482  // if collecting statistics on how many iterations one may need for a certain Ipopt problem,
1483  // do not use a tight iteration limit (but use some limit to get finished)
1484  param.iterlimit = 1000;
1485 #endif
1486 
1487  // change status info to unsolved, just in case
1488  invalidateSolved(problem);
1489 
1490  // ensure a starting point is available
1491  // also disables param.warmstart if no warmstart available
1492  SCIP_CALL( ensureStartingPoint(scip, problem, param.warmstart) );
1493 
1494  // tell NLP that we are about to start a new solve
1495  problem->nlp->initializeSolve(problem, param);
1496 
1497  // set Ipopt parameters
1498  SCIP_CALL( handleNlpParam(scip, nlpidata, problem, param) );
1499 
1500  try
1501  {
1502 #ifdef PROTECT_SOLVE_BY_MUTEX
1503  /* lock solve_mutex if Ipopt is going to use Mumps as linear solver
1504  * unlocking will happen in the destructor of guard, which is called when this block is left
1505  */
1506  std::unique_lock<std::mutex> guard(solve_mutex, std::defer_lock); /*lint !e{728}*/
1507  std::string linsolver;
1508  (void) problem->ipopt->Options()->GetStringValue("linear_solver", linsolver, "");
1509  if( linsolver == "mumps" )
1510  guard.lock();
1511 #endif
1512 
1513  if( problem->firstrun )
1514  {
1516 
1518 
1519  /* if the expression interpreter or some user expression do not support function values and gradients and Hessians,
1520  * change NLP parameters or give an error
1521  */
1523  {
1526  {
1527  SCIPerrorMessage("Do not have expression interpreter that can compute function values and gradients. Cannot solve NLP with Ipopt.\n");
1528  problem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
1529  problem->termstat = SCIP_NLPTERMSTAT_OTHER;
1530  return SCIP_OKAY;
1531  }
1532 
1533  /* enable Hessian approximation if we are nonquadratic and the expression interpreter or user expression do not support Hessians */
1534  if( !(cap & SCIP_EXPRINTCAPABILITY_HESSIAN) )
1535  {
1536  (void) problem->ipopt->Options()->SetStringValueIfUnset("hessian_approximation", "limited-memory");
1537  problem->nlp->approxhessian = true;
1538  }
1539  else
1540  problem->nlp->approxhessian = false;
1541  }
1542 
1543 #ifdef SCIP_DEBUG
1544  problem->ipopt->Options()->SetStringValue("derivative_test", problem->nlp->approxhessian ? "first-order" : "second-order");
1545 #endif
1546 
1547  status = problem->ipopt->OptimizeTNLP(GetRawPtr(problem->nlp));
1548  }
1549  else
1550  {
1551  // TODO to be strict, we should check whether the eval capability has been changed and the Hessian approximation needs to be enabled (in which case we should call OptimizeTNLP instead)
1552  problem->ipopt->Options()->SetStringValue("warm_start_same_structure", problem->samestructure ? "yes" : "no");
1553  status = problem->ipopt->ReOptimizeTNLP(GetRawPtr(problem->nlp));
1554  }
1555 
1556  // catch the very bad status codes
1557  switch( status ) {
1558  // everything better than Not_Enough_Degrees_Of_Freedom is a non-serious error
1559  case Solve_Succeeded:
1560  case Solved_To_Acceptable_Level:
1561  case Infeasible_Problem_Detected:
1562  case Search_Direction_Becomes_Too_Small:
1563  case Diverging_Iterates:
1564  case User_Requested_Stop:
1565  case Feasible_Point_Found:
1566  case Maximum_Iterations_Exceeded:
1567  case Restoration_Failed:
1568  case Error_In_Step_Computation:
1569  case Maximum_CpuTime_Exceeded:
1570 #if IPOPT_VERSION_MAJOR > 3 || IPOPT_VERSION_MINOR >= 14
1571  case Maximum_WallTime_Exceeded: // new in Ipopt 3.14
1572  // if Ipopt >= 3.14, finalize_solution should always have been called if we get these status codes
1573  // this should have left us with some solution (unless we ran out of memory in finalize_solution)
1574  assert(problem->solprimalvalid || problem->termstat == SCIP_NLPTERMSTAT_OUTOFMEMORY);
1575  assert(problem->soldualvalid || problem->termstat == SCIP_NLPTERMSTAT_OUTOFMEMORY);
1576 #endif
1577  problem->firstrun = false;
1578  problem->samestructure = true;
1579  break;
1580 
1581  case Not_Enough_Degrees_Of_Freedom:
1582  assert(problem->termstat == SCIP_NLPTERMSTAT_OTHER);
1583  assert(problem->solstat == SCIP_NLPSOLSTAT_UNKNOWN);
1584  SCIPdebugMsg(scip, "NLP has too few degrees of freedom.\n");
1585  break;
1586 
1587  case Invalid_Number_Detected:
1588  SCIPdebugMsg(scip, "Ipopt failed because of an invalid number in function or derivative value\n");
1590  /* Ipopt may or may not have called finalize solution
1591  * if it didn't, then we should still have SCIP_NLPSOLSTAT_UNKNOWN as set in the invalidateSolved() call above
1592  * if it did, then finalize_solution will have set SCIP_NLPSOLSTAT_UNKNOWN or SCIP_NLPSOLSTAT_FEASIBLE
1593  */
1594  assert(problem->solstat == SCIP_NLPSOLSTAT_UNKNOWN || problem->solstat == SCIP_NLPSOLSTAT_FEASIBLE);
1595  break;
1596 
1597  case Insufficient_Memory:
1598  assert(problem->termstat == SCIP_NLPTERMSTAT_OTHER);
1599  assert(problem->solstat == SCIP_NLPSOLSTAT_UNKNOWN);
1600  SCIPerrorMessage("Ipopt returned with status \"Insufficient Memory\"\n");
1601  return SCIP_NOMEMORY;
1602 
1603  // really bad ones that could be something very unexpected going wrong within Ipopt
1604  case Unrecoverable_Exception:
1605  case Internal_Error:
1606  assert(problem->termstat == SCIP_NLPTERMSTAT_OTHER);
1607  assert(problem->solstat == SCIP_NLPSOLSTAT_UNKNOWN);
1608  SCIPerrorMessage("Ipopt returned with application return status %d\n", status);
1609  break;
1610 
1611  // the really bad ones that indicate rather a programming error
1612  case Invalid_Problem_Definition:
1613  case Invalid_Option:
1614  case NonIpopt_Exception_Thrown:
1615  assert(problem->termstat == SCIP_NLPTERMSTAT_OTHER);
1616  assert(problem->solstat == SCIP_NLPSOLSTAT_UNKNOWN);
1617  SCIPerrorMessage("Ipopt returned with application return status %d\n", status);
1618  return SCIP_ERROR;
1619  }
1620 
1621 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
1622  SmartPtr<SolveStatistics> stats = problem->ipopt->Statistics();
1623  /* Ipopt does not provide access to the statistics if there was a serious error */
1624  if( IsValid(stats) )
1625  {
1626  problem->lastniter = stats->IterationCount();
1627  problem->lasttime = stats->TotalWallclockTime();
1628 
1629 #ifdef COLLECT_SOLVESTATS
1630  collectStatistic(scip, status, problem, stats);
1631 #endif
1632  }
1633 #else
1634  SmartPtr<IpoptData> ip_data = problem->ipopt->IpoptDataObject();
1635  /* I don't think that there is a situation where ip_data is NULL, but check here anyway */
1636  if( IsValid(ip_data) )
1637  {
1638  problem->lastniter = ip_data->iter_count();
1639  problem->lasttime = ip_data->TimingStats().OverallAlgorithm().TotalWallclockTime();
1640  }
1641 #endif
1642  else
1643  {
1644  problem->lastniter = 0;
1645  problem->lasttime = 0.0;
1646  }
1647  }
1648  catch( IpoptException& except )
1649  {
1650  SCIPerrorMessage("Ipopt returned with exception: %s\n", except.Message().c_str());
1651  return SCIP_ERROR;
1652  }
1653 
1654  return SCIP_OKAY;
1655 }
1656 
1657 /** gives solution status */
1658 static
1659 SCIP_DECL_NLPIGETSOLSTAT(nlpiGetSolstatIpopt)
1660 {
1661  assert(nlpi != NULL);
1662  assert(problem != NULL);
1663 
1664  return problem->solstat;
1665 }
1666 
1667 /** gives termination reason */
1668 static
1669 SCIP_DECL_NLPIGETTERMSTAT(nlpiGetTermstatIpopt)
1670 {
1671  assert(nlpi != NULL);
1672  assert(problem != NULL);
1673 
1674  return problem->termstat;
1675 }
1676 
1677 /** gives primal and dual solution values */
1678 static
1679 SCIP_DECL_NLPIGETSOLUTION(nlpiGetSolutionIpopt)
1680 {
1681  assert(nlpi != NULL);
1682  assert(problem != NULL);
1683 
1684  if( primalvalues != NULL )
1685  *primalvalues = problem->solprimals;
1686 
1687  if( consdualvalues != NULL )
1688  *consdualvalues = problem->soldualcons;
1689 
1690  if( varlbdualvalues != NULL )
1691  *varlbdualvalues = problem->soldualvarlb;
1692 
1693  if( varubdualvalues != NULL )
1694  *varubdualvalues = problem->soldualvarub;
1695 
1696  if( objval != NULL )
1697  *objval = problem->solobjval;
1698 
1699  return SCIP_OKAY;
1700 }
1701 
1702 /** gives solve statistics */
1703 static
1704 SCIP_DECL_NLPIGETSTATISTICS(nlpiGetStatisticsIpopt)
1705 {
1706  assert(nlpi != NULL);
1707  assert(problem != NULL);
1708  assert(statistics != NULL);
1709 
1710  statistics->niterations = problem->lastniter;
1711  statistics->totaltime = problem->lasttime;
1712  statistics->evaltime = SCIPnlpiOracleGetEvalTime(scip, problem->oracle);
1713  statistics->consviol = problem->solconsviol;
1714  statistics->boundviol = problem->solboundviol;
1715 
1716  return SCIP_OKAY;
1717 }
1718 
1719 /** create solver interface for Ipopt solver and includes it into SCIP, if Ipopt is available */
1721  SCIP* scip /**< SCIP data structure */
1722  )
1723 {
1724  SCIP_NLPIDATA* nlpidata;
1725  SCIP_Bool advanced = FALSE;
1726 
1727  assert(scip != NULL);
1728 
1729  SCIP_ALLOC( nlpidata = new SCIP_NLPIDATA() ); /*lint !e774*/
1730 
1731  SCIP_CALL( SCIPincludeNlpi(scip,
1733  nlpiCopyIpopt, nlpiFreeIpopt, nlpiGetSolverPointerIpopt,
1734  nlpiCreateProblemIpopt, nlpiFreeProblemIpopt, nlpiGetProblemPointerIpopt,
1735  nlpiAddVarsIpopt, nlpiAddConstraintsIpopt, nlpiSetObjectiveIpopt,
1736  nlpiChgVarBoundsIpopt, nlpiChgConsSidesIpopt, nlpiDelVarSetIpopt, nlpiDelConstraintSetIpopt,
1737  nlpiChgLinearCoefsIpopt, nlpiChgExprIpopt,
1738  nlpiChgObjConstantIpopt, nlpiSetInitialGuessIpopt, nlpiSolveIpopt, nlpiGetSolstatIpopt, nlpiGetTermstatIpopt,
1739  nlpiGetSolutionIpopt, nlpiGetStatisticsIpopt,
1740  nlpidata) );
1741 
1743 
1744  SCIP_CALL( SCIPaddStringParam(scip, "nlpi/" NLPI_NAME "/optfile", "name of Ipopt options file",
1745  &nlpidata->optfile, FALSE, "", NULL, NULL) );
1746 
1747  SCIP_CALL( SCIPaddRealParam(scip, "nlpi/" NLPI_NAME "/warm_start_push", "amount (relative and absolute) by which starting point is moved away from bounds in warmstarts",
1748  &nlpidata->warm_start_push, FALSE, 1e-9, 0.0, 1.0, NULL, NULL) );
1749 
1750  SmartPtr<RegisteredOptions> reg_options = new RegisteredOptions();
1751  IpoptApplication::RegisterAllIpoptOptions(reg_options);
1752 
1753  for( size_t i = 0; i < sizeof(ipopt_string_params) / sizeof(const char*); ++i )
1754  {
1755  SmartPtr<const RegisteredOption> option = reg_options->GetOption(ipopt_string_params[i]);
1756 
1757  // skip options not available with this build of Ipopt
1758  if( !IsValid(option) )
1759  continue;
1760 
1761  assert(option->Type() == OT_String);
1762 
1763  // prefix parameter name with nlpi/ipopt
1764  std::string paramname("nlpi/" NLPI_NAME "/");
1765  paramname += option->Name();
1766 
1767  // initialize description with short description from Ipopt
1768  std::stringstream descr;
1769  descr << option->ShortDescription();
1770 
1771  // add valid values to description, if there are more than one
1772  // the only case where there are less than 2 valid strings should be when anything is valid (in which case there is one valid string with value "*")
1773  std::vector<RegisteredOption::string_entry> validvals = option->GetValidStrings();
1774  if( validvals.size() > 1 )
1775  {
1776  descr << " Valid values if not empty:";
1777  for( std::vector<RegisteredOption::string_entry>::iterator val = validvals.begin(); val != validvals.end(); ++val )
1778  descr << ' ' << val->value_;
1779  }
1780 
1781 #if IPOPT_VERSION_MAJOR > 3 || IPOPT_VERSION_MINOR >= 14
1782  // since Ipopt 3.14, Ipopt options have an advanced flag
1783  advanced = option->Advanced();
1784 #endif
1785 
1786  // we use the empty string as default to recognize later whether the user set has set the option
1787  SCIP_CALL( SCIPaddStringParam(scip, paramname.c_str(), descr.str().c_str(), NULL, advanced, "", NULL, NULL) );
1788  }
1789 
1790  for( size_t i = 0; i < sizeof(ipopt_int_params) / sizeof(const char*); ++i )
1791  {
1792  assert(i > 0 || strcmp(ipopt_int_params[0], "print_level") == 0); // we assume print_level at index 0
1793 
1794  SmartPtr<const RegisteredOption> option = reg_options->GetOption(ipopt_int_params[i]);
1795 
1796  // skip options not available with this build of Ipopt
1797  if( !IsValid(option) )
1798  continue;
1799 
1800  assert(option->Type() == OT_Integer);
1801 
1802  // prefix parameter name with nlpi/ipopt
1803  std::string paramname("nlpi/" NLPI_NAME "/");
1804  paramname += option->Name();
1805 
1806  int lower = option->LowerInteger();
1807  int upper = option->UpperInteger();
1808 
1809  // we use value lower-1 as signal that the option was not modified by the user
1810  // for that, we require a finite lower bound
1811  assert(lower > INT_MIN);
1812 
1813  // initialize description with short description from Ipopt
1814  std::stringstream descr;
1815  descr << option->ShortDescription();
1816  descr << ' ' << (lower-1) << " to use NLPI or Ipopt default.";
1817 
1818 #if IPOPT_VERSION_MAJOR > 3 || IPOPT_VERSION_MINOR >= 14
1819  // since Ipopt 3.14, Ipopt options have an advanced flag
1820  advanced = option->Advanced();
1821 #endif
1822 
1823  // we use the empty string as default to recognize later whether the user set has set the option
1824  SCIP_CALL( SCIPaddIntParam(scip, paramname.c_str(), descr.str().c_str(),
1825  i == 0 ? &nlpidata->print_level : NULL, advanced,
1826  lower-1, lower-1, upper, NULL, NULL) );
1827  }
1828 
1829  return SCIP_OKAY;
1830 } /*lint !e429 */
1831 
1832 /** gets string that identifies Ipopt (version number) */
1833 const char* SCIPgetSolverNameIpopt(void)
1834 {
1835  return "Ipopt " IPOPT_VERSION;
1836 }
1837 
1838 /** gets string that describes Ipopt */
1839 const char* SCIPgetSolverDescIpopt(void)
1840 {
1841  return "Interior Point Optimizer developed by A. Waechter et.al. (github.com/coin-or/Ipopt)";
1842 }
1843 
1844 /** returns whether Ipopt is available, i.e., whether it has been linked in */
1846 {
1847  return TRUE;
1848 }
1849 
1850 /** gives a pointer to the NLPIORACLE object stored in Ipopt-NLPI's NLPI problem data structure */
1852  SCIP_NLPIPROBLEM* nlpiproblem /**< NLP problem of Ipopt-NLPI */
1853  )
1854 {
1855  assert(nlpiproblem != NULL);
1856 
1857  return nlpiproblem->oracle;
1858 }
1859 
1860 /** Method to return some info about the nlp */
1861 bool ScipNLP::get_nlp_info(
1862  Index& n, /**< place to store number of variables */
1863  Index& m, /**< place to store number of constraints */
1864  Index& nnz_jac_g, /**< place to store number of nonzeros in jacobian */
1865  Index& nnz_h_lag, /**< place to store number of nonzeros in hessian */
1866  IndexStyleEnum& index_style /**< place to store used index style (0-based or 1-based) */
1867  )
1868 {
1869  const int* offset;
1870  SCIP_RETCODE retcode;
1871 
1872  assert(nlpiproblem != NULL);
1873  assert(nlpiproblem->oracle != NULL);
1874 
1875  n = SCIPnlpiOracleGetNVars(nlpiproblem->oracle);
1876  m = SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle);
1877 
1878  retcode = SCIPnlpiOracleGetJacobianSparsity(scip, nlpiproblem->oracle, &offset, NULL);
1879  if( retcode != SCIP_OKAY )
1880  return false;
1881  assert(offset != NULL);
1882  nnz_jac_g = offset[m];
1883 
1884  if( !approxhessian )
1885  {
1886  retcode = SCIPnlpiOracleGetHessianLagSparsity(scip, nlpiproblem->oracle, &offset, NULL);
1887  if( retcode != SCIP_OKAY )
1888  return false;
1889  assert(offset != NULL);
1890  nnz_h_lag = offset[n];
1891  }
1892  else
1893  {
1894  nnz_h_lag = 0;
1895  }
1896 
1897  index_style = TNLP::C_STYLE;
1898 
1899  return true;
1900 }
1901 
1902 /** Method to return the bounds for my problem */
1903 bool ScipNLP::get_bounds_info(
1904  Index n, /**< number of variables */
1905  Number* x_l, /**< buffer to store lower bounds on variables */
1906  Number* x_u, /**< buffer to store upper bounds on variables */
1907  Index m, /**< number of constraints */
1908  Number* g_l, /**< buffer to store lower bounds on constraints */
1909  Number* g_u /**< buffer to store lower bounds on constraints */
1910  )
1911 {
1912  const int* varlincounts;
1913  const int* varnlcounts;
1914 
1915  assert(nlpiproblem != NULL);
1916  assert(nlpiproblem->oracle != NULL);
1917 
1918  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
1919  assert(m == SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle));
1920 
1921  assert(SCIPnlpiOracleGetVarLbs(nlpiproblem->oracle) != NULL || n == 0);
1922  assert(SCIPnlpiOracleGetVarUbs(nlpiproblem->oracle) != NULL || n == 0);
1923 
1924  BMScopyMemoryArray(x_l, SCIPnlpiOracleGetVarLbs(nlpiproblem->oracle), n);
1925  BMScopyMemoryArray(x_u, SCIPnlpiOracleGetVarUbs(nlpiproblem->oracle), n);
1926 #ifndef NDEBUG
1927  for( int i = 0; i < n; ++i )
1928  assert(x_l[i] <= x_u[i]);
1929 #endif
1930 
1931  SCIPnlpiOracleGetVarCounts(scip, nlpiproblem->oracle, &varlincounts, &varnlcounts);
1932 
1933  /* Ipopt performs better when unused variables do not appear, which we can achieve by fixing them,
1934  * since Ipopts TNLPAdapter will hide them from Ipopts NLP. In the dual solution, bound multipliers (z_L, z_U)
1935  * for these variables should have value 0.0 (they are set to -grad Lagrangian).
1936  */
1937  for( int i = 0; i < n; ++i )
1938  {
1939  if( varlincounts[i] == 0 && varnlcounts[i] == 0 )
1940  {
1941  SCIPdebugMsg(scip, "fix unused variable x%d [%g,%g] to 0.0 or bound\n", i, x_l[i], x_u[i]);
1942  assert(x_l[i] <= x_u[i]);
1943  x_l[i] = x_u[i] = MAX(MIN(x_u[i], 0.0), x_l[i]);
1944  }
1945  }
1946 
1947  for( int i = 0; i < m; ++i )
1948  {
1949  g_l[i] = SCIPnlpiOracleGetConstraintLhs(nlpiproblem->oracle, i);
1950  g_u[i] = SCIPnlpiOracleGetConstraintRhs(nlpiproblem->oracle, i);
1951  assert(g_l[i] <= g_u[i]);
1952  }
1953 
1954  return true;
1955 }
1956 
1957 /** Method to return the starting point for the algorithm */ /*lint -e{715}*/
1958 bool ScipNLP::get_starting_point(
1959  Index n, /**< number of variables */
1960  bool init_x, /**< whether initial values for primal values are requested */
1961  Number* x, /**< buffer to store initial primal values */
1962  bool init_z, /**< whether initial values for dual values of variable bounds are requested */
1963  Number* z_L, /**< buffer to store dual values for variable lower bounds */
1964  Number* z_U, /**< buffer to store dual values for variable upper bounds */
1965  Index m, /**< number of constraints */
1966  bool init_lambda, /**< whether initial values for dual values of constraints are required */
1967  Number* lambda /**< buffer to store dual values of constraints */
1968  )
1969 { /*lint --e{715} */
1970  assert(nlpiproblem != NULL);
1971  assert(nlpiproblem->oracle != NULL);
1972 
1973  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
1974  assert(m == SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle));
1975 
1976  if( init_x )
1977  {
1978  assert(nlpiproblem->solprimalvalid);
1979  assert(nlpiproblem->solprimals != NULL);
1980  BMScopyMemoryArray(x, nlpiproblem->solprimals, n);
1981  }
1982 
1983  if( init_z )
1984  {
1985  assert(nlpiproblem->soldualvalid);
1986  assert(nlpiproblem->soldualvarlb != NULL);
1987  assert(nlpiproblem->soldualvarub != NULL);
1988  BMScopyMemoryArray(z_L, nlpiproblem->soldualvarlb, n);
1989  BMScopyMemoryArray(z_U, nlpiproblem->soldualvarub, n);
1990  }
1991 
1992  if( init_lambda )
1993  {
1994  assert(nlpiproblem->soldualvalid);
1995  assert(nlpiproblem->soldualcons != NULL);
1996  BMScopyMemoryArray(lambda, nlpiproblem->soldualcons, m);
1997  }
1998 
1999  return true;
2000 }
2001 
2002 /** Method to return the number of nonlinear variables. */
2003 Index ScipNLP::get_number_of_nonlinear_variables()
2004 {
2005  int count;
2006  int n;
2007 
2008  assert(nlpiproblem != NULL);
2009  assert(nlpiproblem->oracle != NULL);
2010 
2011  n = SCIPnlpiOracleGetNVars(nlpiproblem->oracle);
2012 
2013  count = 0;
2014  for( int i = 0; i < n; ++i )
2015  if( SCIPnlpiOracleIsVarNonlinear(scip, nlpiproblem->oracle, i) )
2016  ++count;
2017 
2018  return count;
2019 }
2020 
2021 /** Method to return the indices of the nonlinear variables */
2022 bool ScipNLP::get_list_of_nonlinear_variables(
2023  Index num_nonlin_vars, /**< number of nonlinear variables */
2024  Index* pos_nonlin_vars /**< array to fill with indices of nonlinear variables */
2025  )
2026 {
2027  int count;
2028  int n;
2029 
2030  assert(nlpiproblem != NULL);
2031  assert(nlpiproblem->oracle != NULL);
2032 
2033  n = SCIPnlpiOracleGetNVars(nlpiproblem->oracle);
2034 
2035  count = 0;
2036  for( int i = 0; i < n; ++i )
2037  {
2038  if( SCIPnlpiOracleIsVarNonlinear(scip, nlpiproblem->oracle, i) )
2039  {
2040  assert(count < num_nonlin_vars);
2041  pos_nonlin_vars[count++] = i;
2042  }
2043  }
2044 
2045  assert(count == num_nonlin_vars);
2046 
2047  return true;
2048 }
2049 
2050 /** Method to return metadata about variables and constraints */ /*lint -e{715}*/
2051 bool ScipNLP::get_var_con_metadata(
2052  Index n, /**< number of variables */
2053  StringMetaDataMapType& var_string_md, /**< variable meta data of string type */
2054  IntegerMetaDataMapType& var_integer_md,/**< variable meta data of integer type */
2055  NumericMetaDataMapType& var_numeric_md,/**< variable meta data of numeric type */
2056  Index m, /**< number of constraints */
2057  StringMetaDataMapType& con_string_md, /**< constraint meta data of string type */
2058  IntegerMetaDataMapType& con_integer_md,/**< constraint meta data of integer type */
2059  NumericMetaDataMapType& con_numeric_md /**< constraint meta data of numeric type */
2060  )
2061 { /*lint --e{715}*/
2062  assert(nlpiproblem != NULL);
2063  assert(nlpiproblem->oracle != NULL);
2064  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2065  assert(m == SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle));
2066 
2067  char** varnames = SCIPnlpiOracleGetVarNames(nlpiproblem->oracle);
2068  if( varnames != NULL )
2069  {
2070  std::vector<std::string>& varnamesvec(var_string_md["idx_names"]);
2071  varnamesvec.reserve((size_t)n);
2072  for( int i = 0; i < n; ++i )
2073  {
2074  if( varnames[i] != NULL )
2075  {
2076  varnamesvec.push_back(varnames[i]); /*lint !e3701*/
2077  }
2078  else
2079  {
2080  char buffer[20];
2081  (void) sprintf(buffer, "nlpivar%8d", i);
2082  varnamesvec.push_back(buffer);
2083  }
2084  }
2085  }
2086 
2087  std::vector<std::string>& consnamesvec(con_string_md["idx_names"]);
2088  consnamesvec.reserve((size_t)m);
2089  for( int i = 0; i < m; ++i )
2090  {
2091  if( SCIPnlpiOracleGetConstraintName(nlpiproblem->oracle, i) != NULL )
2092  {
2093  consnamesvec.push_back(SCIPnlpiOracleGetConstraintName(nlpiproblem->oracle, i));
2094  }
2095  else
2096  {
2097  char buffer[20];
2098  (void) sprintf(buffer, "nlpicons%8d", i);
2099  consnamesvec.push_back(buffer); /*lint !e3701*/
2100  }
2101  }
2102 
2103  return true;
2104 }
2105 
2106 /** Method to return the objective value */ /*lint -e{715}*/
2107 bool ScipNLP::eval_f(
2108  Index n, /**< number of variables */
2109  const Number* x, /**< point to evaluate */
2110  bool new_x, /**< whether some function evaluation method has been called for this point before */
2111  Number& obj_value /**< place to store objective function value */
2112  )
2113 { /*lint --e{715}*/
2114  assert(nlpiproblem != NULL);
2115  assert(nlpiproblem->oracle != NULL);
2116 
2117  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2118 
2119  if( new_x )
2120  ++current_x;
2121  last_f_eval_x = current_x;
2122 
2123  return SCIPnlpiOracleEvalObjectiveValue(scip, nlpiproblem->oracle, x, &obj_value) == SCIP_OKAY;
2124 }
2125 
2126 /** Method to return the gradient of the objective */ /*lint -e{715}*/
2127 bool ScipNLP::eval_grad_f(
2128  Index n, /**< number of variables */
2129  const Number* x, /**< point to evaluate */
2130  bool new_x, /**< whether some function evaluation method has been called for this point before */
2131  Number* grad_f /**< buffer to store objective gradient */
2132  )
2133 { /*lint --e{715}*/
2134  SCIP_Real dummy;
2135 
2136  assert(nlpiproblem != NULL);
2137  assert(nlpiproblem->oracle != NULL);
2138 
2139  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2140 
2141  if( new_x )
2142  ++current_x;
2143  else
2144  {
2145  // pass new_x = TRUE to objective gradient eval iff we have not evaluated the objective function at this point yet
2146  new_x = last_f_eval_x < current_x;
2147  }
2148  // if we evaluate the objective gradient with new_x = true, then this will also evaluate the objective function
2149  // (and if we do with new_x = false, then we already have last_f_eval_x == current_x anyway)
2150  last_f_eval_x = current_x;
2151 
2152  return SCIPnlpiOracleEvalObjectiveGradient(scip, nlpiproblem->oracle, x, new_x, &dummy, grad_f) == SCIP_OKAY;
2153 }
2154 
2155 /** Method to return the constraint residuals */ /*lint -e{715}*/
2156 bool ScipNLP::eval_g(
2157  Index n, /**< number of variables */
2158  const Number* x, /**< point to evaluate */
2159  bool new_x, /**< whether some function evaluation method has been called for this point before */
2160  Index m, /**< number of constraints */
2161  Number* g /**< buffer to store constraint function values */
2162  )
2163 { /*lint --e{715}*/
2164  assert(nlpiproblem != NULL);
2165  assert(nlpiproblem->oracle != NULL);
2166 
2167  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2168 
2169  if( new_x )
2170  ++current_x;
2171  last_g_eval_x = current_x;
2172 
2173  return SCIPnlpiOracleEvalConstraintValues(scip, nlpiproblem->oracle, x, g) == SCIP_OKAY;
2174 }
2175 
2176 /** Method to return:
2177  * 1) The structure of the jacobian (if "values" is NULL)
2178  * 2) The values of the jacobian (if "values" is not NULL)
2179  */ /*lint -e{715}*/
2180 bool ScipNLP::eval_jac_g(
2181  Index n, /**< number of variables */
2182  const Number* x, /**< point to evaluate */
2183  bool new_x, /**< whether some function evaluation method has been called for this point before */
2184  Index m, /**< number of constraints */
2185  Index nele_jac, /**< number of nonzero entries in jacobian */
2186  Index* iRow, /**< buffer to store row indices of nonzero jacobian entries, or NULL if values are requested */
2187  Index* jCol, /**< buffer to store column indices of nonzero jacobian entries, or NULL if values are requested */
2188  Number* values /**< buffer to store values of nonzero jacobian entries, or NULL if structure is requested */
2189  )
2190 { /*lint --e{715}*/
2191  assert(nlpiproblem != NULL);
2192  assert(nlpiproblem->oracle != NULL);
2193 
2194  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2195  assert(m == SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle));
2196 
2197  if( values == NULL )
2198  { /* Ipopt wants to know sparsity structure */
2199  const int* jacoffset;
2200  const int* jaccol;
2201  int j;
2202  int i;
2203 
2204  assert(iRow != NULL);
2205  assert(jCol != NULL);
2206 
2207  if( SCIPnlpiOracleGetJacobianSparsity(scip, nlpiproblem->oracle, &jacoffset, &jaccol) != SCIP_OKAY )
2208  return false;
2209 
2210  assert(jacoffset[0] == 0);
2211  assert(jacoffset[m] == nele_jac);
2212  j = jacoffset[0];
2213  for( i = 0; i < m; ++i )
2214  for( ; j < jacoffset[i+1]; ++j )
2215  iRow[j] = i;
2216 
2217  BMScopyMemoryArray(jCol, jaccol, nele_jac);
2218  }
2219  else
2220  {
2221  if( new_x )
2222  ++current_x;
2223  else
2224  {
2225  // pass new_x = TRUE to Jacobian eval iff we have not evaluated the constraint functions at this point yet
2226  new_x = last_g_eval_x < current_x;
2227  }
2228  // if we evaluate the Jacobian with new_x = true, then this will also evaluate the constraint functions
2229  // (and if we do with new_x = false, then we already have last_g_eval_x == current_x anyway)
2230  last_f_eval_x = current_x;
2231 
2232  if( SCIPnlpiOracleEvalJacobian(scip, nlpiproblem->oracle, x, new_x, NULL, values) != SCIP_OKAY )
2233  return false;
2234  }
2235 
2236  return true;
2237 }
2238 
2239 /** Method to return:
2240  * 1) The structure of the hessian of the lagrangian (if "values" is NULL)
2241  * 2) The values of the hessian of the lagrangian (if "values" is not NULL)
2242  */ /*lint -e{715}*/
2243 bool ScipNLP::eval_h(
2244  Index n, /**< number of variables */
2245  const Number* x, /**< point to evaluate */
2246  bool new_x, /**< whether some function evaluation method has been called for this point before */
2247  Number obj_factor, /**< weight for objective function */
2248  Index m, /**< number of constraints */
2249  const Number* lambda, /**< weights for constraint functions */
2250  bool new_lambda, /**< whether the hessian has been evaluated for these values of lambda before */
2251  Index nele_hess, /**< number of nonzero entries in hessian */
2252  Index* iRow, /**< buffer to store row indices of nonzero hessian entries, or NULL if values are requested */
2253  Index* jCol, /**< buffer to store column indices of nonzero hessian entries, or NULL if values are requested */
2254  Number* values /**< buffer to store values of nonzero hessian entries, or NULL if structure is requested */
2255  )
2256 { /*lint --e{715}*/
2257  assert(nlpiproblem != NULL);
2258  assert(nlpiproblem->oracle != NULL);
2259 
2260  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2261  assert(m == SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle));
2262 
2263  if( values == NULL )
2264  { /* Ipopt wants to know sparsity structure */
2265  const int* heslagoffset;
2266  const int* heslagcol;
2267  int j;
2268  int i;
2269 
2270  assert(iRow != NULL);
2271  assert(jCol != NULL);
2272 
2273  if( SCIPnlpiOracleGetHessianLagSparsity(scip, nlpiproblem->oracle, &heslagoffset, &heslagcol) != SCIP_OKAY )
2274  return false;
2275 
2276  assert(heslagoffset[0] == 0);
2277  assert(heslagoffset[n] == nele_hess);
2278  j = heslagoffset[0];
2279  for( i = 0; i < n; ++i )
2280  for( ; j < heslagoffset[i+1]; ++j )
2281  iRow[j] = i;
2282 
2283  BMScopyMemoryArray(jCol, heslagcol, nele_hess);
2284  }
2285  else
2286  {
2287  bool new_x_obj = new_x;
2288  bool new_x_cons = new_x;
2289  if( new_x )
2290  ++current_x;
2291  else
2292  {
2293  // pass new_x_obj = TRUE iff we have not evaluated the objective function at this point yet
2294  // pass new_x_cons = TRUE iff we have not evaluated the constraint functions at this point yet
2295  new_x_obj = last_f_eval_x < current_x;
2296  new_x_cons = last_g_eval_x < current_x;
2297  }
2298  // evaluating Hessians with new_x will also evaluate the functions itself
2299  last_f_eval_x = current_x;
2300  last_g_eval_x = current_x;
2301 
2302  if( SCIPnlpiOracleEvalHessianLag(scip, nlpiproblem->oracle, x, new_x_obj, new_x_cons, obj_factor, lambda, values) != SCIP_OKAY )
2303  return false;
2304  }
2305 
2306  return true;
2307 }
2308 
2309 /** Method called by the solver at each iteration.
2310  *
2311  * Checks whether SCIP solve is interrupted, objlimit is reached, or fastfail is triggered.
2312  * Sets solution and termination status accordingly.
2313  */ /*lint -e{715}*/
2314 bool ScipNLP::intermediate_callback(
2315  AlgorithmMode mode, /**< current mode of algorithm */
2316  Index iter, /**< current iteration number */
2317  Number obj_value, /**< current objective value */
2318  Number inf_pr, /**< current primal infeasibility */
2319  Number inf_du, /**< current dual infeasibility */
2320  Number mu, /**< current barrier parameter */
2321  Number d_norm, /**< current gradient norm */
2322  Number regularization_size,/**< current size of regularization */
2323  Number alpha_du, /**< current dual alpha */
2324  Number alpha_pr, /**< current primal alpha */
2325  Index ls_trials, /**< current number of linesearch trials */
2326  const IpoptData* ip_data, /**< pointer to Ipopt Data */
2327  IpoptCalculatedQuantities* ip_cq /**< pointer to current calculated quantities */
2328  )
2329 { /*lint --e{715}*/
2330  if( SCIPisSolveInterrupted(scip) )
2331  {
2332  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2333  nlpiproblem->termstat = SCIP_NLPTERMSTAT_INTERRUPT;
2334  return false;
2335  }
2336 
2337  /* feasible point with objective value below lower objective limit -> stop */
2338  if( obj_value <= param.lobjlimit && inf_pr <= param.feastol )
2339  {
2340  nlpiproblem->solstat = SCIP_NLPSOLSTAT_FEASIBLE;
2341  nlpiproblem->termstat = SCIP_NLPTERMSTAT_LOBJLIMIT;
2342  return false;
2343  }
2344 
2345  /* do convergence test if fastfail is enabled */
2346  if( param.fastfail >= SCIP_NLPPARAM_FASTFAIL_AGGRESSIVE )
2347  {
2348  int i;
2349 
2350  if( iter == 0 )
2351  {
2352  conv_lastrestoiter = -1;
2353  }
2354  else if( mode == RestorationPhaseMode )
2355  {
2356  conv_lastrestoiter = iter;
2357  }
2358  else if( conv_lastrestoiter == iter-1 )
2359  {
2360  /* just switched back from restoration mode, reset dual reduction targets */
2361  for( i = 0; i < convcheck_nchecks; ++i )
2362  conv_dutarget[i] = convcheck_minred[i] * inf_du;
2363  }
2364 
2365  if( iter == convcheck_startiter )
2366  {
2367  /* define initial targets and iteration limits */
2368  for( i = 0; i < convcheck_nchecks; ++i )
2369  {
2370  conv_prtarget[i] = convcheck_minred[i] * inf_pr;
2371  conv_dutarget[i] = convcheck_minred[i] * inf_du;
2372  conv_iterlim[i] = iter + convcheck_maxiter[i];
2373  }
2374  }
2375  else if( iter > convcheck_startiter )
2376  {
2377  /* check if we should stop */
2378  for( i = 0; i < convcheck_nchecks; ++i )
2379  {
2380  if( inf_pr <= conv_prtarget[i] )
2381  {
2382  /* sufficient reduction w.r.t. primal infeasibility target
2383  * reset target w.r.t. current infeasibilities
2384  */
2385  conv_prtarget[i] = convcheck_minred[i] * inf_pr;
2386  conv_dutarget[i] = convcheck_minred[i] * inf_du;
2387  conv_iterlim[i] = iter + convcheck_maxiter[i];
2388  }
2389  else if( iter >= conv_iterlim[i] )
2390  {
2391  /* we hit a limit, should we really stop? */
2392  SCIPdebugMsg(scip, "convcheck %d: inf_pr = %e > target %e; inf_du = %e target %e: ",
2393  i, inf_pr, conv_prtarget[i], inf_du, conv_dutarget[i]);
2394  if( mode == RegularMode && iter <= conv_lastrestoiter + convcheck_startiter )
2395  {
2396  /* if we returned from feasibility restoration recently, we allow some more iterations,
2397  * because Ipopt may go for optimality for some iterations, at the costs of infeasibility
2398  */
2399  SCIPdebugPrintf("continue, because restoration phase only %d iters ago\n", iter - conv_lastrestoiter);
2400  }
2401  else if( mode == RegularMode && inf_du <= conv_dutarget[i] && iter < conv_iterlim[i] + convcheck_maxiter[i] )
2402  {
2403  /* if dual reduction is sufficient, we allow for twice the number of iterations to reach primal infeas reduction */
2404  SCIPdebugPrintf("continue, because dual infeas. red. sufficient and only %d iters above limit\n", iter - conv_iterlim[i]);
2405  }
2406  else
2407  {
2408  SCIPdebugPrintf("abort solve\n");
2409  if( inf_pr <= param.feastol )
2410  nlpiproblem->solstat = SCIP_NLPSOLSTAT_FEASIBLE;
2411  else
2412  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2413  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OKAY;
2414  return false;
2415  }
2416  }
2417  }
2418  }
2419  }
2420 
2421  return true;
2422 }
2423 
2424 /** This method is called when the algorithm is complete so the TNLP can store/write the solution. */ /*lint -e{715}*/
2425 void ScipNLP::finalize_solution(
2426  SolverReturn status, /**< solve and solution status */
2427  Index n, /**< number of variables */
2428  const Number* x, /**< primal solution values */
2429  const Number* z_L, /**< dual values of variable lower bounds */
2430  const Number* z_U, /**< dual values of variable upper bounds */
2431  Index m, /**< number of constraints */
2432  const Number* g, /**< values of constraints */
2433  const Number* lambda, /**< dual values of constraints */
2434  Number obj_value, /**< objective function value */
2435  const IpoptData* data, /**< pointer to Ipopt Data */
2436  IpoptCalculatedQuantities* cq /**< pointer to calculated quantities */
2437  )
2438 { /*lint --e{715}*/
2439  assert(nlpiproblem != NULL);
2440  assert(nlpiproblem->oracle != NULL);
2441 
2442  assert(n == SCIPnlpiOracleGetNVars(nlpiproblem->oracle));
2443  assert(m == SCIPnlpiOracleGetNConstraints(nlpiproblem->oracle));
2444 
2445  bool check_feasibility = false; // whether we should check x for feasibility, if not NULL
2446  switch( status )
2447  {
2448  case SUCCESS:
2449  nlpiproblem->solstat = SCIP_NLPSOLSTAT_LOCOPT;
2450  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OKAY;
2451  assert(x != NULL);
2452  break;
2453 
2454  case STOP_AT_ACCEPTABLE_POINT:
2455  /* if stop at acceptable point, then dual infeasibility can be arbitrary large, so claim only feasibility */
2456  case FEASIBLE_POINT_FOUND:
2457  nlpiproblem->solstat = SCIP_NLPSOLSTAT_FEASIBLE;
2458  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OKAY;
2459  assert(x != NULL);
2460  break;
2461 
2462  case MAXITER_EXCEEDED:
2463  check_feasibility = true;
2464  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2465  nlpiproblem->termstat = SCIP_NLPTERMSTAT_ITERLIMIT;
2466  break;
2467 
2468  case CPUTIME_EXCEEDED:
2469  check_feasibility = true;
2470  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2471  nlpiproblem->termstat = SCIP_NLPTERMSTAT_TIMELIMIT;
2472  break;
2473 
2474  case STOP_AT_TINY_STEP:
2475  case RESTORATION_FAILURE:
2476  case ERROR_IN_STEP_COMPUTATION:
2477  check_feasibility = true;
2478  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2479  nlpiproblem->termstat = SCIP_NLPTERMSTAT_NUMERICERROR;
2480  break;
2481 
2482  case LOCAL_INFEASIBILITY:
2483  nlpiproblem->solstat = SCIP_NLPSOLSTAT_LOCINFEASIBLE;
2484  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OKAY;
2485  break;
2486 
2487  case USER_REQUESTED_STOP:
2488  // status codes already set in intermediate_callback
2489  break;
2490 
2491  case DIVERGING_ITERATES:
2492  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNBOUNDED;
2493  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OKAY;
2494  break;
2495 
2496  // for the following status codes, if we get called here at all,
2497  // then Ipopt passes zeros for duals and activities!
2498  // (see https://github.com/coin-or/Ipopt/blob/stable/3.14/src/Interfaces/IpIpoptApplication.cpp#L885-L934)
2499 
2500  case INVALID_NUMBER_DETECTED:
2501  // we can get this, if functions can still be evaluated, but are not differentiable
2502  // (so Ipopt couldn't check local optimality)
2503  // so we enable the check below for whether the point is feasible
2504  check_feasibility = true;
2505  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2506  nlpiproblem->termstat = SCIP_NLPTERMSTAT_EVALERROR;
2507  break;
2508 
2509  case TOO_FEW_DEGREES_OF_FREEDOM:
2510  case INTERNAL_ERROR:
2511  case INVALID_OPTION:
2512  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2513  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OTHER;
2514  break;
2515 
2516  case OUT_OF_MEMORY:
2517  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2518  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OUTOFMEMORY;
2519  break;
2520 
2521  default:
2522  SCIPerrorMessage("Ipopt returned with unknown solution status %d\n", status);
2523  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2524  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OTHER;
2525  break;
2526  }
2527 
2528  assert(x != NULL);
2529  assert(lambda != NULL);
2530  assert(z_L != NULL);
2531  assert(z_U != NULL);
2532 
2533  assert(nlpiproblem->solprimals != NULL);
2534 
2535  if( nlpiproblem->soldualcons == NULL )
2536  {
2537  (void) SCIPallocBlockMemoryArray(scip, &nlpiproblem->soldualcons, m);
2538  }
2539  if( nlpiproblem->soldualvarlb == NULL )
2540  {
2541  (void) SCIPallocBlockMemoryArray(scip, &nlpiproblem->soldualvarlb, n);
2542  }
2543  if( nlpiproblem->soldualvarub == NULL )
2544  {
2545  (void) SCIPallocBlockMemoryArray(scip, &nlpiproblem->soldualvarub, n);
2546  }
2547  if( nlpiproblem->soldualcons == NULL || nlpiproblem->soldualvarlb == NULL || nlpiproblem->soldualvarub == NULL )
2548  {
2549  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2550  nlpiproblem->termstat = SCIP_NLPTERMSTAT_OUTOFMEMORY;
2551  return;
2552  }
2553 
2554  BMScopyMemoryArray(nlpiproblem->solprimals, x, n);
2555  BMScopyMemoryArray(nlpiproblem->soldualcons, lambda, m);
2556  BMScopyMemoryArray(nlpiproblem->soldualvarlb, z_L, n);
2557  BMScopyMemoryArray(nlpiproblem->soldualvarub, z_U, n);
2558  nlpiproblem->solobjval = obj_value;
2559  nlpiproblem->solprimalvalid = true;
2560  nlpiproblem->solprimalgiven = false;
2561  nlpiproblem->soldualvalid = true;
2562  nlpiproblem->soldualgiven = false;
2563 
2564  // get violations, there could be an evaluation error when doing so
2565 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
2566  nlpiproblem->solboundviol = 0.0; // old Ipopt does not calculate bound violations, but for what it's worth, we have set bound_relax_factor=0 then
2567  if( cq == NULL )
2568  {
2569  // with old Ipopt, finalize_solution may be called with cq == NULL if all variables are fixed; we just skip the rest then
2570  nlpiproblem->solconsviol = 0.0;
2571  return;
2572  }
2573 #else
2574  assert(cq != NULL);
2575  nlpiproblem->solboundviol = cq->unscaled_curr_orig_bounds_violation(Ipopt::NORM_MAX);
2576 #endif
2577  try
2578  {
2579  nlpiproblem->solconsviol = cq->unscaled_curr_nlp_constraint_violation(Ipopt::NORM_MAX);
2580 
2581  if( check_feasibility )
2582  {
2583  // we assume that check_feasibility has not been enabled if Ipopt claimed infeasibility, since we should not change solstatus to unknown then
2584  assert(nlpiproblem->solstat != SCIP_NLPSOLSTAT_LOCINFEASIBLE);
2585  if( MAX(nlpiproblem->solconsviol, nlpiproblem->solboundviol) <= param.feastol )
2586  nlpiproblem->solstat = SCIP_NLPSOLSTAT_FEASIBLE;
2587  else
2588  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2589  }
2590  }
2591  catch( const IpoptNLP::Eval_Error& exc )
2592  {
2593  SCIPdebugMsg(scip, "Eval error when checking constraint viol: %s\n", exc.Message().c_str());
2594  assert(status == INVALID_NUMBER_DETECTED);
2595  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2596  nlpiproblem->solconsviol = SCIP_INVALID;
2597  }
2598 
2599  if( nlpiproblem->solstat == SCIP_NLPSOLSTAT_LOCINFEASIBLE )
2600  {
2601  assert(lambda != NULL);
2602  SCIP_Real tol;
2603  (void) nlpiproblem->ipopt->Options()->GetNumericValue("tol", tol, "");
2604 
2605  // Jakobs paper ZR_20-20 says we should have lambda*g(x) + mu*h(x) > 0
2606  // if the NLP is min f(x) s.t. g(x) <= 0, h(x) = 0
2607  // we check this here and change solution status to unknown if the test fails
2608  bool infreasonable = true;
2609  SCIP_Real infproof = 0.0;
2610  for( int i = 0; i < m && infreasonable; ++i )
2611  {
2612  if( fabs(lambda[i]) < tol )
2613  continue;
2614  SCIP_Real side;
2615  if( lambda[i] < 0.0 )
2616  {
2617  // lhs <= g(x) should be active
2618  // in the NLP above, this should be lhs - g(x) <= 0 with negated dual
2619  // so this contributes -lambda*(lhs-g(x)) = lambda*(g(x)-side)
2620  side = SCIPnlpiOracleGetConstraintLhs(nlpiproblem->oracle, i);
2621  if( SCIPisInfinity(scip, -side) )
2622  {
2623  SCIPdebugMessage("inconsistent dual, lambda = %g, but lhs = %g\n", lambda[i], side);
2624  infreasonable = false;
2625  }
2626  }
2627  else
2628  {
2629  // g(x) <= rhs should be active
2630  // in the NLP above, this should be g(x) - rhs <= 0
2631  // so this contributes lambda*(g(x)-rhs)
2632  side = SCIPnlpiOracleGetConstraintRhs(nlpiproblem->oracle, i);
2633  if( SCIPisInfinity(scip, side) )
2634  {
2635  SCIPdebugMessage("inconsistent dual, lambda = %g, but rhs = %g\n", lambda[i], side);
2636  infreasonable = false;
2637  }
2638  }
2639 
2640  // g(x) <= 0
2641  infproof += lambda[i] * (g[i] - side);
2642  // SCIPdebugMessage("cons %d lambda %g, slack %g\n", i, lambda[i], g[i] - side);
2643  }
2644  if( infreasonable )
2645  {
2646  SCIPdebugMessage("infproof = %g should be positive to be valid\n", infproof);
2647  if( infproof <= 0.0 )
2648  infreasonable = false;
2649  }
2650 
2651  if( !infreasonable )
2652  {
2653  // change status to say we don't know
2654  nlpiproblem->solstat = SCIP_NLPSOLSTAT_UNKNOWN;
2655  }
2656  }
2657 }
2658 
2659 /** Calls Lapacks Dsyev routine to compute eigenvalues and eigenvectors of a dense matrix.
2660  *
2661  * It's here, because we use Ipopt's C interface to Lapack.
2662  */
2664  SCIP_Bool computeeigenvectors,/**< should also eigenvectors should be computed ? */
2665  int N, /**< dimension */
2666  SCIP_Real* a, /**< matrix data on input (size N*N); eigenvectors on output if computeeigenvectors == TRUE */
2667  SCIP_Real* w /**< buffer to store eigenvalues (size N) */
2668  )
2669 {
2670  int info;
2671 
2672 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
2673  IpLapackDsyev((bool)computeeigenvectors, N, a, N, w, info);
2674 #else
2675  IpLapackSyev((bool)computeeigenvectors, N, a, N, w, info);
2676 #endif
2677 
2678  if( info != 0 )
2679  {
2680  SCIPerrorMessage("There was an error when calling DSYEV. INFO = %d\n", info);
2681  return SCIP_ERROR;
2682  }
2683 
2684  return SCIP_OKAY;
2685 }
2686 
2687 /** solves a linear problem of the form Ax = b for a regular matrix 3*3 A */
2688 static
2690  SCIP_Real* A, /**< matrix data on input (size 3*3); filled column-wise */
2691  SCIP_Real* b, /**< right hand side vector (size 3) */
2692  SCIP_Real* x, /**< buffer to store solution (size 3) */
2693  SCIP_Bool* success /**< pointer to store if the solving routine was successful */
2694  )
2695 {
2696  SCIP_Real Acopy[9];
2697  SCIP_Real bcopy[3];
2698  int pivotcopy[3];
2699  const int N = 3;
2700  int info;
2701 
2702  assert(A != NULL);
2703  assert(b != NULL);
2704  assert(x != NULL);
2705  assert(success != NULL);
2706 
2707  BMScopyMemoryArray(Acopy, A, N*N);
2708  BMScopyMemoryArray(bcopy, b, N);
2709 
2710  /* compute the LU factorization */
2711 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
2712  IpLapackDgetrf(N, Acopy, pivotcopy, N, info);
2713 #else
2714  IpLapackGetrf(N, Acopy, pivotcopy, N, info);
2715 #endif
2716 
2717  if( info != 0 )
2718  {
2719  SCIPdebugMessage("There was an error when calling Dgetrf. INFO = %d\n", info);
2720  *success = FALSE;
2721  }
2722  else
2723  {
2724  *success = TRUE;
2725 
2726  /* solve linear problem */
2727 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
2728  IpLapackDgetrs(N, 1, Acopy, N, pivotcopy, bcopy, N);
2729 #else
2730  IpLapackGetrs(N, 1, Acopy, N, pivotcopy, bcopy, N);
2731 #endif
2732 
2733  /* copy the solution */
2734  BMScopyMemoryArray(x, bcopy, N);
2735  }
2736 
2737  return SCIP_OKAY;
2738 }
2739 
2740 /** solves a linear problem of the form Ax = b for a regular matrix A
2741  *
2742  * Calls Lapacks DGETRF routine to calculate a LU factorization and uses this factorization to solve
2743  * the linear problem Ax = b.
2744  *
2745  * It's here, because we use Ipopt's C interface to Lapack.
2746  */
2748  int N, /**< dimension */
2749  SCIP_Real* A, /**< matrix data on input (size N*N); filled column-wise */
2750  SCIP_Real* b, /**< right hand side vector (size N) */
2751  SCIP_Real* x, /**< buffer to store solution (size N) */
2752  SCIP_Bool* success /**< pointer to store if the solving routine was successful */
2753  )
2754 {
2755  SCIP_Real* Acopy;
2756  SCIP_Real* bcopy;
2757  int* pivotcopy;
2758  int info;
2759 
2760  assert(N > 0);
2761  assert(A != NULL);
2762  assert(b != NULL);
2763  assert(x != NULL);
2764  assert(success != NULL);
2765 
2766  /* call solveLinearProb3() for performance reasons */
2767  if( N == 3 )
2768  {
2769  SCIP_CALL( solveLinearProb3(A, b, x, success) );
2770  return SCIP_OKAY;
2771  }
2772 
2773  Acopy = NULL;
2774  bcopy = NULL;
2775  pivotcopy = NULL;
2776 
2777  SCIP_ALLOC( BMSduplicateMemoryArray(&Acopy, A, N*N) );
2778  SCIP_ALLOC( BMSduplicateMemoryArray(&bcopy, b, N) );
2779  SCIP_ALLOC( BMSallocMemoryArray(&pivotcopy, N) );
2780 
2781  /* compute the LU factorization */
2782 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
2783  IpLapackDgetrf(N, Acopy, pivotcopy, N, info);
2784 #else
2785  IpLapackGetrf(N, Acopy, pivotcopy, N, info);
2786 #endif
2787 
2788  if( info != 0 )
2789  {
2790  SCIPdebugMessage("There was an error when calling Dgetrf. INFO = %d\n", info);
2791  *success = FALSE;
2792  }
2793  else
2794  {
2795  *success = TRUE;
2796 
2797  /* solve linear problem */
2798 #if IPOPT_VERSION_MAJOR == 3 && IPOPT_VERSION_MINOR < 14
2799  IpLapackDgetrs(N, 1, Acopy, N, pivotcopy, bcopy, N);
2800 #else
2801  IpLapackGetrs(N, 1, Acopy, N, pivotcopy, bcopy, N);
2802 #endif
2803 
2804  /* copy the solution */
2805  BMScopyMemoryArray(x, bcopy, N);
2806  }
2807 
2808  BMSfreeMemoryArray(&pivotcopy);
2809  BMSfreeMemoryArray(&bcopy);
2810  BMSfreeMemoryArray(&Acopy);
2811 
2812  return SCIP_OKAY;
2813 }
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
SCIP_RETCODE SCIPnlpiOracleEvalObjectiveGradient(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Bool isnewx, SCIP_Real *objval, SCIP_Real *objgrad)
Definition: nlpioracle.c:1959
#define SCIPreallocBlockMemoryArray(scip, ptr, oldnum, newnum)
Definition: scip_mem.h:90
enum SCIP_NlpTermStat SCIP_NLPTERMSTAT
Definition: type_nlpi.h:185
SCIP_NLPIDATA * SCIPnlpiGetData(SCIP_NLPI *nlpi)
Definition: nlpi.c:683
SCIP_RETCODE SCIPnlpiOracleChgLinearCoefs(SCIP *scip, SCIP_NLPIORACLE *oracle, int considx, int nentries, const int *varidxs, const SCIP_Real *newcoefs)
Definition: nlpioracle.c:1548
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:84
void * SCIPgetNlpiOracleIpopt(SCIP_NLPIPROBLEM *nlpiproblem)
void SCIPmessagePrintErrorHeader(const char *sourcefile, int sourceline)
Definition: message.c:768
public methods for SCIP parameter handling
static const char * ipopt_string_params[]
string parameters of Ipopt to make available via SCIP parameters
Definition: nlpi_ipopt.cpp:139
methods to interpret (evaluate) an expression "fast"
SmartPtr< ScipNLP > nlp
Definition: nlpi_ipopt.cpp:171
SCIP_Real opttol
Definition: type_nlpi.h:61
static SCIP_DECL_NLPICREATEPROBLEM(nlpiCreateProblemIpopt)
Definition: nlpi_ipopt.cpp:909
public methods for memory management
SCIP_Bool warmstart
Definition: type_nlpi.h:68
#define SCIP_MAXSTRLEN
Definition: def.h:293
static SCIP_DECL_NLPIFREEPROBLEM(nlpiFreeProblemIpopt)
#define NLPI_NAME
Definition: nlpi_ipopt.cpp:90
public solving methods
SCIP_RETCODE SCIPincludeNlpSolverIpopt(SCIP *scip)
SCIP_Real feastol
Definition: type_nlpi.h:60
static SCIP_RETCODE ensureStartingPoint(SCIP *scip, SCIP_NLPIPROBLEM *problem, SCIP_Bool &warmstart)
Definition: nlpi_ipopt.cpp:501
methods to store an NLP and request function, gradient, and Hessian values
SCIP_RETCODE SCIPnlpiOracleEvalHessianLag(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Bool isnewx_obj, SCIP_Bool isnewx_cons, SCIP_Real objfactor, const SCIP_Real *lambda, SCIP_Real *hessian)
Definition: nlpioracle.c:2371
SCIP_Real lasttime
Definition: nlpi_ipopt.cpp:190
SCIP_RETCODE SCIPcallLapackDsyevIpopt(SCIP_Bool computeeigenvectors, int N, SCIP_Real *a, SCIP_Real *w)
SCIP_NLPIORACLE * oracle
#define FALSE
Definition: def.h:87
unsigned short verblevel
Definition: type_nlpi.h:65
static SCIP_DECL_NLPIFREE(nlpiFreeIpopt)
Definition: nlpi_ipopt.cpp:887
static SCIP_DECL_NLPIGETSOLUTION(nlpiGetSolutionIpopt)
#define NLPI_DESC
Definition: nlpi_ipopt.cpp:91
int SCIPgetSubscipDepth(SCIP *scip)
Definition: scip_copy.c:2591
SCIP_RETCODE SCIPnlpiOracleAddConstraints(SCIP *scip, SCIP_NLPIORACLE *oracle, int nconss, const SCIP_Real *lhss, const SCIP_Real *rhss, const int *nlininds, int *const *lininds, SCIP_Real *const *linvals, SCIP_EXPR **exprs, const char **consnames)
Definition: nlpioracle.c:1158
SCIP_Real solvertol
Definition: type_nlpi.h:62
static SCIP_RETCODE solveLinearProb3(SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
SCIP_Real SCIPinfinity(SCIP *scip)
#define FEASTOLFACTOR
Definition: nlpi_ipopt.cpp:95
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_Real * soldualvarub
Definition: nlpi_ipopt.cpp:185
SCIP_RETCODE SCIPnlpiOracleSetProblemName(SCIP *scip, SCIP_NLPIORACLE *oracle, const char *name)
Definition: nlpioracle.c:1036
static SCIP_RETCODE handleNlpParam(SCIP *scip, SCIP_NLPIDATA *nlpidata, SCIP_NLPIPROBLEM *nlpiproblem, const SCIP_NLPPARAM param)
pass NLP solve parameters to Ipopt
Definition: nlpi_ipopt.cpp:571
SCIP_Real SCIPnlpiOracleGetEvalTime(SCIP *scip, SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:2434
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:116
static SCIP_DECL_NLPISETINITIALGUESS(nlpiSetInitialGuessIpopt)
void SCIPmessageVPrintInfo(SCIP_MESSAGEHDLR *messagehdlr, const char *formatstr, va_list ap)
Definition: message.c:599
SCIP_RETCODE SCIPaddStringParam(SCIP *scip, const char *name, const char *desc, char **valueptr, SCIP_Bool isadvanced, const char *defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:185
#define SCIPdebugMessage
Definition: pub_message.h:87
SCIP_RETCODE SCIPnlpiOracleSetObjective(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real constant, int nlin, const int *lininds, const SCIP_Real *linvals, SCIP_EXPR *expr)
Definition: nlpioracle.c:1219
SCIP_Real solconsviol
Definition: nlpi_ipopt.cpp:187
SCIP_MESSAGEHDLR * SCIPgetMessagehdlr(SCIP *scip)
Definition: scip_message.c:79
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
unsigned int SCIP_EXPRINTCAPABILITY
static SCIP_DECL_NLPIGETPROBLEMPOINTER(nlpiGetProblemPointerIpopt)
static SCIP_DECL_NLPISETOBJECTIVE(nlpiSetObjectiveIpopt)
SCIP_Real SCIPnlpiOracleGetConstraintRhs(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1812
#define SCIPdebugMsg
Definition: scip_message.h:69
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:74
SCIP_VAR ** x
Definition: circlepacking.c:54
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_Real timelimit
Definition: type_nlpi.h:63
public methods for numerical tolerances
SCIP_Real solobjval
Definition: nlpi_ipopt.cpp:186
SCIP_RANDNUMGEN * randnumgen
Definition: nlpi_ipopt.cpp:168
#define DEFAULT_RANDSEED
Definition: nlpi_ipopt.cpp:97
void SCIPmessageVPrintError(const char *formatstr, va_list ap)
Definition: message.c:795
public methods for NLPI solver interfaces
SCIP_VAR * w
Definition: circlepacking.c:58
char * SCIPparamGetString(SCIP_PARAM *param)
Definition: paramset.c:902
public methods for handling parameter settings
static SCIP_DECL_NLPIGETSTATISTICS(nlpiGetStatisticsIpopt)
#define SUCCESS
Definition: portab.h:45
const char * SCIPgetSolverNameIpopt(void)
#define BMSfreeMemoryArray(ptr)
Definition: memory.h:140
static SCIP_DECL_NLPICHGOBJCONSTANT(nlpiChgObjConstantIpopt)
#define SCIPerrorMessage
Definition: pub_message.h:55
#define NLPI_PRIORITY
Definition: nlpi_ipopt.cpp:92
int SCIPnlpiOracleGetNConstraints(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1717
#define SCIPdebugPrintf
Definition: pub_message.h:90
struct SCIP_NlpiData SCIP_NLPIDATA
Definition: type_nlpi.h:43
enum SCIP_NlpSolStat SCIP_NLPSOLSTAT
Definition: type_nlpi.h:159
#define SCIP_NLPPARAM_PRINT(param)
Definition: type_nlpi.h:133
int SCIPnlpiOracleGetNVars(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1707
SCIP_RETCODE SCIPnlpiOracleResetEvalTime(SCIP *scip, SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:2418
char * SCIPnlpiOracleGetConstraintName(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1825
SCIP_RETCODE SCIPnlpiOracleChgVarBounds(SCIP *scip, SCIP_NLPIORACLE *oracle, int nvars, const int *indices, const SCIP_Real *lbs, const SCIP_Real *ubs)
Definition: nlpioracle.c:1248
SCIP_Bool expectinfeas
Definition: type_nlpi.h:67
SCIP_Bool warmstart
SCIP_Real SCIPnlpiOracleGetObjectiveConstant(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1788
#define NULL
Definition: lpi_spx1.cpp:155
SCIP_PARAM * SCIPgetParam(SCIP *scip, const char *name)
Definition: scip_param.c:225
static void invalidateSolution(SCIP_NLPIPROBLEM *problem)
Definition: nlpi_ipopt.cpp:485
static const char * ipopt_int_params[]
integer parameters of Ipopt to make available via SCIP parameters
Definition: nlpi_ipopt.cpp:135
SCIP_RETCODE SCIPnlpiOracleGetJacobianSparsity(SCIP *scip, SCIP_NLPIORACLE *oracle, const int **offset, const int **col)
Definition: nlpioracle.c:2018
public methods for problem copies
SCIP_RETCODE SCIPnlpiOracleAddVars(SCIP *scip, SCIP_NLPIORACLE *oracle, int nvars, const SCIP_Real *lbs, const SCIP_Real *ubs, const char **varnames)
Definition: nlpioracle.c:1072
SCIP_RETCODE SCIPsolveLinearEquationsIpopt(int N, SCIP_Real *A, SCIP_Real *b, SCIP_Real *x, SCIP_Bool *success)
#define SCIP_CALL(x)
Definition: def.h:384
static SCIP_DECL_NLPICHGVARBOUNDS(nlpiChgVarBoundsIpopt)
SCIP_RETCODE SCIPnlpiOracleChgExpr(SCIP *scip, SCIP_NLPIORACLE *oracle, int considx, SCIP_EXPR *expr)
Definition: nlpioracle.c:1644
SCIP_EXPRINTCAPABILITY SCIPexprintGetCapability(void)
static SCIP_DECL_NLPICOPY(nlpiCopyIpopt)
Definition: nlpi_ipopt.cpp:878
SCIP_Real SCIPnlpiOracleGetConstraintLhs(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1799
#define BMSduplicateMemoryArray(ptr, source, num)
Definition: memory.h:136
SCIP_Real * soldualvarlb
Definition: nlpi_ipopt.cpp:184
static const int convcheck_nchecks
Definition: nlpi_ipopt.cpp:129
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
Ipopt NLP interface.
SmartPtr< IpoptApplication > ipopt
Definition: nlpi_ipopt.cpp:170
public data structures and miscellaneous methods
SCIP_RETCODE SCIPnlpiOracleEvalJacobian(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Bool isnewx, SCIP_Real *convals, SCIP_Real *jacobi)
Definition: nlpioracle.c:2150
#define SCIP_EXPRINTCAPABILITY_HESSIAN
SCIP_RETCODE SCIPnlpiOracleDelVarSet(SCIP *scip, SCIP_NLPIORACLE *oracle, int *delstats)
Definition: nlpioracle.c:1320
#define SCIP_Bool
Definition: def.h:84
#define MAXPERTURB
Definition: nlpi_ipopt.cpp:94
const char * SCIPnlpiOracleGetProblemName(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1060
static const char * paramname[]
Definition: lpi_msk.c:4998
static SCIP_DECL_NLPICHGLINEARCOEFS(nlpiChgLinearCoefsIpopt)
static const int convcheck_startiter
Definition: nlpi_ipopt.cpp:130
static SCIP_DECL_NLPIGETSOLSTAT(nlpiGetSolstatIpopt)
#define MAX(x, y)
Definition: tclique_def.h:83
SCIP_RETCODE SCIPnlpiOracleCreate(SCIP *scip, SCIP_NLPIORACLE **oracle)
Definition: nlpioracle.c:974
static SCIP_DECL_NLPIGETTERMSTAT(nlpiGetTermstatIpopt)
const SCIP_Real * SCIPnlpiOracleGetVarLbs(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1727
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:127
static SCIP_DECL_NLPIADDCONSTRAINTS(nlpiAddConstraintsIpopt)
SCIP_Bool SCIPnlpiOracleIsConstraintNonlinear(SCIP_NLPIORACLE *oracle, int considx)
Definition: nlpioracle.c:1838
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
const SCIP_Real * SCIPnlpiOracleGetVarUbs(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1737
SCIP_RETCODE SCIPnlpiOracleFree(SCIP *scip, SCIP_NLPIORACLE **oracle)
Definition: nlpioracle.c:1004
SCIP_Real SCIPrandomGetReal(SCIP_RANDNUMGEN *randnumgen, SCIP_Real minrandval, SCIP_Real maxrandval)
Definition: misc.c:10025
int SCIPparamGetInt(SCIP_PARAM *param)
Definition: paramset.c:725
static SCIP_DECL_NLPIGETSOLVERPOINTER(nlpiGetSolverPointerIpopt)
Definition: nlpi_ipopt.cpp:899
SCIP_VAR ** b
Definition: circlepacking.c:56
SCIP_Bool SCIPnlpiOracleIsVarNonlinear(SCIP *scip, SCIP_NLPIORACLE *oracle, int varidx)
Definition: nlpioracle.c:1757
general public methods
static void invalidateSolved(SCIP_NLPIPROBLEM *problem)
Definition: nlpi_ipopt.cpp:470
#define SCIP_EXPRINTCAPABILITY_FUNCVALUE
SCIP_RETCODE SCIPnlpiOracleEvalObjectiveValue(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Real *objval)
Definition: nlpioracle.c:1878
public methods for random numbers
SCIP_VAR * a
Definition: circlepacking.c:57
static SCIP_DECL_NLPICHGEXPR(nlpiChgExprIpopt)
const char * SCIPgetSolverDescIpopt(void)
static SCIP_DECL_NLPIDELVARSET(nlpiDelVarSetIpopt)
SCIP_RETCODE SCIPnlpiOracleEvalConstraintValues(SCIP *scip, SCIP_NLPIORACLE *oracle, const SCIP_Real *x, SCIP_Real *convals)
Definition: nlpioracle.c:1927
#define SCIP_Real
Definition: def.h:177
SCIP_RETCODE SCIPnlpiOracleGetHessianLagSparsity(SCIP *scip, SCIP_NLPIORACLE *oracle, const int **offset, const int **col)
Definition: nlpioracle.c:2277
int SCIPparamGetIntDefault(SCIP_PARAM *param)
Definition: paramset.c:761
SCIP_EXPRINTCAPABILITY SCIPnlpiOracleGetEvalCapability(SCIP *scip, SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1855
public methods for message handling
#define SCIP_INVALID
Definition: def.h:197
static SCIP_DECL_NLPISOLVE(nlpiSolveIpopt)
static SCIP_DECL_NLPIDELCONSSET(nlpiDelConstraintSetIpopt)
SCIP_RETCODE SCIPincludeExternalCodeInformation(SCIP *scip, const char *name, const char *description)
Definition: scip_general.c:704
SCIP_NLPTERMSTAT termstat
SCIP_NLPSOLSTAT solstat
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:102
SCIP_RETCODE SCIPincludeNlpi(SCIP *scip, const char *name, const char *description, int priority, SCIP_DECL_NLPICOPY((*nlpicopy)), SCIP_DECL_NLPIFREE((*nlpifree)), SCIP_DECL_NLPIGETSOLVERPOINTER((*nlpigetsolverpointer)), SCIP_DECL_NLPICREATEPROBLEM((*nlpicreateproblem)), SCIP_DECL_NLPIFREEPROBLEM((*nlpifreeproblem)), SCIP_DECL_NLPIGETPROBLEMPOINTER((*nlpigetproblempointer)), SCIP_DECL_NLPIADDVARS((*nlpiaddvars)), SCIP_DECL_NLPIADDCONSTRAINTS((*nlpiaddconstraints)), SCIP_DECL_NLPISETOBJECTIVE((*nlpisetobjective)), SCIP_DECL_NLPICHGVARBOUNDS((*nlpichgvarbounds)), SCIP_DECL_NLPICHGCONSSIDES((*nlpichgconssides)), SCIP_DECL_NLPIDELVARSET((*nlpidelvarset)), SCIP_DECL_NLPIDELCONSSET((*nlpidelconsset)), SCIP_DECL_NLPICHGLINEARCOEFS((*nlpichglinearcoefs)), SCIP_DECL_NLPICHGEXPR((*nlpichgexpr)), SCIP_DECL_NLPICHGOBJCONSTANT((*nlpichgobjconstant)), SCIP_DECL_NLPISETINITIALGUESS((*nlpisetinitialguess)), SCIP_DECL_NLPISOLVE((*nlpisolve)), SCIP_DECL_NLPIGETSOLSTAT((*nlpigetsolstat)), SCIP_DECL_NLPIGETTERMSTAT((*nlpigettermstat)), SCIP_DECL_NLPIGETSOLUTION((*nlpigetsolution)), SCIP_DECL_NLPIGETSTATISTICS((*nlpigetstatistics)), SCIP_NLPIDATA *nlpidata)
Definition: scip_nlpi.c:98
static const SCIP_Real convcheck_minred[convcheck_nchecks]
Definition: nlpi_ipopt.cpp:132
#define SCIP_CALL_ABORT(x)
Definition: def.h:363
char ** SCIPnlpiOracleGetVarNames(SCIP_NLPIORACLE *oracle)
Definition: nlpioracle.c:1747
SCIP_Real * solprimals
Definition: nlpi_ipopt.cpp:182
SCIP_Real * soldualcons
Definition: nlpi_ipopt.cpp:183
#define SCIP_ALLOC(x)
Definition: def.h:395
static SCIP_DECL_NLPICHGCONSSIDES(nlpiChgConsSidesIpopt)
SCIP_RETCODE SCIPnlpiOracleChgConsSides(SCIP *scip, SCIP_NLPIORACLE *oracle, int nconss, const int *indices, const SCIP_Real *lhss, const SCIP_Real *rhss)
Definition: nlpioracle.c:1285
SCIP_RETCODE SCIPnlpiOracleDelConsSet(SCIP *scip, SCIP_NLPIORACLE *oracle, int *delstats)
Definition: nlpioracle.c:1462
static SCIP_DECL_NLPIADDVARS(nlpiAddVarsIpopt)
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:130
SCIP_NLPPARAM_FASTFAIL fastfail
Definition: type_nlpi.h:66
SCIP_Real SCIPfloor(SCIP *scip, SCIP_Real val)
#define SCIP_EXPRINTCAPABILITY_GRADIENT
SCIP_Real solboundviol
Definition: nlpi_ipopt.cpp:188
void SCIPnlpiOracleGetVarCounts(SCIP *scip, SCIP_NLPIORACLE *oracle, const int **lincounts, const int **nlcounts)
Definition: nlpioracle.c:1772
static const int convcheck_maxiter[convcheck_nchecks]
Definition: nlpi_ipopt.cpp:131
SCIP_Bool SCIPisSolveInterrupted(SCIP *scip)
Definition: scip_solve.c:3548
SCIP_RETCODE SCIPnlpiOracleChgObjConstant(SCIP *scip, SCIP_NLPIORACLE *oracle, SCIP_Real objconstant)
Definition: nlpioracle.c:1690