1 // (C) Copyright International Business Machines Corporation and
2 // Carnegie Mellon University 2004, 2007
3 //
4 // All Rights Reserved.
5 // This code is published under the Eclipse Public License.
6 //
7 // Authors :
8 // Pierre Bonami, Carnegie Mellon University,
9 // Carl D. Laird, Carnegie Mellon University,
10 // Andreas Waechter, International Business Machines Corporation
11 //
12 // Date : 12/01/2004
13
14 #include "BonminConfig.h"
15
16 #include "BonOsiTMINLPInterface.hpp"
17 #include "BonTMINLP.hpp"
18 #include "BonTMINLP2TNLP.hpp"
19 #include "BonTNLP2FPNLP.hpp"
20 #include "BonTMINLP2OsiLP.hpp"
21 #include "BonTNLPSolver.hpp"
22 #include "CoinTime.hpp"
23 #include <climits>
24 #include <string>
25 #include <sstream>
26 #include "BonAuxInfos.hpp"
27
28 #include "Ipopt/BonIpoptSolver.hpp"
29 #include "Ipopt/BonIpoptWarmStart.hpp"
30 #ifdef COIN_HAS_FILTERSQP
31 #include "Filter/BonFilterSolver.hpp"
32 #include "Filter/BonFilterWarmStart.hpp"
33 //#include "Filter/BonBqpdWarmStart.hpp"
34 #endif
35
36 #include "OsiBranchingObject.hpp"
37 #include "OsiRowCutDebugger.hpp"
38 #include "BonStrongBranchingSolver.hpp"
39
40 //Macros to try and make messages definition less heavy
41 #include "BonMsgUtils.hpp"
42
43 using namespace Ipopt;
44
45
46 extern bool BonminAbortAll;
47 namespace Bonmin {
48 ///Register options
49 static void
register_general_options(SmartPtr<RegisteredOptions> roptions)50 register_general_options
51 (SmartPtr<RegisteredOptions> roptions)
52 {
53 roptions->SetRegisteringCategory("NLP interface", RegisteredOptions::BonminCategory);
54 roptions->AddStringOption3("nlp_solver",
55 "Choice of the solver for local optima of continuous NLP's",
56 "Ipopt",
57 "Ipopt", "Interior Point OPTimizer (https://projects.coin-or.org/Ipopt)",
58 "filterSQP", "Sequential quadratic programming trust region "
59 "algorithm (http://www-unix.mcs.anl.gov/~leyffer/solvers.html)",
60 "all", "run all available solvers at each node",
61 "Note that option will work only if the specified solver has been installed. Ipopt will usually be installed with Bonmin by default. For FilterSQP please see http://www-unix.mcs.anl.gov/~leyffer/solvers.html on how to obtain it and https://projects.coin-or.org/Bonmin/wiki/HintTricks on how to configure Bonmin to use it.");
62 roptions->setOptionExtraInfo("nlp_solver",127);
63
64 roptions->AddStringOption4("warm_start",
65 "Select the warm start method",
66 "none",
67 "none","No warm start, just start NLPs from optimal solution of the root relaxation",
68 "fake_basis", "builds fake basis, useful for cut management in Cbc (warm start is the same as in none)",
69 "optimum","Warm start with direct parent optimum",
70 "interior_point","Warm start with an interior point of direct parent",
71 "This will affect the function getWarmStart(), and as a consequence the warm starting in the various algorithms.");
72 roptions->setOptionExtraInfo("warm_start",8);
73
74 roptions->SetRegisteringCategory("Output and Loglevel", RegisteredOptions::BonminCategory);
75
76 roptions->AddBoundedIntegerOption("nlp_log_level",
77 "specify NLP solver interface log level (independent from ipopt print_level).",
78 0,2,1,
79 "Set the level of output of the OsiTMINLPInterface : "
80 "0 - none, 1 - normal, 2 - verbose"
81 );
82 roptions->setOptionExtraInfo("nlp_log_level",127);
83
84 roptions->AddStringOption2("file_solution",
85 "Write a file bonmin.sol with the solution",
86 "no",
87 "yes","",
88 "no","","");
89 roptions->setOptionExtraInfo("file_solution",127);
90
91 roptions->SetRegisteringCategory("NLP solution robustness", RegisteredOptions::BonminCategory);
92
93 roptions->AddLowerBoundedNumberOption("max_random_point_radius",
94 "Set max value r for coordinate of a random point.",
95 0.,1,1e5,
96 "When picking a random point, coordinate i will be in the interval [min(max(l,-r),u-r), max(min(u,r),l+r)] "
97 "(where l is the lower bound for the variable and u is its upper bound)");
98 roptions->setOptionExtraInfo("max_random_point_radius",8);
99
100 roptions->AddStringOption3("random_point_type","method to choose a random starting point",
101 "Jon",
102 "Jon", "Choose random point uniformly between the bounds",
103 "Andreas", "perturb the starting point of the problem within a prescribed interval",
104 "Claudia", "perturb the starting point using the perturbation radius suffix information",
105 "");
106 roptions->setOptionExtraInfo("random_point_type",8);
107
108 roptions->AddLowerBoundedNumberOption("random_point_perturbation_interval",
109 "Amount by which starting point is perturbed when choosing to pick random point by perturbing starting point",
110 0.,true, 1.,
111 "");
112 roptions->setOptionExtraInfo("random_point_perturbation_interval",8);
113
114
115 roptions->AddLowerBoundedIntegerOption
116 ("num_iterations_suspect",
117 "Number of iterations over which a node is considered \"suspect\" (for debugging purposes only, see detailed documentation).",
118 -1,-1,
119 "When the number of iterations to solve a node is above this number, the subproblem at this"
120 " node is considered to be suspect and it will be written into a file (set to -1 to deactivate this).");
121 roptions->setOptionExtraInfo("num_iterations_suspect",127);
122
123
124
125 roptions->AddLowerBoundedIntegerOption("num_retry_unsolved_random_point",
126 "Number $k$ of times that the algorithm will try to resolve an unsolved NLP with a random starting point "
127 "(we call unsolved an NLP for which Ipopt is not "
128 "able to guarantee optimality within the specified tolerances).",
129 0,0,
130 "When Ipopt fails to solve a continuous NLP sub-problem, if $k > 0$, the algorithm will "
131 "try again to solve the failed NLP with $k$ new randomly chosen starting points "
132 " or until the problem is solved with success.");
133 roptions->setOptionExtraInfo("num_retry_unsolved_random_point",127);
134
135 roptions->AddLowerBoundedNumberOption("resolve_on_small_infeasibility",
136 "If a locally infeasible problem is infeasible by less than this, resolve it "
137 "with initial starting point.",
138 0.,false, 0.,
139 "It is set to 0 by default with Ipopt. "
140 "When using FilterSQP, Bonmin sets it to a small value.");
141 roptions->setOptionExtraInfo("resolve_on_small_infeasibility",8);
142
143
144 roptions->SetRegisteringCategory("Nonconvex problems", RegisteredOptions::BonminCategory);
145
146
147 roptions->AddLowerBoundedIntegerOption("num_resolve_at_root",
148 "Number $k$ of tries to resolve the root node with different starting points.",
149 0,0,
150 "The algorithm will solve the root node with $k$ random starting points"
151 " and will keep the best local optimum found.");
152 roptions->setOptionExtraInfo("num_resolve_at_root",8);
153
154 roptions->AddLowerBoundedIntegerOption("num_resolve_at_node",
155 "Number $k$ of tries to resolve a node (other than the root) of the tree with different starting point.",
156 0,0,
157 "The algorithm will solve all the nodes with $k$ different random starting points "
158 "and will keep the best local optimum found.");
159 roptions->setOptionExtraInfo("num_resolve_at_node",8);
160
161 roptions->AddLowerBoundedIntegerOption("num_resolve_at_infeasibles",
162 "Number $k$ of tries to resolve an infeasible node (other than the root) of the tree with different starting point.",
163 0,0,
164 "The algorithm will solve all the infeasible nodes with $k$ different random starting points "
165 "and will keep the best local optimum found.");
166 roptions->setOptionExtraInfo("num_resolve_at_infeasibles",8);
167
168
169 roptions->AddStringOption2("dynamic_def_cutoff_decr",
170 "Do you want to define the parameter cutoff_decr dynamically?",
171 "no",
172 "no", "",
173 "yes", "");
174 roptions->setOptionExtraInfo("dynamic_def_cutoff_decr",8);
175
176 roptions->AddLowerBoundedNumberOption("coeff_var_threshold",
177 "Coefficient of variation threshold (for dynamic definition of cutoff_decr).",
178 0.0,
179 false,
180 0.1);
181 roptions->setOptionExtraInfo("coeff_var_threshold",8);
182
183 roptions->AddNumberOption("first_perc_for_cutoff_decr",
184 "The percentage used when, the coeff of variance is smaller than the threshold, to compute the cutoff_decr dynamically.",
185 -0.02);
186 roptions->setOptionExtraInfo("first_perc_for_cutoff_decr",8);
187
188 roptions->AddNumberOption("second_perc_for_cutoff_decr",
189 "The percentage used when, the coeff of variance is greater than the threshold, to compute the cutoff_decr dynamically.",
190 -0.05);
191 roptions->setOptionExtraInfo("second_perc_for_cutoff_decr",8);
192
193 }
194
register_OA_options(SmartPtr<RegisteredOptions> roptions)195 static void register_OA_options
196 (SmartPtr<RegisteredOptions> roptions)
197 {
198 //
199
200 roptions->SetRegisteringCategory("Outer Approximation strengthening", RegisteredOptions::UndocumentedCategory);
201 roptions->AddStringOption2("disjunctive_cut_type",
202 "Determine if and what kind of disjunctive cuts should be computed.",
203 "none",
204 "none", "No disjunctive cuts.",
205 "most-fractional", "If discrete variables present, compute disjunction for most-fractional variable");
206 roptions->setOptionExtraInfo("disjunctive_cut_type",119);
207
208
209
210
211 roptions->AddStringOption4("cut_strengthening_type",
212 "Determines if and what kind of cut strengthening should be performed.",
213 "none",
214 "none", "No strengthening of cuts.",
215 "sglobal", "Strengthen global cuts.",
216 "uglobal-slocal", "Unstrengthened global and strengthened local cuts",
217 "sglobal-slocal", "Strengthened global and strengthened local cuts",
218 "");
219 roptions->setOptionExtraInfo("cut_strengthening_type",119);
220
221 roptions->SetRegisteringCategory("Outer Approximation cuts generation", RegisteredOptions::BonminCategory);
222
223 roptions->AddStringOption2("oa_cuts_scope","Specify if OA cuts added are to be set globally or locally valid",
224 "global",
225 "local","Cuts are treated as locally valid",
226 "global", "Cuts are treated as globally valid",
227 "");
228 roptions->setOptionExtraInfo("oa_cuts_scope",119);
229
230 roptions->AddStringOption2("add_only_violated_oa","Do we add all OA cuts or only the ones violated by current point?",
231 "no",
232 "no","Add all cuts",
233 "yes","Add only violated cuts","");
234 roptions->setOptionExtraInfo("add_only_violated_oa",119);
235
236
237 roptions->AddLowerBoundedNumberOption("tiny_element","Value for tiny element in OA cut",
238 -0.,0,1e-08,
239 "We will remove \"cleanly\" (by relaxing cut) an element lower"
240 " than this.");
241 roptions->setOptionExtraInfo("tiny_element",119);
242
243 roptions->AddLowerBoundedNumberOption("very_tiny_element","Value for very tiny element in OA cut",
244 -0.,0,1e-17,
245 "Algorithm will take the risk of neglecting an element lower"
246 " than this.");
247 roptions->setOptionExtraInfo("very_tiny_element",119);
248
249 roptions->AddLowerBoundedNumberOption("oa_rhs_relax","Value by which to relax OA cut",
250 -0.,0,1e-8,
251 "RHS of OA constraints will be relaxed by this amount times the absolute value of the initial rhs if it is >= 1 (otherwise by this amount)."
252 );
253 roptions->setOptionExtraInfo("oa_rhs_relax",119);
254
255 roptions->SetRegisteringCategory("Output and Loglevel", RegisteredOptions::BonminCategory);
256
257 roptions->AddLowerBoundedIntegerOption("oa_cuts_log_level",
258 "level of log when generating OA cuts.",
259 0, 0,
260 "0: outputs nothing,\n"
261 "1: when a cut is generated, its violation and index of row from which it originates,\n"
262 "2: always output violation of the cut.\n"
263 "3: output generated cuts incidence vectors.");
264 roptions->setOptionExtraInfo("oa_cuts_log_level",119);
265
266 }
267
268
269 ///Register options
270 void
registerOptions(SmartPtr<RegisteredOptions> roptions)271 OsiTMINLPInterface::registerOptions
272 (SmartPtr<RegisteredOptions> roptions)
273 {
274 // We try to register the options - if those have been registered
275 // already, we catch the exception and don't need to do it again
276 try {
277 register_general_options(roptions);
278 register_OA_options(roptions);
279 #ifdef COIN_HAS_FILTERSQP
280 FilterSolver::RegisterOptions(roptions);
281 #endif
282 IpoptSolver::RegisterOptions(roptions);
283 }
284 catch(RegisteredOptions::OPTION_ALREADY_REGISTERED) {
285 // skipping
286 }
287
288
289 }
290
291 /** To set some application specific defaults. */
292 void
setAppDefaultOptions(SmartPtr<OptionsList> Options)293 OsiTMINLPInterface::setAppDefaultOptions(SmartPtr<OptionsList> Options)
294 {}
295
Messages()296 OsiTMINLPInterface::Messages::Messages
297 ():CoinMessages((int)OSITMINLPINTERFACE_DUMMY_END)
298 {
299 strcpy(source_ ,"NLP");
300 ADD_MSG(SOLUTION_FOUND, std_m, 2,
301 "After %d tries found a solution of %g (previous best %g).");
302
303 ADD_MSG(INFEASIBLE_SOLUTION_FOUND, std_m, 2,
304 "After %d tries found an solution of %g infeasible problem.");
305
306 ADD_MSG(UNSOLVED_PROBLEM_FOUND, std_m, 2,
307 "After %d tries found an solution of %g unsolved problem.");
308 ADD_MSG(WARN_SUCCESS_WS, warn_m, 2,
309 "Problem not solved with warm start but solved without");
310
311 ADD_MSG(WARNING_RESOLVING, warn_m,2,
312 "Trying to resolve NLP with different starting point (%d attempts).");
313 ADD_MSG(WARN_SUCCESS_RANDOM, warn_m, 1,
314 "Problem initially not solved but solved with a random starting point (success on %d attempt)");
315 ADD_MSG(WARN_CONTINUING_ON_FAILURE, warn_m, 1,
316 "Warning : continuing branching, while there are unrecovered failures at nodes");
317
318 ADD_MSG(SUSPECT_PROBLEM, std_m, 2, "NLP number %d is suspect (see bounds and start file)");
319 ADD_MSG(IPOPT_SUMMARY, std_m, 2, "Ipopt return (for %s): status %2d, iter count %4d, time %g");
320 ADD_MSG(BETTER_SOL, std_m, 2, "Solution of value %g found on %d'th attempt");
321
322 ADD_MSG(LOG_HEAD, std_m, 1,
323 "\n "
324 " Num Status Obj It time Location");
325 ADD_MSG(LOG_LINE, std_m, 1,
326 "%c %8d %11s %g %8d %g %20s");
327
328 ADD_MSG(ALTERNATE_OBJECTIVE, std_m, 1, "Objective value recomputed with alternate objective: %g.");
329
330 ADD_MSG(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, warn_m, 1,
331 "resolve called before any call to initialSol can not use warm starts.");
332 ADD_MSG(ERROR_NO_TNLPSOLVER, warn_m, 1,"Can not parse options when no IpApplication has been created");
333 ADD_MSG(WARNING_NON_CONVEX_OA, warn_m, 1,
334 "OA on non-convex constraint is very experimental.");
335 ADD_MSG(SOLVER_DISAGREE_STATUS, warn_m, 1, "%s says problem %s, %s says %s.");
336 ADD_MSG(SOLVER_DISAGREE_VALUE, warn_m, 1, "%s gives objective %.16g, %s gives %.16g.");
337
338 }
339
340
341 void
print(OsiRowCut & row)342 OsiTMINLPInterface::OaMessageHandler::print(OsiRowCut &row){
343 FILE * fp = filePointer();
344 const int &n = row.row().getNumElements();
345 fprintf(fp,"Row cut has %d elements. Lower bound: %g, upper bound %g.\n", n, row.lb(), row.ub());
346 const int * idx = row.row().getIndices();
347 const double * val = row.row().getElements();
348 for(int i = 0 ; i < n ; i++){
349 fprintf(fp,"%g, x%d",val[i], idx[i]);
350 if(i && i % 7 == 0)
351 fprintf(fp,"\n");
352 }
353 }
354
OaMessages()355 OsiTMINLPInterface::OaMessages::OaMessages(): CoinMessages((int) OA_MESSAGES_DUMMY_END){
356 strcpy(source_,"OaCg");
357 ADD_MSG(VIOLATED_OA_CUT_GENERATED, std_m, 1,"Row %d, cut violation is %g: Outer approximation cut generated.");
358
359 ADD_MSG(CUT_NOT_VIOLATED_ENOUGH, std_m, 2,"Row %d, cut violation is %g: Outer approximation cut not generated.");
360
361 ADD_MSG(OA_CUT_GENERATED, std_m, 1,"Row %d: Outer approximation cut not generated.");
362 }
363 bool OsiTMINLPInterface::hasPrintedOptions=0;
364
365 ////////////////////////////////////////////////////////////////////
366 // Constructors and desctructors //
367 ////////////////////////////////////////////////////////////////////
368 /// Default Constructor
OsiTMINLPInterface()369 OsiTMINLPInterface::OsiTMINLPInterface():
370 OsiSolverInterface(),
371 tminlp_(NULL),
372 problem_(NULL),
373 problem_to_optimize_(NULL),
374 feasibility_mode_(false),
375 app_(NULL),
376 debug_apps_(),
377 testOthers_(false),
378 warmstart_(NULL),
379 rowsense_(NULL),
380 rhs_(NULL),
381 rowrange_(NULL),
382 reducedCosts_(NULL),
383 OsiDualObjectiveLimit_(1e200),
384 hasVarNamesFile_(true),
385 nCallOptimizeTNLP_(0),
386 totalNlpSolveTime_(0),
387 totalIterations_(0),
388 maxRandomRadius_(1e08),
389 randomGenerationType_(0),
390 max_perturbation_(COIN_DBL_MAX),
391 pushValue_(1e-02),
392 numRetryInitial_(-1),
393 numRetryResolve_(-1),
394 numRetryInfeasibles_(-1),
395 numRetryUnsolved_(1),
396 infeasibility_epsilon_(0),
397 dynamicCutOff_(0),
398 coeff_var_threshold_(0.1),
399 first_perc_for_cutoff_decr_(-0.02),
400 second_perc_for_cutoff_decr_(-0.05),
401 messages_(),
402 pretendFailIsInfeasible_(0),
403 pretendSucceededNext_(0),
404 hasContinuedAfterNlpFailure_(false),
405 numIterationSuspect_(-1),
406 hasBeenOptimized_(false),
407 obj_(NULL),
408 feasibilityProblem_(NULL),
409 jRow_(NULL),
410 jCol_(NULL),
411 jValues_(NULL),
412 nnz_jac(0),
413 constTypes_(NULL),
414 nNonLinear_(0),
415 tiny_(1e-8),
416 veryTiny_(1e-20),
417 rhsRelax_(tiny_),
418 infty_(1e100),
419 warmStartMode_(None),
420 firstSolve_(true),
421 cutStrengthener_(NULL),
422 oaMessages_(),
423 oaHandler_(NULL),
424 newCutoffDecr(COIN_DBL_MAX)
425
426 {
427 oaHandler_ = new OaMessageHandler;
428 oaHandler_->setLogLevel(0);
429 }
430
431 void
initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,Ipopt::SmartPtr<Ipopt::OptionsList> options,Ipopt::SmartPtr<Ipopt::Journalist> journalist,const std::string & prefix,Ipopt::SmartPtr<TMINLP> tminlp)432 OsiTMINLPInterface::initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
433 Ipopt::SmartPtr<Ipopt::OptionsList> options,
434 Ipopt::SmartPtr<Ipopt::Journalist> journalist,
435 const std::string & prefix,
436 Ipopt::SmartPtr<TMINLP> tminlp){
437 if(!IsValid(app_))
438 createApplication(roptions, options, journalist, prefix);
439 setModel(tminlp);
440 }
441
setSolver(Ipopt::SmartPtr<TNLPSolver> app)442 void OsiTMINLPInterface::setSolver(Ipopt::SmartPtr<TNLPSolver> app){
443 app_ = app;
444 }
445
446
447 void
createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,Ipopt::SmartPtr<Ipopt::OptionsList> options,Ipopt::SmartPtr<Ipopt::Journalist> journalist,const std::string & prefix)448 OsiTMINLPInterface::createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
449 Ipopt::SmartPtr<Ipopt::OptionsList> options,
450 Ipopt::SmartPtr<Ipopt::Journalist> journalist,
451 const std::string & prefix
452 )
453 {
454 assert(!IsValid(app_));
455 int ival;
456 options->GetEnumValue("nlp_solver", ival, prefix.c_str());
457 Solver s = (Solver) ival;
458 if(s == EFilterSQP){
459 testOthers_ = false;;
460 #ifdef COIN_HAS_FILTERSQP
461 app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
462 #else
463 throw SimpleError("createApplication",
464 "Bonmin not configured to run with FilterSQP.");
465 #endif
466
467 debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix));
468 }
469 else if(s == EIpopt){
470 testOthers_ = false;
471 app_ = new IpoptSolver(roptions, options, journalist, prefix);
472
473 #ifdef COIN_HAS_FILTERSQP
474 debug_apps_.push_back(new Bonmin::FilterSolver(roptions, options, journalist, prefix));
475 #endif
476 }
477 else if(s == EAll){
478 #ifdef COIN_HAS_FILTERSQP
479 app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
480 #else
481 throw SimpleError("createApplication",
482 "Bonmin not configured to run with Ipopt.");
483 #endif
484 debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix));
485 testOthers_ = true;
486 }
487 if (!app_->Initialize("")) {
488 throw CoinError("Error during initialization of app_","createApplication", "OsiTMINLPInterface");
489 }
490 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin() ;
491 i != debug_apps_.end() ; i++){
492 (*i)->Initialize("");
493 }
494 extractInterfaceParams();
495
496 }
497
498 /** Facilitator to allocate a tminlp and an application. */
499 void
setModel(SmartPtr<TMINLP> tminlp)500 OsiTMINLPInterface::setModel(SmartPtr<TMINLP> tminlp)
501 {
502 assert(IsValid(tminlp));
503 tminlp_ = tminlp;
504 problem_ = new TMINLP2TNLP(tminlp_);
505 feasibilityProblem_ = new TNLP2FPNLP
506 (SmartPtr<TNLP>(GetRawPtr(problem_)));
507 if(feasibility_mode_){
508 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
509 }
510 else {
511 problem_to_optimize_ = GetRawPtr(problem_);
512 }
513 }
514
515
516
517 void
readOptionFile(const std::string & fileName)518 OsiTMINLPInterface::readOptionFile(const std::string & fileName)
519 {
520 if(IsValid(app_)){
521 std::ifstream is;
522 if (fileName != "") {
523 try {
524 is.open(fileName.c_str());
525 }
526 catch(std::bad_alloc) {
527 throw CoinError("Not enough memory to open option file.\n", "readOptionFile", "OsiTMINLPInterface");
528 }
529 }
530 options()->ReadFromStream(*app_->journalist(), is);
531 extractInterfaceParams();
532 }
533 }
534
535 /// Copy constructor
OsiTMINLPInterface(const OsiTMINLPInterface & source)536 OsiTMINLPInterface::OsiTMINLPInterface (const OsiTMINLPInterface &source):
537 OsiSolverInterface(source),
538 tminlp_(source.tminlp_),
539 problem_(NULL),
540 problem_to_optimize_(NULL),
541 feasibility_mode_(source.feasibility_mode_),
542 rowsense_(NULL),
543 rhs_(NULL),
544 rowrange_(NULL),
545 reducedCosts_(NULL),
546 OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
547 hasVarNamesFile_(source.hasVarNamesFile_),
548 nCallOptimizeTNLP_(0),
549 totalNlpSolveTime_(0),
550 totalIterations_(0),
551 maxRandomRadius_(source.maxRandomRadius_),
552 randomGenerationType_(source.randomGenerationType_),
553 max_perturbation_(source.max_perturbation_),
554 pushValue_(source.pushValue_),
555 numRetryInitial_(source.numRetryInitial_),
556 numRetryResolve_(source.numRetryResolve_),
557 numRetryInfeasibles_(source.numRetryInfeasibles_),
558 numRetryUnsolved_(source.numRetryUnsolved_),
559 infeasibility_epsilon_(source.infeasibility_epsilon_),
560 dynamicCutOff_(source.dynamicCutOff_),
561 coeff_var_threshold_(source.coeff_var_threshold_),
562 first_perc_for_cutoff_decr_(source.first_perc_for_cutoff_decr_),
563 second_perc_for_cutoff_decr_(source.second_perc_for_cutoff_decr_),
564 messages_(),
565 pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
566 pretendSucceededNext_(source.pretendSucceededNext_),
567 hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
568 numIterationSuspect_(source.numIterationSuspect_),
569 hasBeenOptimized_(source.hasBeenOptimized_),
570 obj_(NULL),
571 feasibilityProblem_(NULL),
572 jRow_(NULL),
573 jCol_(NULL),
574 jValues_(NULL),
575 nnz_jac(source.nnz_jac),
576 constTypes_(NULL),
577 // constTypesNum_(NULL),
578 nNonLinear_(0),
579 tiny_(source.tiny_),
580 veryTiny_(source.veryTiny_),
581 rhsRelax_(source.rhsRelax_),
582 infty_(source.infty_),
583 warmStartMode_(source.warmStartMode_),
584 firstSolve_(true),
585 cutStrengthener_(source.cutStrengthener_),
586 oaMessages_(),
587 oaHandler_(NULL),
588 newCutoffDecr(source.newCutoffDecr),
589 strong_branching_solver_(source.strong_branching_solver_)
590 {
591 if(IsValid(source.tminlp_)) {
592 problem_ = source.problem_->clone();
593 feasibilityProblem_ = new TNLP2FPNLP
594 (SmartPtr<TNLP>(GetRawPtr(problem_)), source.feasibilityProblem_);
595 if(feasibility_mode_)
596 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
597 else
598 problem_to_optimize_ = GetRawPtr(problem_);
599 pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
600 pretendSucceededNext_ = source.pretendSucceededNext_;
601
602 setAuxiliaryInfo(source.getAuxiliaryInfo());
603 // Copy options from old application
604 app_ = source.app_->clone();
605 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::const_iterator i = source.debug_apps_.begin();
606 i != source.debug_apps_.end() ; i++){
607 debug_apps_.push_back((*i)->clone());
608 }
609 testOthers_ = source.testOthers_;
610 }
611 else {
612 throw SimpleError("Don't know how to copy an empty IpoptInterface.",
613 "copy constructor");
614 }
615
616 warmstart_ = source.warmstart_ ? source.warmstart_->clone() : NULL;
617
618 if(source.obj_) {
619 obj_ = new double[source.getNumCols()];
620 CoinCopyN(source.obj_, source.getNumCols(), obj_);
621 }
622
623
624 oaHandler_ = new OaMessageHandler(*source.oaHandler_);;
625
626 }
627
628 OsiSolverInterface *
clone(bool copyData) const629 OsiTMINLPInterface::clone(bool copyData ) const
630 {
631 if(copyData)
632 return new OsiTMINLPInterface(*this);
633 else return new OsiTMINLPInterface;
634 }
635
636 /// Assignment operator
operator =(const OsiTMINLPInterface & rhs)637 OsiTMINLPInterface & OsiTMINLPInterface::operator=(const OsiTMINLPInterface& rhs)
638 {
639 if(this!= &rhs) {
640 OsiSolverInterface::operator=(rhs);
641 OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
642 nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
643 totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
644 totalIterations_ = rhs.totalIterations_;
645 maxRandomRadius_ = rhs.maxRandomRadius_;
646 hasVarNamesFile_ = rhs.hasVarNamesFile_;
647 pushValue_ = rhs.pushValue_;
648
649 delete warmstart_;
650 warmstart_ = NULL;
651
652 if(IsValid(rhs.tminlp_)) {
653
654 tminlp_ = rhs.tminlp_;
655 problem_ = new TMINLP2TNLP(tminlp_);
656 problem_to_optimize_ = GetRawPtr(problem_);
657 app_ = rhs.app_->clone();
658
659 warmstart_ = rhs.warmstart_ ? rhs.warmstart_->clone() : NULL;
660
661 feasibilityProblem_ = new TNLP2FPNLP
662 (SmartPtr<TNLP>(GetRawPtr(problem_)));
663 nnz_jac = rhs.nnz_jac;
664
665 if(constTypes_ != NULL) {
666 delete [] constTypes_;
667 constTypes_ = NULL;
668 }
669 if(rhs.constTypes_ != NULL) {
670 constTypes_ = new TNLP::LinearityType[getNumRows()];
671 CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
672 }
673 /*
674 if(constTypesNum_ != NULL) {
675 delete [] constTypesNum_;
676 constTypesNum_ = NULL;
677 }
678 if(rhs.constTypesNum_ != NULL) {
679 constTypesNum_ = new int[getNumRows()];
680 CoinCopyN(rhs.constTypesNum_, getNumRows(), constTypesNum_);
681 }
682 */
683 if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
684 jValues_ = new double [nnz_jac];
685 jCol_ = new Index [nnz_jac];
686 jRow_ = new Index [nnz_jac];
687 CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
688 CoinCopyN(rhs.jCol_ , nnz_jac,jCol_ );
689 CoinCopyN(rhs.jRow_ , nnz_jac,jRow_ );
690 }
691 else if(nnz_jac > 0) {
692 throw CoinError("Arrays for storing jacobian are inconsistant.",
693 "copy constructor",
694 "IpoptOAInterface");
695 }
696 tiny_ = rhs.tiny_;
697 veryTiny_ = rhs.veryTiny_;
698 rhsRelax_ = rhs.rhsRelax_;
699 infty_ = rhs.infty_;
700 warmStartMode_ = rhs.warmStartMode_;
701 newCutoffDecr = rhs.newCutoffDecr;
702
703 }
704 else {
705 tminlp_ =NULL;
706 problem_ = NULL;
707 app_ = NULL;
708 feasibilityProblem_ = NULL;
709 }
710
711
712 if(obj_) {
713 delete [] obj_;
714 obj_ = NULL;
715 }
716 if(rhs.obj_) {
717 obj_ = new double[rhs.getNumCols()];
718 CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
719 }
720
721 hasVarNamesFile_ = rhs.hasVarNamesFile_;
722
723 nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
724 totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
725 totalIterations_ = rhs.totalIterations_;
726 maxRandomRadius_ = rhs.maxRandomRadius_;
727 pushValue_ = rhs.pushValue_;
728 numRetryInitial_ = rhs.numRetryInitial_;
729 numRetryResolve_ = rhs.numRetryResolve_;
730 numRetryInfeasibles_ = rhs.numRetryInfeasibles_;
731 numRetryUnsolved_ = rhs.numRetryUnsolved_;
732 infeasibility_epsilon_ = rhs.infeasibility_epsilon_;
733 pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
734 pretendSucceededNext_ = rhs.pretendSucceededNext_;
735 numIterationSuspect_ = rhs.numIterationSuspect_;
736
737 hasBeenOptimized_ = rhs.hasBeenOptimized_;
738 cutStrengthener_ = rhs.cutStrengthener_;
739
740 delete oaHandler_;
741 oaHandler_ = new OaMessageHandler(*rhs.oaHandler_);
742 strong_branching_solver_ = rhs.strong_branching_solver_;
743
744 dynamicCutOff_ = rhs.dynamicCutOff_;
745 coeff_var_threshold_ = rhs.coeff_var_threshold_;
746 first_perc_for_cutoff_decr_ = rhs.first_perc_for_cutoff_decr_;
747 second_perc_for_cutoff_decr_ = rhs.second_perc_for_cutoff_decr_;
748
749 freeCachedData();
750 }
751 return *this;
752 }
753
options() const754 const SmartPtr<OptionsList> OsiTMINLPInterface::options() const
755 {
756 if(!IsValid(app_)) {
757 messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
758 return NULL;
759 }
760 else
761 return app_->options();
762 }
763
options()764 SmartPtr<OptionsList> OsiTMINLPInterface::options()
765 {
766 if(!IsValid(app_)) {
767 messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
768 return NULL;
769 }
770 else
771 return app_->options();
772 }
773
774 /// Destructor
~OsiTMINLPInterface()775 OsiTMINLPInterface::~OsiTMINLPInterface ()
776 {
777 freeCachedData();
778 delete [] jRow_;
779 delete [] jCol_;
780 delete [] jValues_;
781 delete [] constTypes_;
782 delete [] obj_;
783 delete oaHandler_;
784 delete warmstart_;
785 }
786
787 void
freeCachedColRim()788 OsiTMINLPInterface::freeCachedColRim()
789 {
790 if(reducedCosts_!=NULL) {
791 delete [] reducedCosts_;
792 reducedCosts_ = NULL;
793 }
794
795 }
796
797 void
freeCachedRowRim()798 OsiTMINLPInterface::freeCachedRowRim()
799 {
800 if(rowsense_!=NULL) {
801 delete [] rowsense_;
802 rowsense_ = NULL;
803 }
804 if(rhs_!=NULL) {
805 delete [] rhs_;
806 rhs_ = NULL;
807 }
808 if(rowrange_!=NULL) {
809 delete [] rowrange_;
810 rowrange_ = NULL;
811 }
812 // if(dualsol_!=NULL)
813 // {
814 // delete [] dualsol_; dualsol_ = NULL;
815 // }
816 }
817
818 void
freeCachedData()819 OsiTMINLPInterface::freeCachedData()
820 {
821 freeCachedColRim();
822 freeCachedRowRim();
823 }
824
825 const char * OsiTMINLPInterface::OPT_SYMB="OPT";
826 const char * OsiTMINLPInterface::FAILED_SYMB="FAILED";
827 const char * OsiTMINLPInterface::UNBOUND_SYMB="UNBOUNDED";
828 const char * OsiTMINLPInterface::INFEAS_SYMB="INFEAS";
829 const char * OsiTMINLPInterface::TIME_SYMB="TIME";
830 ///////////////////////////////////////////////////////////////////
831 // WarmStart Information //
832 ///////////////////////////////////////////////////////////////////
833
834
835 void
resolveForCost(int numsolve,bool keepWarmStart)836 OsiTMINLPInterface::resolveForCost(int numsolve, bool keepWarmStart)
837 {
838 // This method assumes that a problem has just been solved and we try for a
839 // different solution. So disregard (in fact, clear out) any warmstart info
840 // we might have, and acquire a new one before returning.
841 delete warmstart_;
842 warmstart_ = NULL;
843
844 double * of_current = (numsolve > 0) ? new double[numsolve]: NULL;
845 int num_failed, num_infeas;
846 double mean, std_dev, var_coeff;
847 double min = DBL_MAX;
848 double max = -DBL_MAX;
849
850 Coin::SmartPtr<SimpleReferencedPtr<CoinWarmStart> > ws_backup = NULL;
851 if(warmStartMode_ <= Optimum && keepWarmStart){
852 //if warm start is not exposed, we need to store the current starting point to
853 //restore it at the end of the method
854 ws_backup = make_referenced(app_->getUsedWarmStart(problem_));
855 }
856 /** Save current optimum. */
857 vector<double> point(getNumCols()*3+ getNumRows());
858 double bestBound = (isProvenOptimal()) ? getObjValue() : DBL_MAX;
859 CoinCopyN(getColSolution(),
860 getNumCols(), point());
861 CoinCopyN(getRowPrice(),
862 2*getNumCols()+ getNumRows(),
863 point() + getNumCols());
864 TNLPSolver::ReturnStatus savedStatus = optimizationStatus_;
865 if(isProvenOptimal())
866 messageHandler()->message(SOLUTION_FOUND,
867 messages_)
868 <<1<<getObjValue()<<bestBound
869 <<CoinMessageEol;
870 else
871 messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
872 messages_)
873 <<1
874 <<CoinMessageEol;
875
876 num_failed = 0;
877 num_infeas = 0;
878 mean = 0;
879
880 for(int f = 0; f < numsolve ; f++) {
881 messageHandler()->message(WARNING_RESOLVING,
882 messages_)
883 <<f+1<< CoinMessageEol ;
884 randomStartingPoint();
885 solveAndCheckErrors(0,0,"resolve cost");
886
887
888 char c=' ';
889 //Is solution better than previous
890 if(isProvenOptimal() &&
891 getObjValue()<bestBound) {
892 c='*';
893 messageHandler()->message(BETTER_SOL, messages_)<<getObjValue()<<f+1<< CoinMessageEol;
894 CoinCopyN(getColSolution(),
895 getNumCols(), point());
896 CoinCopyN(getRowPrice(),
897 2*getNumCols()+ getNumRows(),
898 point() + getNumCols());
899 bestBound = getObjValue();
900 savedStatus = optimizationStatus_;
901 }
902
903 messageHandler()->message(LOG_LINE, messages_)
904 <<c<<f+1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<app_->CPUTime()<<"resolve cost"<<CoinMessageEol;
905
906 if(isAbandoned()) {
907 num_failed++;
908 }
909 else if(isProvenPrimalInfeasible()) {
910 num_infeas++;
911 }
912
913 else if(isProvenOptimal())
914 messageHandler()->message(SOLUTION_FOUND,
915 messages_)
916 <<f+2<<getObjValue()<<bestBound
917 <<CoinMessageEol;
918 else if(!isAbandoned())
919 messageHandler()->message(UNSOLVED_PROBLEM_FOUND,
920 messages_)
921 <<f+2
922 <<CoinMessageEol;
923 else
924 messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
925 messages_)
926 <<f+2
927 <<CoinMessageEol;
928
929 if(of_current != NULL){
930 if(isProvenOptimal())
931 {
932 of_current[f] = getObjValue();
933 mean=mean+of_current[f];
934 if (of_current[f] < min)
935 min = of_current[f];
936 else if (of_current[f] > max)
937 max = of_current[f];
938 }
939 else
940 {
941 of_current[f] = 0;
942 }
943 }
944 }
945
946
947 if(of_current != NULL){
948 //calculate the mean
949 mean=mean/(numsolve-num_failed-num_infeas);
950
951 std_dev = 0;
952
953 //calculate the std deviation
954 for(int i=0; i<numsolve; i++)
955 {
956 if(of_current[i]!=0)
957 std_dev=std_dev+pow(of_current[i]-mean,2);
958 }
959 std_dev=pow((std_dev/(numsolve-num_failed-num_infeas)),0.5);
960
961 //calculate coeff of variation
962 var_coeff=std_dev/mean;
963
964 if(dynamicCutOff_)
965 {
966 if(var_coeff<0.1)
967 {
968 setNewCutoffDecr(mean*first_perc_for_cutoff_decr_);
969 }
970 else
971 {
972 setNewCutoffDecr(mean*second_perc_for_cutoff_decr_);
973 }
974 }
975 }
976
977 problem_->Set_x_sol(getNumCols(),point());
978 problem_->Set_dual_sol((int) point.size()-getNumCols(), point() + getNumCols());
979 problem_->set_obj_value(bestBound);
980 optimizationStatus_ = savedStatus;
981 hasBeenOptimized_ = true;
982
983 if(warmStartMode_ < Optimum && keepWarmStart) {
984 app_->setWarmStart(ws_backup->ptr(), problem_);
985 }
986 }
987
988 void
resolveForRobustness(int numsolve)989 OsiTMINLPInterface::resolveForRobustness(int numsolve)
990 {
991 // This method assumes that a problem has just been solved and we try for a
992 // different solution. So disregard (in fact, clear out) any warmstart info
993 // we might have, and acquire a new one before returning.
994 delete warmstart_;
995 warmstart_ = NULL;
996
997
998 CoinWarmStart * ws_backup = NULL;
999 if(warmStartMode_ < Optimum){
1000 //if warm start is not exposed, we need to store the current starting point to
1001 //restore it at the end of the method
1002 ws_backup = app_->getUsedWarmStart(problem_);
1003 }
1004 //std::cerr<<"Resolving the problem for robustness"<<std::endl;
1005 //First remove warm start point and resolve
1006 app_->disableWarmStart();
1007 problem()->resetStartingPoint();
1008 messageHandler()->message(WARNING_RESOLVING,
1009 messages_)
1010 <<1<< CoinMessageEol ;
1011 solveAndCheckErrors(0,0,"resolve robustness");
1012
1013
1014 char c='*';
1015 if(isAbandoned()) {
1016 c=' ';
1017 }
1018 messageHandler()->message(LOG_LINE, messages_)
1019 <<c<<1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<
1020 app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
1021
1022
1023 if(!isAbandoned()) {
1024 messageHandler()->message(WARN_SUCCESS_WS, messages_) << CoinMessageEol ;
1025 // re-enable warmstart and get it
1026 app_->enableWarmStart();
1027 if (warmStartMode_ < Optimum) {
1028 app_->setWarmStart(ws_backup, problem_);
1029 delete ws_backup;
1030 }
1031 return; //we won go on
1032 }
1033
1034 //still unsolved try again with different random starting points
1035 for(int f = 0; f < numsolve ; f++) {
1036 messageHandler()->message(WARNING_RESOLVING,
1037 messages_)
1038 <<f+2<< CoinMessageEol ;
1039
1040 randomStartingPoint();
1041 solveAndCheckErrors(0,0,"resolve robustness");
1042
1043
1044 messageHandler()->message(IPOPT_SUMMARY, messages_)
1045 <<"resolveForRobustness"<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
1046
1047
1048 char c='*';
1049 if(isAbandoned()) {
1050 c=' ';
1051 }
1052 messageHandler()->message(LOG_LINE, messages_)
1053 <<c<<f+2<<statusAsString()<<getObjValue()
1054 <<app_->IterationCount()<<app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
1055
1056
1057 if(!isAbandoned()) {
1058 messageHandler()->message(WARN_SUCCESS_RANDOM, messages_)
1059 << f+2 << CoinMessageEol ;
1060 // re-enable warmstart and get it
1061 app_->enableWarmStart();
1062 if (warmStartMode_ < Optimum) {
1063 app_->setWarmStart(ws_backup, problem_);
1064 delete ws_backup;
1065 }
1066
1067 return; //we have found a solution and continue
1068 }
1069 }
1070
1071
1072 if(warmStartMode_ < Optimum){
1073 app_->setWarmStart(ws_backup, problem_);
1074 delete ws_backup;
1075 }
1076 if(pretendFailIsInfeasible_) {
1077 pretendSucceededNext_ = true;
1078 if(pretendFailIsInfeasible_ == 1) {
1079 messageHandler()->message(WARN_CONTINUING_ON_FAILURE,
1080 messages_)
1081 <<CoinMessageEol;
1082 hasContinuedAfterNlpFailure_ = 1;
1083 }
1084 return;
1085 }
1086 else {
1087 std::string probName;
1088 getStrParam(OsiProbName,probName);
1089 throw newUnsolvedError(app_->errorCode(), problem_,
1090 probName);
1091 }
1092 // Do NOT get warmstart in other cases
1093 }
1094
1095 ////////////////////////////////////////////////////////////////////
1096 // Problem information methods //
1097 ////////////////////////////////////////////////////////////////////
1098 /// Get number of columns
getNumCols() const1099 int OsiTMINLPInterface::getNumCols() const
1100 {
1101
1102 return problem_->num_variables();
1103 }
1104
1105
1106 /// Get number of rows
1107 int
getNumRows() const1108 OsiTMINLPInterface::getNumRows() const
1109 {
1110 return problem_->num_constraints();
1111 }
1112
1113 const double *
getColLower() const1114 OsiTMINLPInterface::getColLower() const
1115 {
1116 return problem_->x_l();
1117 }
1118
1119 const double *
getColUpper() const1120 OsiTMINLPInterface::getColUpper() const
1121 {
1122 return problem_->x_u();
1123 }
1124
1125 #if 1
1126
1127
1128 ///get name of variables
1129 const OsiSolverInterface::OsiNameVec&
getVarNames()1130 OsiTMINLPInterface::getVarNames() {
1131 return getColNames();
1132 }
1133 #endif
1134
1135
extractSenseRhsAndRange() const1136 void OsiTMINLPInterface::extractSenseRhsAndRange() const
1137 {
1138 assert(rowsense_==NULL&&rhs_==NULL&&rowrange_==NULL);
1139 int numrows = problem_->num_constraints();
1140 if(numrows == 0) return;
1141 const double * rowLower = getRowLower();
1142 const double * rowUpper = getRowUpper();
1143 rowsense_ = new char [numrows];
1144 rhs_ = new double [numrows];
1145 rowrange_ = new double [numrows];
1146 for(int i = 0 ; i < numrows ; i++) {
1147 rowrange_[i]=0.;
1148 convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
1149 }
1150 }
1151
1152 /** Get pointer to array[getNumRows()] of row constraint senses.
1153 <ul>
1154 <li>'L': <= constraint
1155 <li>'E': = constraint
1156 <li>'G': >= constraint
1157 <li>'R': ranged constraint
1158 <li>'N': free constraint
1159 </ul>
1160 */
1161 const char *
getRowSense() const1162 OsiTMINLPInterface::getRowSense() const
1163 {
1164 if(rowsense_==NULL) {
1165 extractSenseRhsAndRange();
1166 }
1167 return rowsense_;
1168 }
1169
1170 /** Get pointer to array[getNumRows()] of rows right-hand sides
1171 <ul>
1172 <li> if rowsense()[i] == 'L' then rhs()[i] == rowupper()[i]
1173 <li> if rowsense()[i] == 'G' then rhs()[i] == rowlower()[i]
1174 <li> if rowsense()[i] == 'R' then rhs()[i] == rowupper()[i]
1175 <li> if rowsense()[i] == 'N' then rhs()[i] == 0.0
1176 </ul>
1177 */
1178 const double *
getRightHandSide() const1179 OsiTMINLPInterface::getRightHandSide() const
1180 {
1181 if(rhs_==NULL) {
1182 extractSenseRhsAndRange();
1183 }
1184 return rhs_;
1185 }
1186
1187 /** Get pointer to array[getNumRows()] of row ranges.
1188 <ul>
1189 <li> if rowsense()[i] == 'R' then
1190 rowrange()[i] == rowupper()[i] - rowlower()[i]
1191 <li> if rowsense()[i] != 'R' then
1192 rowrange()[i] is 0.0
1193 </ul>
1194 */
1195 const double *
getRowRange() const1196 OsiTMINLPInterface::getRowRange() const
1197 {
1198 if(rowrange_==NULL) {
1199 extractSenseRhsAndRange();
1200 }
1201 return rowrange_;
1202 }
1203
1204 const double *
getRowLower() const1205 OsiTMINLPInterface::getRowLower() const
1206 {
1207 return problem_->g_l();
1208 }
1209
1210 const double *
getRowUpper() const1211 OsiTMINLPInterface::getRowUpper() const
1212 {
1213 return problem_->g_u();
1214 }
1215
1216 /// Return true if column is continuous
1217 bool
isContinuous(int colNumber) const1218 OsiTMINLPInterface::isContinuous(int colNumber) const
1219 {
1220 return (problem_->var_types()[colNumber]==TMINLP::CONTINUOUS);
1221 }
1222
1223 /// Return true if column is binary
1224 bool
isBinary(int colNumber) const1225 OsiTMINLPInterface::isBinary(int colNumber) const
1226 {
1227 return (problem_->var_types()[colNumber]==TMINLP::BINARY);
1228 }
1229
1230 /** Return true if column is integer.
1231 Note: This function returns true if the the column
1232 is binary or a general integer.
1233 */
1234 bool
isInteger(int colNumber) const1235 OsiTMINLPInterface::isInteger(int colNumber) const
1236 {
1237 return ((problem_->var_types()[colNumber]==TMINLP::BINARY)||
1238 (problem_->var_types()[colNumber]==TMINLP::INTEGER));
1239 }
1240
1241 /// Return true if column is general integer
1242 bool
isIntegerNonBinary(int colNumber) const1243 OsiTMINLPInterface::isIntegerNonBinary(int colNumber) const
1244 {
1245 return (problem_->var_types()[colNumber]==TMINLP::INTEGER);
1246 }
1247 /// Return true if column is binary and not fixed at either bound
1248 bool
isFreeBinary(int colNumber) const1249 OsiTMINLPInterface::isFreeBinary(int colNumber) const
1250 {
1251 return ((problem_->var_types()[colNumber]==TMINLP::BINARY)
1252 &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
1253 }
1254
1255 /// Get solver's value for infinity
1256 double
getInfinity() const1257 OsiTMINLPInterface::getInfinity() const
1258 {
1259 return COIN_DBL_MAX;
1260 }
1261
1262 /// Get pointer to array[getNumCols()] of primal solution vector
1263 const double *
getColSolution() const1264 OsiTMINLPInterface::getColSolution() const
1265 {
1266 if(hasBeenOptimized_)
1267 return problem_->x_sol();
1268 else
1269 return problem_->x_init();
1270 }
1271
1272 /// Get pointer to array[getNumRows()] of dual prices
1273 const double *
getRowPrice() const1274 OsiTMINLPInterface::getRowPrice() const
1275 {
1276 if(hasBeenOptimized_)
1277 return problem_->duals_sol();
1278 else
1279 return problem_->duals_init();
1280 }
1281
1282 /// Get a pointer to array[getNumCols()] of reduced costs
1283 const double *
getReducedCost() const1284 OsiTMINLPInterface::getReducedCost() const
1285 {
1286 (*handler_)<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<CoinMessageEol;
1287 if(reducedCosts_==NULL) {
1288 reducedCosts_ = new double [getNumCols()];
1289 CoinFillN(reducedCosts_,getNumCols(),0.);
1290 }
1291 return reducedCosts_;
1292 }
1293
1294 /** Get pointer to array[getNumRows()] of row activity levels (constraint
1295 matrix times the solution vector */
1296 const double *
getRowActivity() const1297 OsiTMINLPInterface::getRowActivity() const
1298 {
1299 return problem_->g_sol();
1300 }
1301
1302 /** Get how many iterations it took to solve the problem (whatever
1303 "iteration" mean to the solver.
1304 */
1305 int
getIterationCount() const1306 OsiTMINLPInterface::getIterationCount() const
1307 {
1308 return app_->IterationCount();
1309 }
1310
1311
1312 /** Set a single column lower bound.
1313 Use -getInfinity() for -infinity. */
1314 void
setColLower(int elementIndex,double elementValue)1315 OsiTMINLPInterface::setColLower( int elementIndex, double elementValue )
1316 {
1317 // if(fabs(problem_->x_l()[elementIndex]-elementValue)>1e-06)
1318 problem_->SetVariableLowerBound(elementIndex,elementValue);
1319 hasBeenOptimized_ = false;
1320 }
1321
1322 /** Set a single column upper bound.
1323 Use getInfinity() for infinity. */
1324 void
setColUpper(int elementIndex,double elementValue)1325 OsiTMINLPInterface::setColUpper( int elementIndex, double elementValue )
1326 {
1327 // if(fabs(problem_->x_u()[elementIndex]-elementValue)>1e-06)
1328 problem_->SetVariableUpperBound(elementIndex,elementValue);
1329 hasBeenOptimized_ = false;
1330 }
1331
1332 /** Set the lower bounds for all columns
1333 Use -getInfinity() for -infinity. */
1334 void
setColLower(const double * array)1335 OsiTMINLPInterface::setColLower( const double* array )
1336 {
1337 problem_->SetVariablesLowerBounds(problem_->num_variables(),
1338 array);
1339 hasBeenOptimized_ = false;
1340 }
1341
1342 /** Set Set the upper bounds for all columns
1343 Use getInfinity() for infinity. */
1344 void
setColUpper(const double * array)1345 OsiTMINLPInterface::setColUpper( const double* array )
1346 {
1347 problem_->SetVariablesUpperBounds(problem_->num_variables(),
1348 array);
1349 hasBeenOptimized_ = false;
1350 }
1351
1352 /** Set a single row lower bound.
1353 Use -getInfinity() for -infinity. */
1354 void
setRowLower(int elementIndex,double elementValue)1355 OsiTMINLPInterface::setRowLower( int elementIndex, double elementValue )
1356 {
1357 throw SimpleError("Not implemented yet but should be if necessary.",
1358 "setRowLower");
1359 hasBeenOptimized_ = false;
1360 }
1361
1362 /** Set a single row upper bound.
1363 Use getInfinity() for infinity. */
1364 void
setRowUpper(int elementIndex,double elementValue)1365 OsiTMINLPInterface::setRowUpper( int elementIndex, double elementValue )
1366 {
1367 throw SimpleError("Not implemented yet but should be if necessary.",
1368 "setRowUpper");
1369 hasBeenOptimized_ = false;
1370 }
1371
1372 /** Set the type of a single row */
1373 void
setRowType(int index,char sense,double rightHandSide,double range)1374 OsiTMINLPInterface::setRowType(int index, char sense, double rightHandSide,
1375 double range)
1376 {
1377 throw SimpleError("Not implemented yet but should be if necessary.",
1378 "setRowType");
1379 hasBeenOptimized_ = false;
1380 }
1381
1382
1383 /// Set the objective function sense.
1384 /// (1 for min (default), -1 for max)
1385 void
setObjSense(double s)1386 OsiTMINLPInterface::setObjSense(double s)
1387 {
1388 throw SimpleError("Can not change objective sense of an Ipopt problem.",
1389 "setObjSense");
1390 hasBeenOptimized_ = false;
1391 }
1392
1393 /** Set the primal solution variable values
1394
1395 colsol[getNumCols()] is an array of values for the primal variables.
1396 These values are copied to memory owned by the solver interface object
1397 or the solver. They will be returned as the result of getColSolution()
1398 until changed by another call to setColSolution() or by a call to any
1399 solver routine. Whether the solver makes use of the solution in any
1400 way is solver-dependent.
1401 */
1402 void
setColSolution(const double * colsol)1403 OsiTMINLPInterface::setColSolution(const double *colsol)
1404 {
1405 if(colsol != NULL)
1406 problem_->setxInit(getNumCols(), colsol);
1407 else
1408 problem_->resetStartingPoint();
1409 hasBeenOptimized_ = false;
1410 }
1411
1412 /** Set dual solution variable values
1413
1414 rowprice[getNumRows()] is an array of values for the dual
1415 variables. These values are copied to memory owned by the solver
1416 interface object or the solver. They will be returned as the result of
1417 getRowPrice() until changed by another call to setRowPrice() or by a
1418 call to any solver routine. Whether the solver makes use of the
1419 solution in any way is solver-dependent.
1420 */
1421
1422 void
setRowPrice(const double * rowprice)1423 OsiTMINLPInterface::setRowPrice(const double * rowprice)
1424 {
1425 problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
1426 hasBeenOptimized_ = false;
1427 }
1428
1429 /*! \brief Get an empty warm start object
1430
1431 This routine returns an empty CoinWarmStartBasis object. Its purpose is
1432 to provide a way to give a client a warm start basis object of the
1433 appropriate type, which can resized and modified as desired.
1434 */
1435 CoinWarmStart *
getEmptyWarmStart() const1436 OsiTMINLPInterface::getEmptyWarmStart () const
1437 {return app_->getEmptyWarmStart();}
1438
1439 /** Get warmstarting information */
1440 CoinWarmStart*
getWarmStart() const1441 OsiTMINLPInterface::getWarmStart() const
1442 {
1443 if (warmStartMode_ >= Optimum) {
1444 return internal_getWarmStart();;
1445 }
1446 else if(warmStartMode_ == FakeBasis){
1447 return build_fake_basis();
1448 }
1449 else {
1450 return getEmptyWarmStart();
1451 }
1452 }
1453 /** Set warmstarting information. Return true/false depending on whether
1454 the warmstart information was accepted or not. */
1455 bool
setWarmStart(const CoinWarmStart * ws)1456 OsiTMINLPInterface::setWarmStart(const CoinWarmStart* ws)
1457 {
1458 if (warmStartMode_ >= Optimum) {
1459 return internal_setWarmStart(ws);
1460 }
1461 else {
1462 return true;
1463 }
1464 }
1465
1466
1467 #define EPSILON 1e-4
1468
1469 CoinWarmStart*
build_fake_basis() const1470 OsiTMINLPInterface::build_fake_basis() const{
1471 CoinWarmStartBasis * r_val = new CoinWarmStartBasis();
1472 int n_con = problem_->num_constraints();
1473 r_val->setSize(problem_->num_variables(), n_con);
1474 const double * act = problem_->g_sol();
1475 const double * lb = problem_->g_l();
1476 const double * ub = problem_->g_u();
1477 for(int i = 0 ; i < n_con ; i++){
1478 if(lb[i] > ub[i] - EPSILON){
1479 r_val->setArtifStatus(i, CoinWarmStartBasis::isFree);
1480 }
1481 if(act[i] > ub[i] - EPSILON){
1482 r_val->setArtifStatus(i, CoinWarmStartBasis::atLowerBound);
1483 }
1484 else if(act[i] < lb[i] + EPSILON){
1485 r_val->setArtifStatus(i, CoinWarmStartBasis::atLowerBound);
1486 }
1487 else {
1488 r_val->setArtifStatus(i, CoinWarmStartBasis::basic);
1489 }
1490 }
1491 return r_val;
1492 }
1493
1494 /** Get warmstarting information */
1495 CoinWarmStart*
internal_getWarmStart() const1496 OsiTMINLPInterface::internal_getWarmStart() const
1497 {
1498 if (warmStartMode_ >= Optimum && warmstart_) {
1499 return warmstart_->clone();
1500 }
1501 else {
1502 return getEmptyWarmStart();
1503 }
1504 }
1505 /** Set warmstarting information. Return true/false depending on whether
1506 the warmstart information was accepted or not. */
1507 bool
internal_setWarmStart(const CoinWarmStart * ws)1508 OsiTMINLPInterface::internal_setWarmStart(const CoinWarmStart* ws)
1509 {
1510 delete warmstart_;
1511 warmstart_ = NULL;
1512 hasBeenOptimized_ = false;
1513 if (warmStartMode_ >= Optimum) {
1514 if (ws == NULL) {
1515 return true;
1516 }
1517 if(app_->warmStartIsValid(ws)) {
1518 warmstart_ = ws->clone();
1519 return true;
1520 }
1521 // See if it is anything else than the CoinWarmStartBasis that all others
1522 // derive from
1523 const CoinWarmStartPrimalDual* pdws =
1524 dynamic_cast<const CoinWarmStartPrimalDual*>(ws);
1525 if (pdws) {
1526 // Must be an IpoptWarmStart, since the others do not derive from this.
1527 // Accept it.
1528 warmstart_ = new IpoptWarmStart(*pdws);
1529 return true;
1530 }
1531 return false;
1532 }
1533 else {
1534 return true;
1535 }
1536 }
1537
1538 /** Set the index-th variable to be a continuous variable */
1539 void
setContinuous(int index)1540 OsiTMINLPInterface::setContinuous(int index)
1541 {
1542 problem_->SetVariableType(index, TMINLP::CONTINUOUS);
1543 hasBeenOptimized_ = false;
1544 }
1545 /** Set the index-th variable to be an integer variable */
1546 void
setInteger(int index)1547 OsiTMINLPInterface::setInteger(int index)
1548 {
1549 problem_->SetVariableType(index, TMINLP::INTEGER);
1550 hasBeenOptimized_ = false;
1551 }
1552
1553 /// Get objective function value (can't use default)
1554 double
getObjValue() const1555 OsiTMINLPInterface::getObjValue() const
1556 {
1557 return problem_->obj_value();
1558 }
1559
1560 //#############################################################################
1561 // Parameter related methods
1562 //#############################################################################
1563
1564 bool
setIntParam(OsiIntParam key,int value)1565 OsiTMINLPInterface::setIntParam(OsiIntParam key, int value)
1566 {
1567 // debugMessage("OsiCpxSolverInterface::setIntParam(%d, %d)\n", key, value);
1568
1569 bool retval = false;
1570 switch (key) {
1571 case OsiMaxNumIteration:
1572 retval = false;
1573 break;
1574 case OsiMaxNumIterationHotStart:
1575 if( value >= 0 ) {
1576 retval = false;
1577 }
1578 else
1579 retval = false;
1580 break;
1581 case OsiLastIntParam:
1582 retval = false;
1583 break;
1584 default:
1585 retval = false;
1586 (*handler_)<< "Unhandled case in setIntParam\n"<<CoinMessageEol;
1587 break;
1588 }
1589 return retval;
1590 }
1591
1592 //-----------------------------------------------------------------------------
1593
1594 bool
setDblParam(OsiDblParam key,double value)1595 OsiTMINLPInterface::setDblParam(OsiDblParam key, double value)
1596 {
1597 // debugMessage("OsiTMINLPInterface::setDblParam(%d, %g)\n", key, value);
1598
1599 bool retval = false;
1600 switch (key) {
1601 case OsiDualObjectiveLimit:
1602 OsiDualObjectiveLimit_ = value;
1603 retval = true;
1604 break;
1605 case OsiPrimalObjectiveLimit:
1606 (*handler_)<<"Can not set primal objective limit parameter"<<CoinMessageEol;
1607 retval = false;
1608 break;
1609 case OsiDualTolerance:
1610 (*handler_)<<"Can not set dual tolerance parameter"<<CoinMessageEol;
1611 retval = false;
1612 break;
1613 case OsiPrimalTolerance:
1614 (*handler_)<<"Can not set primal tolerance parameter"<<CoinMessageEol;
1615 retval = false;
1616 case OsiObjOffset:
1617 retval = OsiSolverInterface::setDblParam(key,value);
1618 break;
1619 case OsiLastDblParam:
1620 retval = false;
1621 break;
1622 default:
1623 retval = false;
1624 (*handler_) << "Unhandled case in setDblParam"<<CoinMessageEol;
1625 break;
1626 }
1627 return retval;
1628 }
1629
1630
1631 //-----------------------------------------------------------------------------
1632
1633 bool
setStrParam(OsiStrParam key,const std::string & value)1634 OsiTMINLPInterface::setStrParam(OsiStrParam key, const std::string & value)
1635 {
1636 // debugMessage("OsiTMINLPInterface::setStrParam(%d, %s)\n", key, value.c_str());
1637
1638 bool retval=false;
1639 switch (key) {
1640 case OsiProbName:
1641 OsiSolverInterface::setStrParam(key,value);
1642 return retval = true;
1643 case OsiSolverName:
1644 return false;
1645 case OsiLastStrParam:
1646 return false;
1647 }
1648 return false;
1649 }
1650
1651 //-----------------------------------------------------------------------------
1652
1653 bool
getIntParam(OsiIntParam key,int & value) const1654 OsiTMINLPInterface::getIntParam(OsiIntParam key, int& value) const
1655 {
1656 // debugMessage("OsiTMINLPInterface::getIntParam(%d)\n", key);
1657
1658 value = -COIN_INT_MAX; // Give a dummy value
1659 bool retval = false;
1660 switch (key) {
1661 case OsiMaxNumIteration:
1662 retval = false;
1663 break;
1664 case OsiMaxNumIterationHotStart:
1665 retval = false;
1666 break;
1667 case OsiLastIntParam:
1668 retval = false;
1669 break;
1670 default:
1671 retval = false;
1672 (*handler_) << "Unhandled case in setIntParam"<<CoinMessageEol;
1673 }
1674 return retval;
1675 }
1676
1677 //-----------------------------------------------------------------------------
1678
1679 bool
getDblParam(OsiDblParam key,double & value) const1680 OsiTMINLPInterface::getDblParam(OsiDblParam key, double& value) const
1681 {
1682 // debugMessage("OsiTMINLPInterface::getDblParam(%d)\n", key);
1683
1684 bool retval = false;
1685 switch (key) {
1686 case OsiDualObjectiveLimit:
1687 value = OsiDualObjectiveLimit_;
1688 retval = true;
1689 break;
1690 case OsiPrimalObjectiveLimit:
1691 value = getInfinity();
1692 retval = true;
1693 break;
1694 case OsiDualTolerance:
1695 retval = false;
1696 break;
1697 case OsiPrimalTolerance:
1698 options()->GetNumericValue("tol", value,"");
1699 value = 1e-07;
1700 retval = true;
1701 break;
1702 case OsiObjOffset:
1703 retval = OsiSolverInterface::getDblParam(key, value);
1704 break;
1705 case OsiLastDblParam:
1706 retval = false;
1707 break;
1708 }
1709 return retval;
1710 }
1711
1712
1713 //-----------------------------------------------------------------------------
1714
1715 bool
getStrParam(OsiStrParam key,std::string & value) const1716 OsiTMINLPInterface::getStrParam(OsiStrParam key, std::string & value) const
1717 {
1718 // debugMessage("OsiTMINLPInterface::getStrParam(%d)\n", key);
1719
1720 switch (key) {
1721 case OsiProbName:
1722 OsiSolverInterface::getStrParam(key, value);
1723 break;
1724 case OsiSolverName:
1725 value = "Ipopt";
1726 break;
1727 case OsiLastStrParam:
1728 return false;
1729 }
1730
1731 return true;
1732 }
1733
1734 void
set_linearizer(Ipopt::SmartPtr<TMINLP2OsiLP> linearizer)1735 OsiTMINLPInterface::set_linearizer(Ipopt::SmartPtr<TMINLP2OsiLP> linearizer)
1736 {
1737 linearizer_ = linearizer->clone();
1738 linearizer_->set_tols(tiny_, veryTiny_, rhsRelax_, infty_);
1739 linearizer_->set_model(GetRawPtr(problem_));
1740 }
1741
1742 Ipopt::SmartPtr<TMINLP2OsiLP>
linearizer()1743 OsiTMINLPInterface::linearizer(){
1744 return linearizer_;
1745 }
1746
1747
1748 void
randomStartingPoint()1749 OsiTMINLPInterface::randomStartingPoint()
1750 {
1751 int numcols = getNumCols();
1752 const double * colLower = getColLower();
1753 const double * colUpper = getColUpper();
1754 double * sol = new double[numcols];
1755 const Number * x_init = problem_->x_init_user();
1756 const double* perturb_radius = NULL;
1757 if (randomGenerationType_ == perturb_suffix) {
1758 const TMINLP::PerturbInfo* pertubinfo = tminlp_->perturbInfo();
1759 if (pertubinfo) {
1760 perturb_radius = pertubinfo->GetPerturbationArray();
1761 }
1762 if (!perturb_radius) {
1763 throw SimpleError("Can't use perturb_radius if no radii are given.",
1764 "randomStartingPoint");
1765 }
1766 }
1767 for(int i = 0 ; i < numcols ; i++) {
1768 if(randomGenerationType_ == uniform || x_init[i] < colLower[i] || x_init[i] > colUpper[i]) {
1769 double lower = std::min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
1770 lower = std::max(colLower[i], lower);
1771 double upper = std::max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
1772 upper = std::min(colUpper[i],upper);
1773 lower = std::min(upper,lower);
1774 upper = std::max(upper, lower);
1775 double interval = upper - lower;
1776 sol[i] = CoinDrand48()*(interval) + lower;
1777 }
1778 else if (randomGenerationType_ == perturb) {
1779 const double lower = std::max(x_init[i] - max_perturbation_, colLower[i]);
1780 const double upper = std::min(x_init[i] + max_perturbation_, colUpper[i]);
1781 const double interval = upper - lower;
1782 sol[i] = lower + CoinDrand48()*(interval);
1783 }
1784 else if (randomGenerationType_ == perturb_suffix) {
1785 const double radius = perturb_radius[i];
1786 const double lower = std::max(x_init[i] - radius*max_perturbation_, colLower[i]);
1787 const double upper = std::min(x_init[i] + radius*max_perturbation_, colUpper[i]);
1788 const double interval = upper - lower;
1789 sol[i] = lower + CoinDrand48()*(interval);
1790 }
1791 }
1792 app_->disableWarmStart();
1793 setColSolution(sol);
1794 delete [] sol;
1795 }
1796
1797
1798
1799 /** This methods initialiaze arrays for storing the jacobian */
initializeJacobianArrays()1800 int OsiTMINLPInterface::initializeJacobianArrays()
1801 {
1802 Index n, m, nnz_h_lag;
1803 TNLP::IndexStyleEnum index_style;
1804 problem_to_optimize_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
1805
1806 if(jRow_ != NULL) delete [] jRow_;
1807 if(jCol_ != NULL) delete [] jCol_;
1808 if(jValues_ != NULL) delete [] jValues_;
1809
1810 jRow_ = new Index[nnz_jac];
1811 jCol_ = new Index[nnz_jac];
1812 jValues_ = new Number[nnz_jac];
1813 problem_to_optimize_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
1814 if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
1815 {
1816 for(int i = 0 ; i < nnz_jac ; i++){
1817 jRow_[i]--;
1818 jCol_[i]--;
1819 }
1820 }
1821
1822 if(constTypes_ != NULL) delete [] constTypes_;
1823 // if(constTypesNum_ != NULL) delete [] constTypesNum_;
1824
1825 constTypes_ = new TNLP::LinearityType[getNumRows()];
1826 problem_to_optimize_->get_constraints_linearity(getNumRows(), constTypes_);
1827 // constTypesNum_ = new int[getNumRows()];
1828 for(int i = 0; i < getNumRows() ; i++) {
1829 if(constTypes_[i]==TNLP::NON_LINEAR) {
1830 //constTypesNum_[i] =
1831 nNonLinear_++;
1832 }
1833 }
1834 return nnz_jac;
1835 }
1836
1837
1838 double
getConstraintsViolation(const double * x,double & obj)1839 OsiTMINLPInterface::getConstraintsViolation(const double *x, double &obj)
1840 {
1841 int numcols = getNumCols();
1842 int numrows = getNumRows();
1843 double * g = new double[numrows];
1844 tminlp_->eval_g(numcols, x, 1, numrows, g);
1845 const double * rowLower = getRowLower();
1846 const double * rowUpper = getRowUpper();
1847
1848
1849 double norm = 0;
1850 for(int i = 0; i< numrows ; i++) {
1851 if(!constTypes_ || constTypes_[i] == TNLP::NON_LINEAR) {
1852 double rowViolation = 0;
1853 if(rowLower[i] > -1e10)
1854 rowViolation = std::max(0.,rowLower[i] - g[i]);
1855
1856 if(rowUpper[i] < 1e10);
1857 rowViolation = std::max(rowViolation, g[i] - rowUpper[i]);
1858
1859 norm = rowViolation > norm ? rowViolation : norm;
1860 }
1861 }
1862 tminlp_->eval_f(numcols, x, 1, obj);
1863 delete [] g;
1864 return norm;
1865 }
1866
1867 /** get infinity norm of constraint violation of a point and objective error*/
1868 double
getNonLinearitiesViolation(const double * x,const double obj)1869 OsiTMINLPInterface::getNonLinearitiesViolation(const double *x, const double obj)
1870 {
1871 double f;
1872 double norm = getConstraintsViolation(x, f);
1873 assert((f - obj) > -1e-08);
1874 norm = (f - obj) > norm ? f - obj : norm;
1875 return norm;
1876 }
1877
1878
1879
1880 //A procedure to try to remove small coefficients in OA cuts (or make it non small
1881 static inline
cleanNnz(double & value,double colLower,double colUpper,double rowLower,double rowUpper,double colsol,double & lb,double & ub,double tiny,double veryTiny,double infty)1882 bool cleanNnz(double &value, double colLower, double colUpper,
1883 double rowLower, double rowUpper, double colsol,
1884 double & lb, double &ub, double tiny, double veryTiny,
1885 double infty)
1886 {
1887 //return 1;
1888 if(fabs(value)>= tiny) return 1;
1889
1890 if(fabs(value)<veryTiny) return 0;//Take the risk?
1891
1892 //try and remove
1893 bool colUpBounded = colUpper < infty;
1894 bool colLoBounded = colLower > -infty;
1895 bool rowNotLoBounded = rowLower <= - infty;
1896 bool rowNotUpBounded = rowUpper >= infty;
1897 bool pos = value > 0;
1898 if(colUpBounded && pos && rowNotUpBounded) {
1899 lb += value * (colsol - colUpper);
1900 return 0;
1901 }
1902 else
1903 if(colUpBounded && !pos && rowNotLoBounded) {
1904 ub += value * (colsol - colUpper);
1905 return 0;
1906 }
1907 else
1908 if(colLoBounded && !pos && rowNotUpBounded) {
1909 lb += value * (colsol - colLower);
1910 return 0;
1911 }
1912 else
1913 if(colLoBounded && pos && rowNotLoBounded) {
1914 ub += value * (colsol - colLower);
1915 return 0;
1916 }
1917 //can not remove coefficient
1918 return 1;
1919 }
1920
1921 /** Get the outer approximation constraints at the point x.
1922 */
1923 void
getOuterApproximation(OsiCuts & cs,const double * x,int getObj,const double * x2,double theta,bool global)1924 OsiTMINLPInterface::getOuterApproximation(OsiCuts &cs, const double * x,
1925 int getObj, const double * x2,
1926 double theta, bool global)
1927 {
1928 if(IsValid(linearizer_) && x2 == NULL){
1929 linearizer_->get_oas(cs, x, getObj, global);
1930 return;
1931 }
1932 int n,m, nnz_jac_g, nnz_h_lag;
1933 TNLP::IndexStyleEnum index_style;
1934 problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
1935 if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
1936 initializeJacobianArrays();
1937 assert(jRow_ != NULL);
1938 assert(jCol_ != NULL);
1939 vector<double> g(m);
1940 problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
1941 problem_to_optimize_->eval_g(n,x,1,m,g());
1942 //As jacobian is stored by cols fill OsiCuts with cuts
1943 vector<CoinPackedVector> cuts(nNonLinear_ + 1);
1944 vector<double> lb(nNonLinear_ + 1);
1945 vector<double> ub(nNonLinear_ + 1);
1946
1947 vector<int> row2cutIdx(m,-1);//store correspondance between index of row and index of cut (some cuts are not generated for rows because linear, or not binding). -1 if constraint does not generate a cut, otherwise index in cuts.
1948 int numCuts = 0;
1949
1950 const double * rowLower = getRowLower();
1951 const double * rowUpper = getRowUpper();
1952 const double * colLower = getColLower();
1953 const double * colUpper = getColUpper();
1954 const double * duals = getRowPrice() + 2 * n;
1955 double infty = getInfinity();
1956
1957 for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
1958 if(constTypes_[rowIdx] == TNLP::NON_LINEAR) {
1959 row2cutIdx[rowIdx] = numCuts;
1960 if(rowLower[rowIdx] > - infty_)
1961 lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
1962 else
1963 lb[numCuts] = - infty;
1964 if(rowUpper[rowIdx] < infty_)
1965 ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
1966 else
1967 ub[numCuts] = infty;
1968 if(rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty)
1969 {
1970 if(duals[rowIdx] >= 0)// <= inequality
1971 lb[numCuts] = - infty;
1972 if(duals[rowIdx] <= 0)// >= inequality
1973 ub[numCuts] = infty;
1974 }
1975
1976 numCuts++;
1977 }
1978 }
1979
1980
1981 for(int i = 0 ; i < nnz_jac_g ; i++) {
1982 const int &rowIdx = jRow_[i];
1983 const int & cutIdx = row2cutIdx[ rowIdx ];
1984 if(cutIdx != -1) {
1985 const int &colIdx = jCol_[i];
1986 //"clean" coefficient
1987 if(cleanNnz(jValues_[i],colLower[colIdx], colUpper[colIdx],
1988 rowLower[rowIdx], rowUpper[rowIdx],
1989 x[colIdx],
1990 lb[cutIdx],
1991 ub[cutIdx], tiny_, veryTiny_, infty_)) {
1992 cuts[cutIdx].insert(colIdx,jValues_[i]);
1993 if(lb[cutIdx] > - infty)
1994 lb[cutIdx] += jValues_[i] * x[colIdx];
1995 if(ub[cutIdx] < infty)
1996 ub[cutIdx] += jValues_[i] * x[colIdx];
1997 }
1998 }
1999 }
2000
2001 vector<int> cut2rowIdx(0);
2002 if (IsValid(cutStrengthener_) || oaHandler_->logLevel() > 0) {
2003 cut2rowIdx.resize(numCuts);// Store correspondance between indices of cut and indices of rows. For each cut
2004 for(int rowIdx = 0 ; rowIdx < m ; rowIdx++){
2005 if(row2cutIdx[rowIdx] >= 0){
2006 cut2rowIdx[row2cutIdx[rowIdx]] = rowIdx;
2007 }
2008 }
2009 }
2010
2011 for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
2012 //Compute cut violation
2013 if(x2 != NULL) {
2014 double rhs = cuts[cutIdx].dotProduct(x2);
2015 double violation = 0.;
2016 violation = std::max(violation, rhs - ub[cutIdx]);
2017 violation = std::max(violation, lb[cutIdx] - rhs);
2018 if(violation < theta && oaHandler_->logLevel() > 0) {
2019 oaHandler_->message(CUT_NOT_VIOLATED_ENOUGH, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
2020 continue;}
2021 if(oaHandler_->logLevel() > 0)
2022 oaHandler_->message(VIOLATED_OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
2023 }
2024 OsiRowCut newCut;
2025 // if(lb[i]>-1e20) assert (ub[i]>1e20);
2026
2027 if (IsValid(cutStrengthener_)) {
2028 const int& rowIdx = cut2rowIdx[cutIdx];
2029 bool retval =
2030 cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
2031 GetRawPtr(problem_), rowIdx,
2032 cuts[cutIdx], lb[cutIdx], ub[cutIdx], g[rowIdx],
2033 rowLower[rowIdx], rowUpper[rowIdx],
2034 n, x, infty);
2035 if (!retval) {
2036 (*messageHandler()) << "error in cutStrengthener_->ComputeCuts\n";
2037 //exit(-2);
2038 }
2039 }
2040 if(global) {
2041 newCut.setGloballyValidAsInteger(1);
2042 }
2043 if(lb[cutIdx] > infty) lb[cutIdx] -= rhsRelax_*std::max(fabs(lb[cutIdx]), 1.);
2044 if(ub[cutIdx] < infty) ub[cutIdx] += rhsRelax_*std::max(fabs(ub[cutIdx]), 1.);
2045 newCut.setLb(lb[cutIdx]);
2046 newCut.setUb(ub[cutIdx]);
2047 newCut.setRow(cuts[cutIdx]);
2048 if(oaHandler_->logLevel()>2){
2049 oaHandler_->print(newCut);}
2050 cs.insert(newCut);
2051 }
2052
2053 if(getObj == 2 || (getObj && !problem_->hasLinearObjective())) { // Get the objective cuts
2054 vector<double> obj(n);
2055 problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2056 double f;
2057 problem_to_optimize_->eval_f(n, x, 1, f);
2058
2059 CoinPackedVector v;
2060 v.reserve(n);
2061 lb[nNonLinear_] = -f;
2062 ub[nNonLinear_] = -f;
2063 //double minCoeff = 1e50;
2064 for(int i = 0; i<n ; i++) {
2065 if(cleanNnz(obj[i],colLower[i], colUpper[i],
2066 -getInfinity(), 0,
2067 x[i],
2068 lb[nNonLinear_],
2069 ub[nNonLinear_],tiny_, 1e-15, infty_)) {
2070 v.insert(i,obj[i]);
2071 lb[nNonLinear_] += obj[i] * x[i];
2072 ub[nNonLinear_] += obj[i] * x[i];
2073 }
2074 }
2075 v.insert(n,-1);
2076 //Compute cut violation
2077 bool genCut = true;
2078 if(x2 != NULL) {
2079 double rhs = v.dotProduct(x2);
2080 double violation = std::max(0., rhs - ub[nNonLinear_]);
2081 if(violation < theta) genCut = false;
2082 }
2083 if(genCut) {
2084 if (IsValid(cutStrengthener_)) {
2085 lb[nNonLinear_] = -infty;
2086 bool retval =
2087 cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
2088 GetRawPtr(problem_), -1,
2089 v, lb[nNonLinear_], ub[nNonLinear_],
2090 ub[nNonLinear_], -infty, 0.,
2091 n, x, infty);
2092 if (!retval) {
2093 (*handler_)<< "error in cutStrengthener_->ComputeCuts"<<CoinMessageEol;
2094 //exit(-2);
2095 }
2096 }
2097 OsiRowCut newCut;
2098 if(global)
2099 newCut.setGloballyValidAsInteger(1);
2100 //newCut.setEffectiveness(99.99e99);
2101 newCut.setRow(v);
2102 newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
2103 newCut.setUb(ub[nNonLinear_]);
2104 cs.insert(newCut);
2105 }
2106 }
2107 }
2108
2109 /** Get a benders cut from solution.*/
2110 void
getBendersCut(OsiCuts & cs,bool global)2111 OsiTMINLPInterface::getBendersCut(OsiCuts &cs,
2112 bool global){
2113 int n,m, nnz_jac_g, nnz_h_lag;
2114 TNLP::IndexStyleEnum index_style;
2115 problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2116 if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2117 initializeJacobianArrays();
2118 assert(jRow_ != NULL);
2119 assert(jCol_ != NULL);
2120 vector<double> g(m);
2121 const double * x = getColSolution();
2122 problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
2123 problem_to_optimize_->eval_g(n,x,1,m,g());
2124 //As jacobian is stored by cols fill OsiCuts with cuts
2125 vector<double> cut(n+1,0.);
2126 vector<bool> keep(m+1,false);
2127 double lb = 0;
2128 double ub = 0;
2129
2130 const double * rowLower = getRowLower();
2131 const double * rowUpper = getRowUpper();
2132 const double * colLower = getColLower();
2133 const double * colUpper = getColUpper();
2134 const double * duals = getRowPrice() + 2 * n;
2135 //double infty = getInfinity();
2136
2137 for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
2138 if(constTypes_[rowIdx] == TNLP::NON_LINEAR && fabs(duals[rowIdx]) > 1e-06)
2139 {
2140 const double & lam = duals[rowIdx];
2141 keep[rowIdx] = true;
2142 assert(lam < 0 || rowUpper[rowIdx] < 1e10);
2143 assert(lam > 0 || rowLower[rowIdx] > -1e10);
2144 if(lam < 0){
2145 assert(rowLower[rowIdx] > -1e10);
2146 ub += lam*(rowLower[rowIdx] -g[rowIdx]);
2147 }
2148 else {
2149 assert(rowUpper[rowIdx] < 1e10);
2150 ub += lam*(rowUpper[rowIdx] -g[rowIdx]);
2151 }
2152 }
2153 }
2154
2155
2156 for(int i = 0 ; i < nnz_jac_g ; i++) {
2157 const int &rowIdx = jRow_[i];
2158 if (!keep[rowIdx]) continue;
2159 const int &colIdx = jCol_[i];
2160 //"clean" coefficient
2161 const double & lam = duals[rowIdx];
2162 double coeff = lam*jValues_[i];
2163 if(cleanNnz(coeff,colLower[colIdx], colUpper[colIdx],
2164 rowLower[rowIdx], rowUpper[rowIdx], x[colIdx], lb,
2165 ub, tiny_, veryTiny_, infty_)) {
2166 cut[colIdx] += coeff;
2167 ub += coeff * x[colIdx];
2168 }
2169 }
2170
2171 CoinPackedVector v;
2172 if(!problem_->hasLinearObjective() &&
2173 isProvenOptimal()) { // Get the objective cuts
2174 vector<double> obj(n);
2175 problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2176 double f;
2177 problem_to_optimize_->eval_f(n, x, 1, f);
2178
2179 ub = -f;
2180 //double minCoeff = 1e50;
2181 for(int i = 0; i<n ; i++) {
2182 if(cleanNnz(obj[i],colLower[i], colUpper[i], -getInfinity(), 0,
2183 x[i], lb, ub,tiny_, 1e-15, infty_)) {
2184 cut[i] += obj[i];
2185 ub += obj[i] * x[i];
2186 }
2187 }
2188 v.insert(n,-1);
2189 }
2190 for(int i = 0 ; i < n ; i++){
2191 if(fabs(cut[i])>1e-020){
2192 v.insert(i, cut[i]);
2193 }
2194 }
2195 OsiRowCut newCut;
2196 if(global)
2197 newCut.setGloballyValidAsInteger(1);
2198 newCut.setLb(-COIN_DBL_MAX/*Infinity*/);
2199 newCut.setUb(ub);
2200 newCut.setRow(v);
2201 cs.insert(newCut);
2202 }
2203
2204 /** Get the outer approximation of a single constraint at the point x.
2205 */
2206 void
getConstraintOuterApproximation(OsiCuts & cs,int rowIdx,const double * x,const double * x2,bool global)2207 OsiTMINLPInterface::getConstraintOuterApproximation(OsiCuts &cs, int rowIdx,
2208 const double * x,
2209 const double * x2, bool global)
2210 {
2211 double g;
2212 int * indices = new int[getNumCols()];
2213 double * values = new double[getNumCols()];
2214 int nnz;
2215 problem_->eval_grad_gi(getNumCols(), x, 1, rowIdx, nnz, indices, values);
2216 problem_->eval_gi(getNumCols(),x,1, rowIdx, g);
2217
2218 CoinPackedVector cut;
2219 double lb;
2220 double ub;
2221
2222
2223 const double rowLower = getRowLower()[rowIdx];
2224 const double rowUpper = getRowUpper()[rowIdx];
2225 const double * colLower = getColLower();
2226 const double * colUpper = getColUpper();
2227 const double dual = (getRowPrice() + 2 * getNumCols())[rowIdx];
2228 double infty = getInfinity();
2229
2230 if(rowLower > - infty_)
2231 lb = rowLower - g;
2232 else
2233 lb = - infty;
2234 if(rowUpper < infty_)
2235 ub = rowUpper - g;
2236 else
2237 ub = infty;
2238 if(rowLower > -infty && rowUpper < infty)
2239 {
2240 if(dual >= 0)// <= inequality
2241 lb = - infty;
2242 if(dual <= 0)// >= inequality
2243 ub = infty;
2244 }
2245
2246 for(int i = 0 ; i < nnz; i++) {
2247 const int &colIdx = indices[i];
2248 //"clean" coefficient
2249 if(cleanNnz(values[i],colLower[colIdx], colUpper[colIdx],
2250 rowLower, rowUpper,
2251 x[colIdx],
2252 lb,
2253 ub, tiny_, veryTiny_, infty_)) {
2254 cut.insert(colIdx,values[i]);
2255 if(lb > - infty)
2256 lb += values[i] * x[colIdx];
2257 if(ub < infty)
2258 ub += values[i] * x[colIdx];
2259 }
2260 }
2261
2262 OsiRowCut newCut;
2263
2264 if(global) {
2265 newCut.setGloballyValidAsInteger(1);
2266 }
2267 //newCut.setEffectiveness(99.99e99);
2268 newCut.setLb(lb);
2269 newCut.setUb(ub);
2270 newCut.setRow(cut);
2271 cs.insert(newCut);
2272
2273 delete [] indices;
2274 delete [] values;
2275 }
2276
2277 void
switchToFeasibilityProblem(size_t n,const double * x_bar,const int * inds,double a,double s,int L)2278 OsiTMINLPInterface::switchToFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
2279 double a, double s, int L){
2280 if(! IsValid(feasibilityProblem_)) {
2281 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2282 }
2283 feasibilityProblem_->set_use_feasibility_pump_objective(true);
2284 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
2285 feasibilityProblem_->setLambda(a);
2286 feasibilityProblem_->setSigma(s);
2287 feasibilityProblem_->setNorm(L);
2288 feasibilityProblem_->set_use_cutoff_constraint(false);
2289 feasibilityProblem_->set_use_local_branching_constraint(false);
2290 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
2291 feasibility_mode_ = true;
2292 }
2293
2294 void
switchToFeasibilityProblem(size_t n,const double * x_bar,const int * inds,double rhs_local_branching_constraint)2295 OsiTMINLPInterface::switchToFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
2296 double rhs_local_branching_constraint){
2297 if(! IsValid(feasibilityProblem_)) {
2298 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2299 }
2300 feasibilityProblem_->set_use_feasibility_pump_objective(false);
2301 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
2302 feasibilityProblem_->set_use_cutoff_constraint(false);
2303 feasibilityProblem_->set_use_local_branching_constraint(true);
2304 feasibilityProblem_->set_rhs_local_branching_constraint(rhs_local_branching_constraint);
2305 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
2306 feasibility_mode_ = true;
2307 }
2308
2309 void
switchToOriginalProblem()2310 OsiTMINLPInterface::switchToOriginalProblem(){
2311 problem_to_optimize_ = GetRawPtr(problem_);
2312 feasibility_mode_ = false;
2313 }
2314
2315 double
solveFeasibilityProblem(size_t n,const double * x_bar,const int * inds,double a,double s,int L)2316 OsiTMINLPInterface::solveFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
2317 double a, double s, int L)
2318 {
2319 if(! IsValid(feasibilityProblem_)) {
2320 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2321 }
2322 feasibilityProblem_->set_use_feasibility_pump_objective(true);
2323 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
2324 feasibilityProblem_->setLambda(a);
2325 feasibilityProblem_->setSigma(s);
2326 feasibilityProblem_->setNorm(L);
2327 feasibilityProblem_->set_use_cutoff_constraint(false);
2328 feasibilityProblem_->set_use_local_branching_constraint(false);
2329 nCallOptimizeTNLP_++;
2330 totalNlpSolveTime_-=CoinCpuTime();
2331 SmartPtr<TNLPSolver> app2 = app_->clone();
2332 app2->options()->SetIntegerValue("print_level", (Index) 0);
2333 optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2334 totalNlpSolveTime_+=CoinCpuTime();
2335 hasBeenOptimized_=true;
2336 return getObjValue();
2337 }
2338
2339 double
solveFeasibilityProblem(size_t n,const double * x_bar,const int * inds,int L,double cutoff)2340 OsiTMINLPInterface::solveFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
2341 int L, double cutoff)
2342 {
2343 if(! IsValid(feasibilityProblem_)) {
2344 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
2345 }
2346 feasibilityProblem_->set_use_feasibility_pump_objective(true);
2347 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
2348 feasibilityProblem_->setLambda(1.0);
2349 feasibilityProblem_->setSigma(0.0);
2350 feasibilityProblem_->setNorm(L);
2351 feasibilityProblem_->set_use_cutoff_constraint(true);
2352 feasibilityProblem_->set_cutoff(cutoff);
2353 feasibilityProblem_->set_use_local_branching_constraint(false);
2354 nCallOptimizeTNLP_++;
2355 totalNlpSolveTime_-=CoinCpuTime();
2356 SmartPtr<TNLPSolver> app2 = app_->clone();
2357 app2->options()->SetIntegerValue("print_level", (Index) 0);
2358 optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
2359 totalNlpSolveTime_+=CoinCpuTime();
2360 hasBeenOptimized_=true;
2361 return getObjValue();
2362 }
2363
2364 #if 0
2365 double
2366 OsiTMINLPInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs, bool addOnlyViolated, bool global)
2367 {
2368 double ret_val = solveFeasibilityProblem(n, x_bar, inds, 1, 0, 2);
2369 getOuterApproximation(cs, getColSolution(), 0, (addOnlyViolated? x_bar:NULL)
2370 , global);
2371 return ret_val;
2372 }
2373 #endif
2374
2375 static bool WarnedForNonConvexOa=false;
2376
2377 void
extractLinearRelaxation(OsiSolverInterface & si,const double * x,bool getObj)2378 OsiTMINLPInterface::extractLinearRelaxation(OsiSolverInterface &si,
2379 const double * x, bool getObj
2380 )
2381 {
2382 if(IsValid(linearizer_)){
2383 linearizer_->extract(&si, x, getObj);
2384 return;
2385 }
2386 int n;
2387 int m;
2388 int nnz_jac_g;
2389 int nnz_h_lag;
2390 TNLP::IndexStyleEnum index_style;
2391 //Get problem information
2392 problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
2393
2394 //if not allocated allocate spaced for stroring jacobian
2395 //if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
2396 initializeJacobianArrays();
2397
2398 //get Jacobian
2399 problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
2400
2401
2402 vector<double> g(m);
2403 problem_to_optimize_->eval_g(n, x, 1, m, g());
2404
2405 vector<double> rowLow(m);
2406 vector<double> rowUp(m);
2407 vector<int> nonBindings(m);;//store non binding constraints (which are to be removed from OA)
2408 int numNonBindings = 0;
2409 const double * rowLower = getRowLower();
2410 const double * rowUpper = getRowUpper();
2411 vector<double> colLower(n);
2412 vector<double> colUpper(n);
2413 std::copy(getColLower(), getColLower() + n, colLower.begin());
2414 std::copy(getColUpper(), getColUpper() + n, colUpper.begin());
2415 const double * duals = getRowPrice() + 2*n;
2416 assert(m==getNumRows());
2417 double infty = si.getInfinity();
2418 for(int i = 0 ; i < m ; i++) {
2419 if(constTypes_[i] == TNLP::NON_LINEAR) {
2420 //If constraint is range not binding prepare to remove it
2421 if(rowLower[i] > -infty_ && rowUpper[i] < infty_ && fabs(duals[i]) == 0.)
2422 {
2423 nonBindings[numNonBindings++] = i;
2424 continue;
2425 }
2426 else
2427 if(rowLower[i] > - infty_){
2428 rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
2429 if(! WarnedForNonConvexOa && rowUpper[i] < infty_){
2430 messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
2431 WarnedForNonConvexOa = true;
2432 }
2433 }
2434 else
2435 rowLow[i] = - infty;
2436 if(rowUpper[i] < infty_)
2437 rowUp[i] = (rowUpper[i] - g[i]) + 1e-07;
2438 else
2439 rowUp[i] = infty;
2440
2441 //If equality or ranged constraint only add one side by looking at sign of dual multiplier
2442 if(rowLower[i] > -infty_ && rowUpper[i] < infty_)
2443 {
2444 if(duals[i] >= 0.)// <= inequality
2445 rowLow[i] = - infty;
2446 if(duals[i] <= 0.)// >= inequality
2447 rowUp[i] = infty;
2448 }
2449 }
2450 else {
2451 if(rowLower[i] > -infty_){
2452 rowLow[i] = (rowLower[i]);
2453 }
2454 else
2455 rowLow[i] = - infty;
2456 if(rowUpper[i] < infty_){
2457 rowUp[i] = (rowUpper[i]);
2458 }
2459 else
2460 rowUp[i] = infty;
2461 }
2462 }
2463
2464
2465
2466 //Then convert everything to a CoinPackedMatrix
2467 //Go through values, clean coefficients and fix bounds
2468 for(int i = 0 ; i < nnz_jac_g ; i++) {
2469 if(constTypes_[jRow_[i]] != TNLP::LINEAR){//For linear just copy is fine.
2470 if(//For other clean tinys
2471 cleanNnz(jValues_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
2472 rowLower[jRow_[i]], rowUpper[jRow_[i]],
2473 x[jCol_[i]],
2474 rowLow[jRow_[i]],
2475 rowUp[jRow_[i]], tiny_, veryTiny_, infty_)) {
2476 if(rowLow[jRow_[i]] > - infty)
2477 rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
2478 if(rowUp[jRow_[i]] < infty)
2479 rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
2480 }
2481 }
2482 else {
2483 //double value = jValues_[i] * x[jCol_[i]];
2484 //rowLow[jRow_[i]] += value;
2485 //rowUp[jRow_[i]] += value;
2486 }
2487 }
2488
2489
2490 CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
2491 mat.setDimensions(m,n); // In case matrix was empty, this should be enough
2492
2493 //remove non-bindings equality constraints
2494 mat.deleteRows(numNonBindings, nonBindings());
2495
2496 int numcols=getNumCols();
2497 vector<double> obj(numcols);
2498 for(int i = 0 ; i < numcols ; i++){
2499 obj[i] = 0.;
2500 if(colLower[i] <= - infty_) colLower[i] = - infty;
2501 if(colUpper[i] >= infty_) colUpper[i] = infty;
2502 }
2503
2504
2505 //Relax rhs's by tiny_
2506 for(size_t i = 0 ; i < rowLow.size() ; i++){
2507 if(rowLow[i] > infty) rowLow[i] -= rhsRelax_*std::max(fabs(rowLow[i]), 1.);
2508 if(rowUp[i] < infty) rowUp[i] += rhsRelax_*std::max(fabs(rowUp[i]), 1.);
2509 }
2510
2511 si.loadProblem(mat, colLower(), colUpper(), obj(), rowLow(), rowUp());
2512 for(int i = 0 ; i < numcols ; i++) {
2513 if(isInteger(i))
2514 si.setInteger(i);
2515 }
2516 if(getObj) {
2517 bool addObjVar = false;
2518 if(problem_->hasLinearObjective()){
2519 double zero;
2520 vector<double> x0(n,0.);
2521 problem_to_optimize_->eval_f(n, x0(), 1, zero);
2522 si.setDblParam(OsiObjOffset, -zero);
2523 //Copy the linear objective and don't create a dummy variable.
2524 problem_to_optimize_->eval_grad_f(n, x, 1,obj());
2525 si.setObjective(obj());
2526 }
2527 else {
2528 addObjVar = true;
2529 }
2530
2531 if(addObjVar){
2532 addObjectiveFunction(si, x);
2533 }
2534 }
2535
2536 //si.writeMpsNative("Toto.mps", NULL, NULL, 1);
2537 }
2538
2539 void
addObjectiveFunction(OsiSolverInterface & si,const double * x)2540 OsiTMINLPInterface::addObjectiveFunction(OsiSolverInterface &si,
2541 const double *x){
2542 const double * colLower = getColLower();
2543 const double * colUpper = getColUpper();
2544 int numcols = getNumCols();
2545 assert(numcols == si.getNumCols() );
2546 vector<double> obj(numcols);
2547 problem_to_optimize_->eval_grad_f(numcols, x, 1,obj());
2548 //add variable alpha
2549 //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
2550 CoinPackedVector a;
2551 si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
2552
2553 // Now get the objective cuts
2554 // get the gradient, pack it and add the cut
2555 double ub;
2556 problem_to_optimize_->eval_f(numcols, x, 1, ub);
2557 ub*=-1;
2558 double lb = -1e300;
2559 CoinPackedVector objCut;
2560 CoinPackedVector * v = &objCut;
2561 v->reserve(numcols+1);
2562 for(int i = 0; i<numcols ; i++) {
2563 if(si.getNumRows())
2564 {
2565 if(cleanNnz(obj[i],colLower[i], colUpper[i],
2566 -getInfinity(), 0,
2567 x[i],
2568 lb,
2569 ub, tiny_, veryTiny_, infty_)) {
2570 v->insert(i,obj[i]);
2571 lb += obj[i] * x[i];
2572 ub += obj[i] * x[i];
2573 }
2574 }
2575 else //Unconstrained problem can not put clean coefficient
2576 {
2577 if(cleanNnz(obj[i],colLower[i], colUpper[i],
2578 -getInfinity(), 0,
2579 x[i],
2580 lb,
2581 ub, 1e-03, 1e-08, infty_)) {
2582 v->insert(i,obj[i]);
2583 lb += obj[i] * x[i];
2584 ub += obj[i] * x[i];
2585 }
2586 }
2587 }
2588 v->insert(numcols,-1);
2589 si.addRow(objCut, lb, ub);
2590 }
2591
2592 /** Add a collection of linear cuts to problem formulation.*/
2593 void
applyRowCuts(int numberCuts,const OsiRowCut * cuts)2594 OsiTMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
2595 {
2596 if(numberCuts)
2597 freeCachedRowRim();
2598 const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
2599 for(int i = 0 ; i < numberCuts ; i++)
2600 {
2601 cutsPtrs[i] = &cuts[i];
2602 }
2603 problem_->addCuts(numberCuts, cutsPtrs);
2604 delete [] cutsPtrs;
2605 }
2606
2607 void
solveAndCheckErrors(bool warmStarted,bool throwOnFailure,const char * whereFrom)2608 OsiTMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
2609 const char * whereFrom)
2610 {
2611 if (BonminAbortAll == true) return;
2612 totalNlpSolveTime_-=CoinCpuTime();
2613 if(warmStarted)
2614 optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
2615 else
2616 optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
2617 totalNlpSolveTime_+=CoinCpuTime();
2618 nCallOptimizeTNLP_++;
2619 hasBeenOptimized_ = true;
2620
2621 if(getRowCutDebugger()){
2622 //printf("On the optimal path %g < %g?\n", getObjValue(), getRowCutDebugger()->optimalValue());
2623 if(! (isProvenOptimal() && getObjValue() <= getRowCutDebugger()->optimalValue())){
2624 std::string probName;
2625 getStrParam(OsiProbName, probName);
2626 fprintf(stderr, "Problem on optimal path is infeasible!\n");
2627 throw newUnsolvedError(app_->errorCode(), problem_, probName);
2628 }
2629 }
2630
2631 //Options should have been printed if not done already turn off Ipopt output
2632 if(!hasPrintedOptions) {
2633 hasPrintedOptions = 1;
2634 //app_->Options()->SetIntegerValue("print_level",0, true, true);
2635 app_->options()->SetStringValue("print_user_options","no", false, true);
2636 }
2637
2638 bool otherDisagree = false ;
2639 if(!app_->isRecoverable(optimizationStatus_))//Solver failed and the error can not be recovered, throw it
2640 {
2641 std::string probName;
2642 getStrParam(OsiProbName, probName);
2643 throw newUnsolvedError(app_->errorCode(), problem_, probName);
2644 }
2645 else if(testOthers_ && !app_->isError(optimizationStatus_)){
2646 Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
2647 //Try other solvers and see if they agree
2648 int f =1;
2649 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2650 i != debug_apps_.end() ; i++){
2651 TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
2652 messageHandler()->message(LOG_LINE, messages_)
2653 <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
2654 <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
2655 if(!(*i)->isError(otherStatus)){
2656 CoinRelFltEq eq(1e-05);
2657 if(otherStatus != optimizationStatus_){
2658 otherDisagree = true;
2659 messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
2660 <<app_->solverName()<<statusAsString()
2661 <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol;
2662 }
2663 else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
2664 {
2665 otherDisagree = true;
2666 messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
2667 <<app_->solverName()<<problem_->obj_value()
2668 <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol;
2669 }
2670 }
2671 }
2672 }
2673 else if(app_->isError(optimizationStatus_) && ! debug_apps_.empty()){
2674 int f =1;
2675 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
2676 i != debug_apps_.end() && app_->isError(optimizationStatus_) ; i++){
2677 optimizationStatus_ = (*i)->OptimizeTNLP(GetRawPtr(problem_));
2678 messageHandler()->message(LOG_LINE, messages_)
2679 <<'d'<<f++<<statusAsString(optimizationStatus_)<<problem_->obj_value()
2680 <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
2681 }
2682 }
2683 try{
2684 totalIterations_ += app_->IterationCount();
2685 }
2686 catch(SimpleError &E)
2687 {
2688 if (throwOnFailure)//something failed throw
2689 {
2690 throw SimpleError("No statistics available from Ipopt",whereFrom);
2691 }
2692 else {
2693 return;
2694 }
2695 }
2696 if(problem_->hasUpperBoundingObjective()){//Check if solution is integer and recompute objective value using alternative objective function
2697 const double * sol = getColSolution();
2698 bool integerSol = true;
2699 double intTol = 1e-08;
2700 if(objects()){
2701 int nObjects = numberObjects();
2702 OsiObject ** object = objects();
2703 for(int i = 0 ; i< nObjects ; i++){
2704 int dummy;
2705 if(object[i]->infeasibility(this,dummy) > intTol)
2706 {
2707 integerSol=false;
2708 break;
2709 }
2710 }
2711 }
2712 else{//Only works for integer constraints
2713 int numcols = getNumCols();
2714 for(int i = 0 ; i < numcols ; i++){
2715 if(isInteger(i) || isBinary(i)){
2716 if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
2717 integerSol = false;
2718 break;
2719 }
2720 }
2721 }
2722 }
2723 if(integerSol&&isProvenOptimal()){
2724 double help= problem_->evaluateUpperBoundingFunction(sol);
2725
2726
2727 OsiAuxInfo * auxInfo = getAuxiliaryInfo();
2728 Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
2729 if(bonInfo!=0)
2730 {
2731
2732 if(help<bonInfo->bestObj2())
2733 {
2734 bonInfo->setBestObj2(help);
2735 bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
2736
2737 messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
2738 <<help<<CoinMessageEol;
2739 }
2740 }
2741 else {
2742 printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
2743 }
2744 }
2745 }
2746
2747 messageHandler()->message(IPOPT_SUMMARY, messages_)
2748 <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
2749
2750 if((nCallOptimizeTNLP_ % 20) == 1)
2751 messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
2752
2753
2754 if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
2755 ( otherDisagree )){
2756 messageHandler()->message(SUSPECT_PROBLEM,
2757 messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
2758 std::string subProbName;
2759 getStrParam(OsiProbName, subProbName);
2760 std::ostringstream os;
2761 os<<"_"<<nCallOptimizeTNLP_;
2762 subProbName+=os.str();
2763 problem_->outputDiffs(subProbName, NULL/*getVarNames()*/);
2764 }
2765
2766 }
2767
2768 ////////////////////////////////////////////////////////////////////
2769 // Solve Methods //
2770 ////////////////////////////////////////////////////////////////////
initialSolve()2771 void OsiTMINLPInterface::initialSolve()
2772 {
2773 initialSolve("");
2774 }
2775
2776 ////////////////////////////////////////////////////////////////////
2777 // Solve Methods //
2778 ////////////////////////////////////////////////////////////////////
initialSolve(const char * whereFrom)2779 void OsiTMINLPInterface::initialSolve(const char * whereFrom)
2780 {
2781 assert(IsValid(app_));
2782 assert(IsValid(problem_));
2783
2784 if (BonminAbortAll == true) return;
2785 // Discard warmstart_ if we had one
2786 delete warmstart_;
2787 warmstart_ = NULL;
2788
2789 if(!hasPrintedOptions) {
2790 int printOptions;
2791 app_->options()->GetEnumValue("print_user_options",printOptions,app_->prefix());
2792 if(printOptions)
2793 app_->options()->SetStringValue("print_user_options","yes",true,true);
2794 }
2795 if(warmStartMode_ >= Optimum)
2796 app_->disableWarmStart();
2797 solveAndCheckErrors(0,1,"initialSolve");
2798
2799 //Options should have been printed if not done already turn off Ipopt output
2800 if(!hasPrintedOptions) {
2801 hasPrintedOptions = 1;
2802 app_->options()->SetStringValue("print_user_options","no");
2803 app_->options()->SetIntegerValue("print_level",0);
2804 }
2805
2806 messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
2807 <<statusAsString()
2808 <<getObjValue()
2809 <<app_->IterationCount()
2810 <<app_->CPUTime()
2811 <<whereFrom<<CoinMessageEol;
2812
2813 if(BonminAbortAll){
2814 return;
2815 }
2816 int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
2817 if(isAbandoned() ||
2818 ( isProvenPrimalInfeasible() && getObjValue() < infeasibility_epsilon_)) {
2819 resolveForRobustness(numRetryUnsolved_);
2820 }
2821 else if(numRetry)
2822 {
2823 resolveForCost(numRetry, numRetryInitial_ > 0);
2824 /** Only do it once at the root.*/
2825 numRetryInitial_ = 0;
2826 }
2827 firstSolve_ = false;
2828
2829 // if warmstart_ is not empty then had to use resolveFor... and that created
2830 // the warmstart at the end, and we have nothing to do here. Otherwise...
2831 if (! warmstart_ && ! isAbandoned()) {
2832 if (warmStartMode_ >= Optimum) {
2833 warmstart_ = app_->getWarmStart(problem_);
2834 }
2835 }
2836 }
2837
2838 /** Resolve the continuous relaxation after problem modification. */
2839 void
resolve()2840 OsiTMINLPInterface::resolve(){
2841 resolve("");
2842 }
2843 /** Resolve the continuous relaxation after problem modification.
2844 * \note for Ipopt, same as resolve */
2845 void
resolve(const char * whereFrom)2846 OsiTMINLPInterface::resolve(const char * whereFrom)
2847 {
2848 assert(IsValid(app_));
2849 assert(IsValid(problem_));
2850 if (BonminAbortAll == true) return;
2851 int has_warmstart = warmstart_ == NULL ? 0 : 1;
2852 if(warmstart_ == NULL) has_warmstart = 0;
2853 else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
2854 else has_warmstart = 2;
2855 if (has_warmstart < 2) {
2856 // empty (or unrecognized) warmstart
2857 initialSolve(whereFrom);
2858 return;
2859 }
2860 app_->setWarmStart(warmstart_, problem_);
2861 delete warmstart_;
2862 warmstart_ = NULL;
2863
2864 if (INT_BIAS > 0.) {
2865 app_->options()->SetStringValue("warm_start_same_structure", "yes");
2866 }
2867 else {
2868 app_->options()->SetStringValue("warm_start_same_structure", "no");
2869 }
2870
2871 if(problem_->duals_init() != NULL)
2872 app_->enableWarmStart();
2873 else app_->disableWarmStart();
2874 solveAndCheckErrors(1,1,"resolve");
2875
2876 messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
2877 <<statusAsString()
2878 <<getObjValue()
2879 <<app_->IterationCount()
2880 <<app_->CPUTime()
2881 <<whereFrom
2882 <<"totot"
2883 <<CoinMessageEol;
2884
2885 if(isAbandoned() ||
2886 ( (getObjValue() < 1e-06) && isProvenPrimalInfeasible())) {
2887 resolveForRobustness(numRetryUnsolved_);
2888 }
2889 else if(numRetryResolve_ ||
2890 (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
2891 resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
2892
2893 // if warmstart_ is not empty then had to use resolveFor... and that created
2894 // the warmstart at the end, and we have nothing to do here. Otherwise...
2895 if (! warmstart_ && ! isAbandoned()) {
2896 if (warmStartMode_ >= Optimum) {
2897 warmstart_ = app_->getWarmStart(problem_);
2898 }
2899 }
2900 }
2901
2902
2903 ////////////////////////////////////////////////////////////////
2904 // Methods returning info on how the solution process terminated //
2905 ////////////////////////////////////////////////////////////////
2906 /// Are there a numerical difficulties?
isAbandoned() const2907 bool OsiTMINLPInterface::isAbandoned() const
2908 {
2909 if (pretendSucceededNext_)
2910 return false;
2911 pretendSucceededNext_ = false;
2912 return (
2913 (optimizationStatus_==TNLPSolver::iterationLimit)||
2914 (optimizationStatus_==TNLPSolver::computationError)||
2915 (optimizationStatus_==TNLPSolver::illDefinedProblem)||
2916 (optimizationStatus_==TNLPSolver::illegalOption)||
2917 (optimizationStatus_==TNLPSolver::externalException)|
2918 (optimizationStatus_==TNLPSolver::exception)
2919 );
2920 }
2921
2922 /// Is optimality proven?
isProvenOptimal() const2923 bool OsiTMINLPInterface::isProvenOptimal() const
2924 {
2925 return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
2926 (optimizationStatus_==TNLPSolver::solvedOptimalTol);
2927 }
2928 /// Is primal infeasiblity proven?
isProvenPrimalInfeasible() const2929 bool OsiTMINLPInterface::isProvenPrimalInfeasible() const
2930 {
2931 return (optimizationStatus_ == TNLPSolver::provenInfeasible);
2932 }
2933 /// Is dual infeasiblity proven?
isProvenDualInfeasible() const2934 bool OsiTMINLPInterface::isProvenDualInfeasible() const
2935 {
2936 return (optimizationStatus_ == TNLPSolver::unbounded);
2937 }
2938 /// Is the given primal objective limit reached?
isPrimalObjectiveLimitReached() const2939 bool OsiTMINLPInterface::isPrimalObjectiveLimitReached() const
2940 {
2941 (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2942 return 0;
2943 }
2944 /// Is the given dual objective limit reached?
isDualObjectiveLimitReached() const2945 bool OsiTMINLPInterface::isDualObjectiveLimitReached() const
2946 {
2947 // (*messageHandler_)<<"Warning : isDualObjectiveLimitReached not implemented yet"<<CoinMessageEol;
2948 return (optimizationStatus_==TNLPSolver::unbounded);
2949
2950 }
2951 /// Iteration limit reached?
isIterationLimitReached() const2952 bool OsiTMINLPInterface::isIterationLimitReached() const
2953 {
2954 return (optimizationStatus_==TNLPSolver::iterationLimit);
2955 }
2956
2957 void
extractInterfaceParams()2958 OsiTMINLPInterface::extractInterfaceParams()
2959 {
2960 if (IsValid(app_)) {
2961 int logLevel;
2962 app_->options()->GetIntegerValue("nlp_log_level", logLevel,app_->prefix());
2963 messageHandler()->setLogLevel(logLevel);
2964
2965 #ifdef COIN_HAS_FILTERSQP
2966 FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
2967
2968 bool is_given =
2969 #endif
2970 app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,app_->prefix());
2971
2972 #ifdef COIN_HAS_FILTERSQP
2973 if(filter && !is_given){
2974 // overwriting default for filter
2975 maxRandomRadius_ = 10.;
2976 }
2977 #endif
2978
2979 int oaCgLogLevel = 0;
2980 app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,app_->prefix());
2981 oaHandler_->setLogLevel(oaCgLogLevel);
2982
2983 int buffy;
2984 app_->options()->GetEnumValue("warm_start", buffy, app_->prefix());
2985 warmStartMode_ = (WarmStartModes) buffy;
2986
2987 app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,app_->prefix());
2988 app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,app_->prefix());
2989 app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,app_->prefix());
2990 app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,app_->prefix());
2991 app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,app_->prefix());
2992 app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,app_->prefix());
2993 app_->options()->GetNumericValue
2994 ("warm_start_bound_frac" ,pushValue_,app_->prefix());
2995 app_->options()->GetNumericValue("tiny_element",tiny_,app_->prefix());
2996 app_->options()->GetNumericValue("very_tiny_element",veryTiny_,app_->prefix());
2997 app_->options()->GetNumericValue("oa_rhs_relax",rhsRelax_,app_->prefix());
2998 app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,app_->prefix());
2999 app_->options()->GetEnumValue("random_point_type",randomGenerationType_,app_->prefix());
3000 int cut_strengthening_type;
3001 app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,app_->prefix());
3002 double lo_inf, up_inf;
3003 app_->options()->GetNumericValue("nlp_lower_bound_inf",lo_inf, app_->prefix());
3004 app_->options()->GetNumericValue("nlp_upper_bound_inf",up_inf, app_->prefix());
3005 infty_ = std::min(fabs(lo_inf), fabs(up_inf));
3006
3007 #ifdef COIN_HAS_FILTERSQP
3008 is_given =
3009 #endif
3010 app_->options()->GetNumericValue("resolve_on_small_infeasibility", infeasibility_epsilon_, app_->prefix());
3011
3012 #ifdef COIN_HAS_FILTERSQP
3013 if(filter && !is_given){
3014 // overwriting default for filter
3015 infeasibility_epsilon_ = 1e-06;
3016 std::string o_name = app_->prefix();
3017 o_name += "resolve_on_small_infeasibility";
3018 app_->options()->SetNumericValue(o_name.c_str(), infeasibility_epsilon_, true, true);
3019 }
3020 #endif
3021 if (cut_strengthening_type != CS_None) {
3022 // TNLP solver to be used in the cut strengthener
3023 cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
3024 }
3025 }
3026 }
3027
3028 void
SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)3029 OsiTMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
3030 {
3031 strong_branching_solver_ = strong_branching_solver;
3032 }
3033
3034 //#define STRONG_COMPARE
3035 #ifdef STRONG_COMPARE
3036 static double objorig;
3037 #endif
3038
3039 void
markHotStart()3040 OsiTMINLPInterface::markHotStart()
3041 {
3042 if (IsValid(strong_branching_solver_)) {
3043 #ifdef STRONG_COMPARE
3044 // AWDEBUG
3045 OsiSolverInterface::markHotStart();
3046 objorig = getObjValue();
3047 #endif
3048 optimizationStatusBeforeHotStart_ = optimizationStatus_;
3049 strong_branching_solver_->markHotStart(this);
3050 }
3051 else {
3052 // Default Implementation
3053 OsiSolverInterface::markHotStart();
3054 }
3055 }
3056
3057 void
solveFromHotStart()3058 OsiTMINLPInterface::solveFromHotStart()
3059 {
3060 if (IsValid(strong_branching_solver_)) {
3061 #ifdef STRONG_COMPARE
3062 // AWDEBUG
3063 OsiSolverInterface::solveFromHotStart();
3064 double obj_nlp = getObjValue() - objorig;
3065 #endif
3066 optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
3067 hasBeenOptimized_ = true;
3068 #ifdef STRONG_COMPARE
3069 double obj_other = getObjValue() - objorig;
3070 printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
3071 obj_nlp, obj_other);
3072 #endif
3073 }
3074 else {
3075 // Default Implementation
3076 OsiSolverInterface::solveFromHotStart();
3077 }
3078 }
3079
3080 void
unmarkHotStart()3081 OsiTMINLPInterface::unmarkHotStart()
3082 {
3083 if (IsValid(strong_branching_solver_)) {
3084 #ifdef STRONG_COMPARE
3085 // AWDEBUG
3086 OsiSolverInterface::unmarkHotStart();
3087 #endif
3088 strong_branching_solver_->unmarkHotStart(this);
3089 optimizationStatus_ = optimizationStatusBeforeHotStart_;
3090 }
3091 else {
3092 // Default Implementation
3093 OsiSolverInterface::unmarkHotStart();
3094 }
3095 }
3096
getObjCoefficients() const3097 const double * OsiTMINLPInterface::getObjCoefficients() const
3098 {
3099 const int n = getNumCols();
3100 delete [] obj_;
3101 obj_ = NULL;
3102 obj_ = new double[n];
3103
3104 bool new_x = true;
3105 const double* x_sol = problem_->x_sol();
3106 bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
3107
3108 if (!retval) {
3109 // Let's see if that happens - it will cause a crash
3110 fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in OsiTMINLPInterface::getObjCoefficients()\n");
3111 delete [] obj_;
3112 obj_ = NULL;
3113 }
3114
3115 return obj_;
3116 }
3117
3118 /** Sets the TMINLP2TNLP to be used by the interface.*/
3119 void
use(Ipopt::SmartPtr<TMINLP2TNLP> tminlp2tnlp)3120 OsiTMINLPInterface::use(Ipopt::SmartPtr<TMINLP2TNLP> tminlp2tnlp){
3121 problem_ = tminlp2tnlp;
3122 problem_to_optimize_ = GetRawPtr(problem_);
3123 feasibilityProblem_->use(GetRawPtr(tminlp2tnlp));}
3124
3125 }/** end namespace Bonmin*/
3126
3127