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