00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "BonminConfig.h"
00015
00016 #include "BonOsiTMINLPInterface.hpp"
00017 #include "BonTMINLP.hpp"
00018 #include "BonTMINLP2TNLP.hpp"
00019 #include "BonTNLP2FPNLP.hpp"
00020 #include "BonTMINLP2OsiLP.hpp"
00021 #include "BonTNLPSolver.hpp"
00022 #include "CoinTime.hpp"
00023 #include <climits>
00024 #include <string>
00025 #include <sstream>
00026 #include "BonAuxInfos.hpp"
00027
00028 #include "Ipopt/BonIpoptSolver.hpp"
00029 #include "Ipopt/BonIpoptWarmStart.hpp"
00030 #ifdef COIN_HAS_FILTERSQP
00031 #include "Filter/BonFilterSolver.hpp"
00032 #include "Filter/BonFilterWarmStart.hpp"
00033
00034 #endif
00035
00036 #include "OsiBranchingObject.hpp"
00037 #include "OsiRowCutDebugger.hpp"
00038 #include "BonStrongBranchingSolver.hpp"
00039
00040
00041 #include "BonMsgUtils.hpp"
00042
00043 using namespace Ipopt;
00044
00045
00046 extern bool BonminAbortAll;
00047 namespace Bonmin {
00049 static void
00050 register_general_options
00051 (SmartPtr<RegisteredOptions> roptions)
00052 {
00053 roptions->SetRegisteringCategory("NLP interface", RegisteredOptions::BonminCategory);
00054 roptions->AddStringOption3("nlp_solver",
00055 "Choice of the solver for local optima of continuous NLP's",
00056 "Ipopt",
00057 "Ipopt", "Interior Point OPTimizer (https://projects.coin-or.org/Ipopt)",
00058 "filterSQP", "Sequential quadratic programming trust region "
00059 "algorithm (http://www-unix.mcs.anl.gov/~leyffer/solvers.html)",
00060 "all", "run all available solvers at each node",
00061 "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.");
00062 roptions->setOptionExtraInfo("nlp_solver",127);
00063
00064 roptions->AddStringOption4("warm_start",
00065 "Select the warm start method",
00066 "none",
00067 "none","No warm start, just start NLPs from optimal solution of the root relaxation",
00068 "fake_basis", "builds fake basis, useful for cut management in Cbc (warm start is the same as in none)",
00069 "optimum","Warm start with direct parent optimum",
00070 "interior_point","Warm start with an interior point of direct parent",
00071 "This will affect the function getWarmStart(), and as a consequence the warm starting in the various algorithms.");
00072 roptions->setOptionExtraInfo("warm_start",8);
00073
00074 roptions->SetRegisteringCategory("Output and Loglevel", RegisteredOptions::BonminCategory);
00075
00076 roptions->AddBoundedIntegerOption("nlp_log_level",
00077 "specify NLP solver interface log level (independent from ipopt print_level).",
00078 0,2,1,
00079 "Set the level of output of the OsiTMINLPInterface : "
00080 "0 - none, 1 - normal, 2 - verbose"
00081 );
00082 roptions->setOptionExtraInfo("nlp_log_level",127);
00083
00084 roptions->AddStringOption2("file_solution",
00085 "Write a file bonmin.sol with the solution",
00086 "no",
00087 "yes","",
00088 "no","","");
00089 roptions->setOptionExtraInfo("file_solution",127);
00090
00091 roptions->SetRegisteringCategory("NLP solution robustness", RegisteredOptions::BonminCategory);
00092
00093 roptions->AddLowerBoundedNumberOption("max_random_point_radius",
00094 "Set max value r for coordinate of a random point.",
00095 0.,1,1e5,
00096 "When picking a random point, coordinate i will be in the interval [min(max(l,-r),u-r), max(min(u,r),l+r)] "
00097 "(where l is the lower bound for the variable and u is its upper bound)");
00098 roptions->setOptionExtraInfo("max_random_point_radius",8);
00099
00100 roptions->AddStringOption3("random_point_type","method to choose a random starting point",
00101 "Jon",
00102 "Jon", "Choose random point uniformly between the bounds",
00103 "Andreas", "perturb the starting point of the problem within a prescribed interval",
00104 "Claudia", "perturb the starting point using the perturbation radius suffix information",
00105 "");
00106 roptions->setOptionExtraInfo("random_point_type",8);
00107
00108 roptions->AddLowerBoundedNumberOption("random_point_perturbation_interval",
00109 "Amount by which starting point is perturbed when choosing to pick random point by perturbing starting point",
00110 0.,true, 1.,
00111 "");
00112 roptions->setOptionExtraInfo("random_point_perturbation_interval",8);
00113
00114
00115 roptions->AddLowerBoundedIntegerOption
00116 ("num_iterations_suspect",
00117 "Number of iterations over which a node is considered \"suspect\" (for debugging purposes only, see detailed documentation).",
00118 -1,-1,
00119 "When the number of iterations to solve a node is above this number, the subproblem at this"
00120 " node is considered to be suspect and it will be written into a file (set to -1 to deactivate this).");
00121 roptions->setOptionExtraInfo("num_iterations_suspect",127);
00122
00123
00124
00125 roptions->AddLowerBoundedIntegerOption("num_retry_unsolved_random_point",
00126 "Number $k$ of times that the algorithm will try to resolve an unsolved NLP with a random starting point "
00127 "(we call unsolved an NLP for which Ipopt is not "
00128 "able to guarantee optimality within the specified tolerances).",
00129 0,0,
00130 "When Ipopt fails to solve a continuous NLP sub-problem, if $k > 0$, the algorithm will "
00131 "try again to solve the failed NLP with $k$ new randomly chosen starting points "
00132 " or until the problem is solved with success.");
00133 roptions->setOptionExtraInfo("num_retry_unsolved_random_point",127);
00134
00135 roptions->AddLowerBoundedNumberOption("resolve_on_small_infeasibility",
00136 "If a locally infeasible problem is infeasible by less than this, resolve it "
00137 "with initial starting point.",
00138 0.,false, 0.,
00139 "It is set to 0 by default with Ipopt. "
00140 "When using FilterSQP, Bonmin sets it to a small value.");
00141 roptions->setOptionExtraInfo("resolve_on_small_infeasibility",8);
00142
00143
00144 roptions->SetRegisteringCategory("Nonconvex problems", RegisteredOptions::BonminCategory);
00145
00146
00147 roptions->AddLowerBoundedIntegerOption("num_resolve_at_root",
00148 "Number $k$ of tries to resolve the root node with different starting points.",
00149 0,0,
00150 "The algorithm will solve the root node with $k$ random starting points"
00151 " and will keep the best local optimum found.");
00152 roptions->setOptionExtraInfo("num_resolve_at_root",8);
00153
00154 roptions->AddLowerBoundedIntegerOption("num_resolve_at_node",
00155 "Number $k$ of tries to resolve a node (other than the root) of the tree with different starting point.",
00156 0,0,
00157 "The algorithm will solve all the nodes with $k$ different random starting points "
00158 "and will keep the best local optimum found.");
00159 roptions->setOptionExtraInfo("num_resolve_at_node",8);
00160
00161 roptions->AddLowerBoundedIntegerOption("num_resolve_at_infeasibles",
00162 "Number $k$ of tries to resolve an infeasible node (other than the root) of the tree with different starting point.",
00163 0,0,
00164 "The algorithm will solve all the infeasible nodes with $k$ different random starting points "
00165 "and will keep the best local optimum found.");
00166 roptions->setOptionExtraInfo("num_resolve_at_infeasibles",8);
00167
00168
00169 roptions->AddStringOption2("dynamic_def_cutoff_decr",
00170 "Do you want to define the parameter cutoff_decr dynamically?",
00171 "no",
00172 "no", "",
00173 "yes", "");
00174 roptions->setOptionExtraInfo("dynamic_def_cutoff_decr",8);
00175
00176 roptions->AddLowerBoundedNumberOption("coeff_var_threshold",
00177 "Coefficient of variation threshold (for dynamic definition of cutoff_decr).",
00178 0.0,
00179 false,
00180 0.1);
00181 roptions->setOptionExtraInfo("coeff_var_threshold",8);
00182
00183 roptions->AddNumberOption("first_perc_for_cutoff_decr",
00184 "The percentage used when, the coeff of variance is smaller than the threshold, to compute the cutoff_decr dynamically.",
00185 -0.02);
00186 roptions->setOptionExtraInfo("first_perc_for_cutoff_decr",8);
00187
00188 roptions->AddNumberOption("second_perc_for_cutoff_decr",
00189 "The percentage used when, the coeff of variance is greater than the threshold, to compute the cutoff_decr dynamically.",
00190 -0.05);
00191 roptions->setOptionExtraInfo("second_perc_for_cutoff_decr",8);
00192
00193 }
00194
00195 static void register_OA_options
00196 (SmartPtr<RegisteredOptions> roptions)
00197 {
00198
00199
00200 roptions->SetRegisteringCategory("Outer Approximation strengthening", RegisteredOptions::UndocumentedCategory);
00201 roptions->AddStringOption2("disjunctive_cut_type",
00202 "Determine if and what kind of disjunctive cuts should be computed.",
00203 "none",
00204 "none", "No disjunctive cuts.",
00205 "most-fractional", "If discrete variables present, compute disjunction for most-fractional variable");
00206 roptions->setOptionExtraInfo("disjunctive_cut_type",119);
00207
00208
00209
00210
00211 roptions->AddStringOption4("cut_strengthening_type",
00212 "Determines if and what kind of cut strengthening should be performed.",
00213 "none",
00214 "none", "No strengthening of cuts.",
00215 "sglobal", "Strengthen global cuts.",
00216 "uglobal-slocal", "Unstrengthened global and strengthened local cuts",
00217 "sglobal-slocal", "Strengthened global and strengthened local cuts",
00218 "");
00219 roptions->setOptionExtraInfo("cut_strengthening_type",119);
00220
00221 roptions->SetRegisteringCategory("Outer Approximation cuts generation", RegisteredOptions::BonminCategory);
00222
00223 roptions->AddStringOption2("oa_cuts_scope","Specify if OA cuts added are to be set globally or locally valid",
00224 "global",
00225 "local","Cuts are treated as locally valid",
00226 "global", "Cuts are treated as globally valid",
00227 "");
00228 roptions->setOptionExtraInfo("oa_cuts_scope",119);
00229
00230 roptions->AddStringOption2("add_only_violated_oa","Do we add all OA cuts or only the ones violated by current point?",
00231 "no",
00232 "no","Add all cuts",
00233 "yes","Add only violated cuts","");
00234 roptions->setOptionExtraInfo("add_only_violated_oa",119);
00235
00236
00237 roptions->AddLowerBoundedNumberOption("tiny_element","Value for tiny element in OA cut",
00238 -0.,0,1e-08,
00239 "We will remove \"cleanly\" (by relaxing cut) an element lower"
00240 " than this.");
00241 roptions->setOptionExtraInfo("tiny_element",119);
00242
00243 roptions->AddLowerBoundedNumberOption("very_tiny_element","Value for very tiny element in OA cut",
00244 -0.,0,1e-17,
00245 "Algorithm will take the risk of neglecting an element lower"
00246 " than this.");
00247 roptions->setOptionExtraInfo("very_tiny_element",119);
00248
00249 roptions->AddLowerBoundedNumberOption("oa_rhs_relax","Value by which to relax OA cut",
00250 -0.,0,1e-8,
00251 "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)."
00252 );
00253 roptions->setOptionExtraInfo("oa_rhs_relax",119);
00254
00255 roptions->SetRegisteringCategory("Output and Loglevel", RegisteredOptions::BonminCategory);
00256
00257 roptions->AddLowerBoundedIntegerOption("oa_cuts_log_level",
00258 "level of log when generating OA cuts.",
00259 0, 0,
00260 "0: outputs nothing,\n"
00261 "1: when a cut is generated, its violation and index of row from which it originates,\n"
00262 "2: always output violation of the cut.\n"
00263 "3: output generated cuts incidence vectors.");
00264 roptions->setOptionExtraInfo("oa_cuts_log_level",119);
00265
00266 }
00267
00268
00270 void
00271 OsiTMINLPInterface::registerOptions
00272 (SmartPtr<RegisteredOptions> roptions)
00273 {
00274
00275
00276 try {
00277 register_general_options(roptions);
00278 register_OA_options(roptions);
00279 #ifdef COIN_HAS_FILTERSQP
00280 FilterSolver::RegisterOptions(roptions);
00281 #endif
00282 IpoptSolver::RegisterOptions(roptions);
00283 }
00284 catch(RegisteredOptions::OPTION_ALREADY_REGISTERED) {
00285
00286 }
00287
00288
00289 }
00290
00292 void
00293 OsiTMINLPInterface::setAppDefaultOptions(SmartPtr<OptionsList> Options)
00294 {}
00295
00296 OsiTMINLPInterface::Messages::Messages
00297 ():CoinMessages((int)OSITMINLPINTERFACE_DUMMY_END)
00298 {
00299 strcpy(source_ ,"NLP");
00300 ADD_MSG(SOLUTION_FOUND, std_m, 2,
00301 "After %d tries found a solution of %g (previous best %g).");
00302
00303 ADD_MSG(INFEASIBLE_SOLUTION_FOUND, std_m, 2,
00304 "After %d tries found an solution of %g infeasible problem.");
00305
00306 ADD_MSG(UNSOLVED_PROBLEM_FOUND, std_m, 2,
00307 "After %d tries found an solution of %g unsolved problem.");
00308 ADD_MSG(WARN_SUCCESS_WS, warn_m, 2,
00309 "Problem not solved with warm start but solved without");
00310
00311 ADD_MSG(WARNING_RESOLVING, warn_m,2,
00312 "Trying to resolve NLP with different starting point (%d attempts).");
00313 ADD_MSG(WARN_SUCCESS_RANDOM, warn_m, 1,
00314 "Problem initially not solved but solved with a random starting point (success on %d attempt)");
00315 ADD_MSG(WARN_CONTINUING_ON_FAILURE, warn_m, 1,
00316 "Warning : continuing branching, while there are unrecovered failures at nodes");
00317
00318 ADD_MSG(SUSPECT_PROBLEM, std_m, 2, "NLP number %d is suspect (see bounds and start file)");
00319 ADD_MSG(IPOPT_SUMMARY, std_m, 2, "Ipopt return (for %s): status %2d, iter count %4d, time %g");
00320 ADD_MSG(BETTER_SOL, std_m, 2, "Solution of value %g found on %d'th attempt");
00321
00322 ADD_MSG(LOG_HEAD, std_m, 1,
00323 "\n "
00324 " Num Status Obj It time Location");
00325 ADD_MSG(LOG_LINE, std_m, 1,
00326 "%c %8d %11s %g %8d %g %20s");
00327
00328 ADD_MSG(ALTERNATE_OBJECTIVE, std_m, 1, "Objective value recomputed with alternate objective: %g.");
00329
00330 ADD_MSG(WARN_RESOLVE_BEFORE_INITIAL_SOLVE, warn_m, 1,
00331 "resolve called before any call to initialSol can not use warm starts.");
00332 ADD_MSG(ERROR_NO_TNLPSOLVER, warn_m, 1,"Can not parse options when no IpApplication has been created");
00333 ADD_MSG(WARNING_NON_CONVEX_OA, warn_m, 1,
00334 "OA on non-convex constraint is very experimental.");
00335 ADD_MSG(SOLVER_DISAGREE_STATUS, warn_m, 1, "%s says problem %s, %s says %s.");
00336 ADD_MSG(SOLVER_DISAGREE_VALUE, warn_m, 1, "%s gives objective %.16g, %s gives %.16g.");
00337
00338 }
00339
00340
00341 void
00342 OsiTMINLPInterface::OaMessageHandler::print(OsiRowCut &row){
00343 FILE * fp = filePointer();
00344 const int &n = row.row().getNumElements();
00345 fprintf(fp,"Row cut has %d elements. Lower bound: %g, upper bound %g.\n", n, row.lb(), row.ub());
00346 const int * idx = row.row().getIndices();
00347 const double * val = row.row().getElements();
00348 for(int i = 0 ; i < n ; i++){
00349 fprintf(fp,"%g, x%d",val[i], idx[i]);
00350 if(i && i % 7 == 0)
00351 fprintf(fp,"\n");
00352 }
00353 }
00354
00355 OsiTMINLPInterface::OaMessages::OaMessages(): CoinMessages((int) OA_MESSAGES_DUMMY_END){
00356 strcpy(source_,"OaCg");
00357 ADD_MSG(VIOLATED_OA_CUT_GENERATED, std_m, 1,"Row %d, cut violation is %g: Outer approximation cut generated.");
00358
00359 ADD_MSG(CUT_NOT_VIOLATED_ENOUGH, std_m, 2,"Row %d, cut violation is %g: Outer approximation cut not generated.");
00360
00361 ADD_MSG(OA_CUT_GENERATED, std_m, 1,"Row %d: Outer approximation cut not generated.");
00362 }
00363 bool OsiTMINLPInterface::hasPrintedOptions=0;
00364
00366
00369 OsiTMINLPInterface::OsiTMINLPInterface():
00370 OsiSolverInterface(),
00371 tminlp_(NULL),
00372 problem_(NULL),
00373 problem_to_optimize_(NULL),
00374 feasibility_mode_(false),
00375 app_(NULL),
00376 debug_apps_(),
00377 testOthers_(false),
00378 warmstart_(NULL),
00379 rowsense_(NULL),
00380 rhs_(NULL),
00381 rowrange_(NULL),
00382 reducedCosts_(NULL),
00383 OsiDualObjectiveLimit_(1e200),
00384 hasVarNamesFile_(true),
00385 nCallOptimizeTNLP_(0),
00386 totalNlpSolveTime_(0),
00387 totalIterations_(0),
00388 maxRandomRadius_(1e08),
00389 randomGenerationType_(0),
00390 max_perturbation_(COIN_DBL_MAX),
00391 pushValue_(1e-02),
00392 numRetryInitial_(-1),
00393 numRetryResolve_(-1),
00394 numRetryInfeasibles_(-1),
00395 numRetryUnsolved_(1),
00396 infeasibility_epsilon_(0),
00397 dynamicCutOff_(0),
00398 coeff_var_threshold_(0.1),
00399 first_perc_for_cutoff_decr_(-0.02),
00400 second_perc_for_cutoff_decr_(-0.05),
00401 messages_(),
00402 pretendFailIsInfeasible_(0),
00403 pretendSucceededNext_(0),
00404 hasContinuedAfterNlpFailure_(false),
00405 numIterationSuspect_(-1),
00406 hasBeenOptimized_(false),
00407 obj_(NULL),
00408 feasibilityProblem_(NULL),
00409 jRow_(NULL),
00410 jCol_(NULL),
00411 jValues_(NULL),
00412 nnz_jac(0),
00413 constTypes_(NULL),
00414 nNonLinear_(0),
00415 tiny_(1e-8),
00416 veryTiny_(1e-20),
00417 rhsRelax_(tiny_),
00418 infty_(1e100),
00419 warmStartMode_(None),
00420 firstSolve_(true),
00421 cutStrengthener_(NULL),
00422 oaMessages_(),
00423 oaHandler_(NULL),
00424 newCutoffDecr(COIN_DBL_MAX)
00425
00426 {
00427 oaHandler_ = new OaMessageHandler;
00428 oaHandler_->setLogLevel(0);
00429 }
00430
00431 void
00432 OsiTMINLPInterface::initialize(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00433 Ipopt::SmartPtr<Ipopt::OptionsList> options,
00434 Ipopt::SmartPtr<Ipopt::Journalist> journalist,
00435 const std::string & prefix,
00436 Ipopt::SmartPtr<TMINLP> tminlp){
00437 if(!IsValid(app_))
00438 createApplication(roptions, options, journalist, prefix);
00439 setModel(tminlp);
00440 }
00441
00442 void OsiTMINLPInterface::setSolver(Ipopt::SmartPtr<TNLPSolver> app){
00443 app_ = app;
00444 }
00445
00446
00447 void
00448 OsiTMINLPInterface::createApplication(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00449 Ipopt::SmartPtr<Ipopt::OptionsList> options,
00450 Ipopt::SmartPtr<Ipopt::Journalist> journalist,
00451 const std::string & prefix
00452 )
00453 {
00454 assert(!IsValid(app_));
00455 int ival;
00456 options->GetEnumValue("nlp_solver", ival, prefix.c_str());
00457 Solver s = (Solver) ival;
00458 if(s == EFilterSQP){
00459 testOthers_ = false;;
00460 #ifdef COIN_HAS_FILTERSQP
00461 app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
00462 #else
00463 throw SimpleError("createApplication",
00464 "Bonmin not configured to run with FilterSQP.");
00465 #endif
00466
00467 debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix));
00468 }
00469 else if(s == EIpopt){
00470 testOthers_ = false;
00471 app_ = new IpoptSolver(roptions, options, journalist, prefix);
00472
00473 #ifdef COIN_HAS_FILTERSQP
00474 debug_apps_.push_back(new Bonmin::FilterSolver(roptions, options, journalist, prefix));
00475 #endif
00476 }
00477 else if(s == EAll){
00478 #ifdef COIN_HAS_FILTERSQP
00479 app_ = new Bonmin::FilterSolver(roptions, options, journalist, prefix);
00480 #else
00481 throw SimpleError("createApplication",
00482 "Bonmin not configured to run with Ipopt.");
00483 #endif
00484 debug_apps_.push_back(new IpoptSolver(roptions, options, journalist, prefix));
00485 testOthers_ = true;
00486 }
00487 if (!app_->Initialize("")) {
00488 throw CoinError("Error during initialization of app_","createApplication", "OsiTMINLPInterface");
00489 }
00490 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin() ;
00491 i != debug_apps_.end() ; i++){
00492 (*i)->Initialize("");
00493 }
00494 extractInterfaceParams();
00495
00496 }
00497
00499 void
00500 OsiTMINLPInterface::setModel(SmartPtr<TMINLP> tminlp)
00501 {
00502 assert(IsValid(tminlp));
00503 tminlp_ = tminlp;
00504 problem_ = new TMINLP2TNLP(tminlp_);
00505 feasibilityProblem_ = new TNLP2FPNLP
00506 (SmartPtr<TNLP>(GetRawPtr(problem_)));
00507 if(feasibility_mode_){
00508 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
00509 }
00510 else {
00511 problem_to_optimize_ = GetRawPtr(problem_);
00512 }
00513 }
00514
00515
00516
00517 void
00518 OsiTMINLPInterface::readOptionFile(const std::string & fileName)
00519 {
00520 if(IsValid(app_)){
00521 std::ifstream is;
00522 if (fileName != "") {
00523 try {
00524 is.open(fileName.c_str());
00525 }
00526 catch(std::bad_alloc) {
00527 throw CoinError("Not enough memory to open option file.\n", "readOptionFile", "OsiTMINLPInterface");
00528 }
00529 }
00530 options()->ReadFromStream(*app_->journalist(), is);
00531 extractInterfaceParams();
00532 }
00533 }
00534
00536 OsiTMINLPInterface::OsiTMINLPInterface (const OsiTMINLPInterface &source):
00537 OsiSolverInterface(source),
00538 tminlp_(source.tminlp_),
00539 problem_(NULL),
00540 problem_to_optimize_(NULL),
00541 feasibility_mode_(source.feasibility_mode_),
00542 rowsense_(NULL),
00543 rhs_(NULL),
00544 rowrange_(NULL),
00545 reducedCosts_(NULL),
00546 OsiDualObjectiveLimit_(source.OsiDualObjectiveLimit_),
00547 hasVarNamesFile_(source.hasVarNamesFile_),
00548 nCallOptimizeTNLP_(0),
00549 totalNlpSolveTime_(0),
00550 totalIterations_(0),
00551 maxRandomRadius_(source.maxRandomRadius_),
00552 randomGenerationType_(source.randomGenerationType_),
00553 max_perturbation_(source.max_perturbation_),
00554 pushValue_(source.pushValue_),
00555 numRetryInitial_(source.numRetryInitial_),
00556 numRetryResolve_(source.numRetryResolve_),
00557 numRetryInfeasibles_(source.numRetryInfeasibles_),
00558 numRetryUnsolved_(source.numRetryUnsolved_),
00559 infeasibility_epsilon_(source.infeasibility_epsilon_),
00560 dynamicCutOff_(source.dynamicCutOff_),
00561 coeff_var_threshold_(source.coeff_var_threshold_),
00562 first_perc_for_cutoff_decr_(source.first_perc_for_cutoff_decr_),
00563 second_perc_for_cutoff_decr_(source.second_perc_for_cutoff_decr_),
00564 messages_(),
00565 pretendFailIsInfeasible_(source.pretendFailIsInfeasible_),
00566 pretendSucceededNext_(source.pretendSucceededNext_),
00567 hasContinuedAfterNlpFailure_(source.hasContinuedAfterNlpFailure_),
00568 numIterationSuspect_(source.numIterationSuspect_),
00569 hasBeenOptimized_(source.hasBeenOptimized_),
00570 obj_(NULL),
00571 feasibilityProblem_(NULL),
00572 jRow_(NULL),
00573 jCol_(NULL),
00574 jValues_(NULL),
00575 nnz_jac(source.nnz_jac),
00576 constTypes_(NULL),
00577
00578 nNonLinear_(0),
00579 tiny_(source.tiny_),
00580 veryTiny_(source.veryTiny_),
00581 rhsRelax_(source.rhsRelax_),
00582 infty_(source.infty_),
00583 warmStartMode_(source.warmStartMode_),
00584 firstSolve_(true),
00585 cutStrengthener_(source.cutStrengthener_),
00586 oaMessages_(),
00587 oaHandler_(NULL),
00588 newCutoffDecr(source.newCutoffDecr),
00589 strong_branching_solver_(source.strong_branching_solver_)
00590 {
00591 if(IsValid(source.tminlp_)) {
00592 problem_ = source.problem_->clone();
00593 feasibilityProblem_ = new TNLP2FPNLP
00594 (SmartPtr<TNLP>(GetRawPtr(problem_)), source.feasibilityProblem_);
00595 if(feasibility_mode_)
00596 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
00597 else
00598 problem_to_optimize_ = GetRawPtr(problem_);
00599 pretendFailIsInfeasible_ = source.pretendFailIsInfeasible_;
00600 pretendSucceededNext_ = source.pretendSucceededNext_;
00601
00602 setAuxiliaryInfo(source.getAuxiliaryInfo());
00603
00604 app_ = source.app_->clone();
00605 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::const_iterator i = source.debug_apps_.begin();
00606 i != source.debug_apps_.end() ; i++){
00607 debug_apps_.push_back((*i)->clone());
00608 }
00609 testOthers_ = source.testOthers_;
00610 }
00611 else {
00612 throw SimpleError("Don't know how to copy an empty IpoptInterface.",
00613 "copy constructor");
00614 }
00615
00616 warmstart_ = source.warmstart_ ? source.warmstart_->clone() : NULL;
00617
00618 if(source.obj_) {
00619 obj_ = new double[source.getNumCols()];
00620 CoinCopyN(source.obj_, source.getNumCols(), obj_);
00621 }
00622
00623
00624 oaHandler_ = new OaMessageHandler(*source.oaHandler_);;
00625
00626 }
00627
00628 OsiSolverInterface *
00629 OsiTMINLPInterface::clone(bool copyData ) const
00630 {
00631 if(copyData)
00632 return new OsiTMINLPInterface(*this);
00633 else return new OsiTMINLPInterface;
00634 }
00635
00637 OsiTMINLPInterface & OsiTMINLPInterface::operator=(const OsiTMINLPInterface& rhs)
00638 {
00639 if(this!= &rhs) {
00640 OsiSolverInterface::operator=(rhs);
00641 OsiDualObjectiveLimit_ = rhs.OsiDualObjectiveLimit_;
00642 nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
00643 totalNlpSolveTime_ = rhs.nCallOptimizeTNLP_;
00644 totalIterations_ = rhs.totalIterations_;
00645 maxRandomRadius_ = rhs.maxRandomRadius_;
00646 hasVarNamesFile_ = rhs.hasVarNamesFile_;
00647 pushValue_ = rhs.pushValue_;
00648
00649 delete warmstart_;
00650 warmstart_ = NULL;
00651
00652 if(IsValid(rhs.tminlp_)) {
00653
00654 tminlp_ = rhs.tminlp_;
00655 problem_ = new TMINLP2TNLP(tminlp_);
00656 problem_to_optimize_ = GetRawPtr(problem_);
00657 app_ = rhs.app_->clone();
00658
00659 warmstart_ = rhs.warmstart_ ? rhs.warmstart_->clone() : NULL;
00660
00661 feasibilityProblem_ = new TNLP2FPNLP
00662 (SmartPtr<TNLP>(GetRawPtr(problem_)));
00663 nnz_jac = rhs.nnz_jac;
00664
00665 if(constTypes_ != NULL) {
00666 delete [] constTypes_;
00667 constTypes_ = NULL;
00668 }
00669 if(rhs.constTypes_ != NULL) {
00670 constTypes_ = new TNLP::LinearityType[getNumRows()];
00671 CoinCopyN(rhs.constTypes_, getNumRows(), constTypes_);
00672 }
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683 if(rhs.jValues_!=NULL && rhs.jRow_ != NULL && rhs.jCol_ != NULL && nnz_jac>0) {
00684 jValues_ = new double [nnz_jac];
00685 jCol_ = new Index [nnz_jac];
00686 jRow_ = new Index [nnz_jac];
00687 CoinCopyN(rhs.jValues_ , nnz_jac,jValues_ );
00688 CoinCopyN(rhs.jCol_ , nnz_jac,jCol_ );
00689 CoinCopyN(rhs.jRow_ , nnz_jac,jRow_ );
00690 }
00691 else if(nnz_jac > 0) {
00692 throw CoinError("Arrays for storing jacobian are inconsistant.",
00693 "copy constructor",
00694 "IpoptOAInterface");
00695 }
00696 tiny_ = rhs.tiny_;
00697 veryTiny_ = rhs.veryTiny_;
00698 rhsRelax_ = rhs.rhsRelax_;
00699 infty_ = rhs.infty_;
00700 warmStartMode_ = rhs.warmStartMode_;
00701 newCutoffDecr = rhs.newCutoffDecr;
00702
00703 }
00704 else {
00705 tminlp_ =NULL;
00706 problem_ = NULL;
00707 app_ = NULL;
00708 feasibilityProblem_ = NULL;
00709 }
00710
00711
00712 if(obj_) {
00713 delete [] obj_;
00714 obj_ = NULL;
00715 }
00716 if(rhs.obj_) {
00717 obj_ = new double[rhs.getNumCols()];
00718 CoinCopyN(rhs.obj_, rhs.getNumCols(), obj_);
00719 }
00720
00721 hasVarNamesFile_ = rhs.hasVarNamesFile_;
00722
00723 nCallOptimizeTNLP_ = rhs.nCallOptimizeTNLP_;
00724 totalNlpSolveTime_ = rhs.totalNlpSolveTime_;
00725 totalIterations_ = rhs.totalIterations_;
00726 maxRandomRadius_ = rhs.maxRandomRadius_;
00727 pushValue_ = rhs.pushValue_;
00728 numRetryInitial_ = rhs.numRetryInitial_;
00729 numRetryResolve_ = rhs.numRetryResolve_;
00730 numRetryInfeasibles_ = rhs.numRetryInfeasibles_;
00731 numRetryUnsolved_ = rhs.numRetryUnsolved_;
00732 infeasibility_epsilon_ = rhs.infeasibility_epsilon_;
00733 pretendFailIsInfeasible_ = rhs.pretendFailIsInfeasible_;
00734 pretendSucceededNext_ = rhs.pretendSucceededNext_;
00735 numIterationSuspect_ = rhs.numIterationSuspect_;
00736
00737 hasBeenOptimized_ = rhs.hasBeenOptimized_;
00738 cutStrengthener_ = rhs.cutStrengthener_;
00739
00740 delete oaHandler_;
00741 oaHandler_ = new OaMessageHandler(*rhs.oaHandler_);
00742 strong_branching_solver_ = rhs.strong_branching_solver_;
00743
00744 dynamicCutOff_ = rhs.dynamicCutOff_;
00745 coeff_var_threshold_ = rhs.coeff_var_threshold_;
00746 first_perc_for_cutoff_decr_ = rhs.first_perc_for_cutoff_decr_;
00747 second_perc_for_cutoff_decr_ = rhs.second_perc_for_cutoff_decr_;
00748
00749 freeCachedData();
00750 }
00751 return *this;
00752 }
00753
00754 const SmartPtr<OptionsList> OsiTMINLPInterface::options() const
00755 {
00756 if(!IsValid(app_)) {
00757 messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
00758 return NULL;
00759 }
00760 else
00761 return app_->options();
00762 }
00763
00764 SmartPtr<OptionsList> OsiTMINLPInterface::options()
00765 {
00766 if(!IsValid(app_)) {
00767 messageHandler()->message(ERROR_NO_TNLPSOLVER, messages_)<<CoinMessageEol;
00768 return NULL;
00769 }
00770 else
00771 return app_->options();
00772 }
00773
00775 OsiTMINLPInterface::~OsiTMINLPInterface ()
00776 {
00777 freeCachedData();
00778 delete [] jRow_;
00779 delete [] jCol_;
00780 delete [] jValues_;
00781 delete [] constTypes_;
00782 delete [] obj_;
00783 delete oaHandler_;
00784 delete warmstart_;
00785 }
00786
00787 void
00788 OsiTMINLPInterface::freeCachedColRim()
00789 {
00790 if(reducedCosts_!=NULL) {
00791 delete [] reducedCosts_;
00792 reducedCosts_ = NULL;
00793 }
00794
00795 }
00796
00797 void
00798 OsiTMINLPInterface::freeCachedRowRim()
00799 {
00800 if(rowsense_!=NULL) {
00801 delete [] rowsense_;
00802 rowsense_ = NULL;
00803 }
00804 if(rhs_!=NULL) {
00805 delete [] rhs_;
00806 rhs_ = NULL;
00807 }
00808 if(rowrange_!=NULL) {
00809 delete [] rowrange_;
00810 rowrange_ = NULL;
00811 }
00812
00813
00814
00815
00816 }
00817
00818 void
00819 OsiTMINLPInterface::freeCachedData()
00820 {
00821 freeCachedColRim();
00822 freeCachedRowRim();
00823 }
00824
00825 const char * OsiTMINLPInterface::OPT_SYMB="OPT";
00826 const char * OsiTMINLPInterface::FAILED_SYMB="FAILED";
00827 const char * OsiTMINLPInterface::UNBOUND_SYMB="UNBOUNDED";
00828 const char * OsiTMINLPInterface::INFEAS_SYMB="INFEAS";
00829 const char * OsiTMINLPInterface::TIME_SYMB="TIME";
00831
00833
00834
00835 void
00836 OsiTMINLPInterface::resolveForCost(int numsolve, bool keepWarmStart)
00837 {
00838
00839
00840
00841 delete warmstart_;
00842 warmstart_ = NULL;
00843
00844 double * of_current = (numsolve > 0) ? new double[numsolve]: NULL;
00845 int num_failed, num_infeas;
00846 double mean, std_dev, var_coeff;
00847 double min = DBL_MAX;
00848 double max = -DBL_MAX;
00849
00850 Coin::SmartPtr<SimpleReferencedPtr<CoinWarmStart> > ws_backup = NULL;
00851 if(warmStartMode_ <= Optimum && keepWarmStart){
00852
00853
00854 ws_backup = make_referenced(app_->getUsedWarmStart(problem_));
00855 }
00857 vector<double> point(getNumCols()*3+ getNumRows());
00858 double bestBound = (isProvenOptimal()) ? getObjValue() : DBL_MAX;
00859 CoinCopyN(getColSolution(),
00860 getNumCols(), point());
00861 CoinCopyN(getRowPrice(),
00862 2*getNumCols()+ getNumRows(),
00863 point() + getNumCols());
00864 TNLPSolver::ReturnStatus savedStatus = optimizationStatus_;
00865 if(isProvenOptimal())
00866 messageHandler()->message(SOLUTION_FOUND,
00867 messages_)
00868 <<1<<getObjValue()<<bestBound
00869 <<CoinMessageEol;
00870 else
00871 messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
00872 messages_)
00873 <<1
00874 <<CoinMessageEol;
00875
00876 num_failed = 0;
00877 num_infeas = 0;
00878 mean = 0;
00879
00880 for(int f = 0; f < numsolve ; f++) {
00881 messageHandler()->message(WARNING_RESOLVING,
00882 messages_)
00883 <<f+1<< CoinMessageEol ;
00884 randomStartingPoint();
00885 solveAndCheckErrors(0,0,"resolve cost");
00886
00887
00888 char c=' ';
00889
00890 if(isProvenOptimal() &&
00891 getObjValue()<bestBound) {
00892 c='*';
00893 messageHandler()->message(BETTER_SOL, messages_)<<getObjValue()<<f+1<< CoinMessageEol;
00894 CoinCopyN(getColSolution(),
00895 getNumCols(), point());
00896 CoinCopyN(getRowPrice(),
00897 2*getNumCols()+ getNumRows(),
00898 point() + getNumCols());
00899 bestBound = getObjValue();
00900 savedStatus = optimizationStatus_;
00901 }
00902
00903 messageHandler()->message(LOG_LINE, messages_)
00904 <<c<<f+1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<app_->CPUTime()<<"resolve cost"<<CoinMessageEol;
00905
00906 if(isAbandoned()) {
00907 num_failed++;
00908 }
00909 else if(isProvenPrimalInfeasible()) {
00910 num_infeas++;
00911 }
00912
00913 else if(isProvenOptimal())
00914 messageHandler()->message(SOLUTION_FOUND,
00915 messages_)
00916 <<f+2<<getObjValue()<<bestBound
00917 <<CoinMessageEol;
00918 else if(!isAbandoned())
00919 messageHandler()->message(UNSOLVED_PROBLEM_FOUND,
00920 messages_)
00921 <<f+2
00922 <<CoinMessageEol;
00923 else
00924 messageHandler()->message(INFEASIBLE_SOLUTION_FOUND,
00925 messages_)
00926 <<f+2
00927 <<CoinMessageEol;
00928
00929 if(of_current != NULL){
00930 if(isProvenOptimal())
00931 {
00932 of_current[f] = getObjValue();
00933 mean=mean+of_current[f];
00934 if (of_current[f] < min)
00935 min = of_current[f];
00936 else if (of_current[f] > max)
00937 max = of_current[f];
00938 }
00939 else
00940 {
00941 of_current[f] = 0;
00942 }
00943 }
00944 }
00945
00946
00947 if(of_current != NULL){
00948
00949 mean=mean/(numsolve-num_failed-num_infeas);
00950
00951 std_dev = 0;
00952
00953
00954 for(int i=0; i<numsolve; i++)
00955 {
00956 if(of_current[i]!=0)
00957 std_dev=std_dev+pow(of_current[i]-mean,2);
00958 }
00959 std_dev=pow((std_dev/(numsolve-num_failed-num_infeas)),0.5);
00960
00961
00962 var_coeff=std_dev/mean;
00963
00964 if(dynamicCutOff_)
00965 {
00966 if(var_coeff<0.1)
00967 {
00968 setNewCutoffDecr(mean*first_perc_for_cutoff_decr_);
00969 }
00970 else
00971 {
00972 setNewCutoffDecr(mean*second_perc_for_cutoff_decr_);
00973 }
00974 }
00975 }
00976
00977 problem_->Set_x_sol(getNumCols(),point());
00978 problem_->Set_dual_sol((int) point.size()-getNumCols(), point() + getNumCols());
00979 problem_->set_obj_value(bestBound);
00980 optimizationStatus_ = savedStatus;
00981 hasBeenOptimized_ = true;
00982
00983 if(warmStartMode_ < Optimum && keepWarmStart) {
00984 app_->setWarmStart(ws_backup->ptr(), problem_);
00985 }
00986 }
00987
00988 void
00989 OsiTMINLPInterface::resolveForRobustness(int numsolve)
00990 {
00991
00992
00993
00994 delete warmstart_;
00995 warmstart_ = NULL;
00996
00997
00998 CoinWarmStart * ws_backup = NULL;
00999 if(warmStartMode_ < Optimum){
01000
01001
01002 ws_backup = app_->getUsedWarmStart(problem_);
01003 }
01004
01005
01006 app_->disableWarmStart();
01007 problem()->resetStartingPoint();
01008 messageHandler()->message(WARNING_RESOLVING,
01009 messages_)
01010 <<1<< CoinMessageEol ;
01011 solveAndCheckErrors(0,0,"resolve robustness");
01012
01013
01014 char c='*';
01015 if(isAbandoned()) {
01016 c=' ';
01017 }
01018 messageHandler()->message(LOG_LINE, messages_)
01019 <<c<<1<<statusAsString()<<getObjValue()<<app_->IterationCount()<<
01020 app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
01021
01022
01023 if(!isAbandoned()) {
01024 messageHandler()->message(WARN_SUCCESS_WS, messages_) << CoinMessageEol ;
01025
01026 app_->enableWarmStart();
01027 if (warmStartMode_ < Optimum) {
01028 app_->setWarmStart(ws_backup, problem_);
01029 delete ws_backup;
01030 }
01031 return;
01032 }
01033
01034
01035 for(int f = 0; f < numsolve ; f++) {
01036 messageHandler()->message(WARNING_RESOLVING,
01037 messages_)
01038 <<f+2<< CoinMessageEol ;
01039
01040 randomStartingPoint();
01041 solveAndCheckErrors(0,0,"resolve robustness");
01042
01043
01044 messageHandler()->message(IPOPT_SUMMARY, messages_)
01045 <<"resolveForRobustness"<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
01046
01047
01048 char c='*';
01049 if(isAbandoned()) {
01050 c=' ';
01051 }
01052 messageHandler()->message(LOG_LINE, messages_)
01053 <<c<<f+2<<statusAsString()<<getObjValue()
01054 <<app_->IterationCount()<<app_->CPUTime()<<"resolve robustness"<<CoinMessageEol;
01055
01056
01057 if(!isAbandoned()) {
01058 messageHandler()->message(WARN_SUCCESS_RANDOM, messages_)
01059 << f+2 << CoinMessageEol ;
01060
01061 app_->enableWarmStart();
01062 if (warmStartMode_ < Optimum) {
01063 app_->setWarmStart(ws_backup, problem_);
01064 delete ws_backup;
01065 }
01066
01067 return;
01068 }
01069 }
01070
01071
01072 if(warmStartMode_ < Optimum){
01073 app_->setWarmStart(ws_backup, problem_);
01074 delete ws_backup;
01075 }
01076 if(pretendFailIsInfeasible_) {
01077 pretendSucceededNext_ = true;
01078 if(pretendFailIsInfeasible_ == 1) {
01079 messageHandler()->message(WARN_CONTINUING_ON_FAILURE,
01080 messages_)
01081 <<CoinMessageEol;
01082 hasContinuedAfterNlpFailure_ = 1;
01083 }
01084 return;
01085 }
01086 else {
01087 std::string probName;
01088 getStrParam(OsiProbName,probName);
01089 throw newUnsolvedError(app_->errorCode(), problem_,
01090 probName);
01091 }
01092
01093 }
01094
01096
01099 int OsiTMINLPInterface::getNumCols() const
01100 {
01101
01102 return problem_->num_variables();
01103 }
01104
01105
01107 int
01108 OsiTMINLPInterface::getNumRows() const
01109 {
01110 return problem_->num_constraints();
01111 }
01112
01113 const double *
01114 OsiTMINLPInterface::getColLower() const
01115 {
01116 return problem_->x_l();
01117 }
01118
01119 const double *
01120 OsiTMINLPInterface::getColUpper() const
01121 {
01122 return problem_->x_u();
01123 }
01124
01125 #if 1
01126
01127
01129 const OsiSolverInterface::OsiNameVec&
01130 OsiTMINLPInterface::getVarNames() {
01131 return getColNames();
01132 }
01133 #endif
01134
01135
01136 void OsiTMINLPInterface::extractSenseRhsAndRange() const
01137 {
01138 assert(rowsense_==NULL&&rhs_==NULL&&rowrange_==NULL);
01139 int numrows = problem_->num_constraints();
01140 if(numrows == 0) return;
01141 const double * rowLower = getRowLower();
01142 const double * rowUpper = getRowUpper();
01143 rowsense_ = new char [numrows];
01144 rhs_ = new double [numrows];
01145 rowrange_ = new double [numrows];
01146 for(int i = 0 ; i < numrows ; i++) {
01147 rowrange_[i]=0.;
01148 convertBoundToSense(rowLower[i], rowUpper[i], rowsense_[i], rhs_[i], rowrange_[i]);
01149 }
01150 }
01151
01161 const char *
01162 OsiTMINLPInterface::getRowSense() const
01163 {
01164 if(rowsense_==NULL) {
01165 extractSenseRhsAndRange();
01166 }
01167 return rowsense_;
01168 }
01169
01178 const double *
01179 OsiTMINLPInterface::getRightHandSide() const
01180 {
01181 if(rhs_==NULL) {
01182 extractSenseRhsAndRange();
01183 }
01184 return rhs_;
01185 }
01186
01195 const double *
01196 OsiTMINLPInterface::getRowRange() const
01197 {
01198 if(rowrange_==NULL) {
01199 extractSenseRhsAndRange();
01200 }
01201 return rowrange_;
01202 }
01203
01204 const double *
01205 OsiTMINLPInterface::getRowLower() const
01206 {
01207 return problem_->g_l();
01208 }
01209
01210 const double *
01211 OsiTMINLPInterface::getRowUpper() const
01212 {
01213 return problem_->g_u();
01214 }
01215
01217 bool
01218 OsiTMINLPInterface::isContinuous(int colNumber) const
01219 {
01220 return (problem_->var_types()[colNumber]==TMINLP::CONTINUOUS);
01221 }
01222
01224 bool
01225 OsiTMINLPInterface::isBinary(int colNumber) const
01226 {
01227 return (problem_->var_types()[colNumber]==TMINLP::BINARY);
01228 }
01229
01234 bool
01235 OsiTMINLPInterface::isInteger(int colNumber) const
01236 {
01237 return ((problem_->var_types()[colNumber]==TMINLP::BINARY)||
01238 (problem_->var_types()[colNumber]==TMINLP::INTEGER));
01239 }
01240
01242 bool
01243 OsiTMINLPInterface::isIntegerNonBinary(int colNumber) const
01244 {
01245 return (problem_->var_types()[colNumber]==TMINLP::INTEGER);
01246 }
01248 bool
01249 OsiTMINLPInterface::isFreeBinary(int colNumber) const
01250 {
01251 return ((problem_->var_types()[colNumber]==TMINLP::BINARY)
01252 &&((getColUpper()[colNumber]-getColLower()[colNumber]) > 1 - 1e-09));
01253 }
01254
01256 double
01257 OsiTMINLPInterface::getInfinity() const
01258 {
01259 return COIN_DBL_MAX;
01260 }
01261
01263 const double *
01264 OsiTMINLPInterface::getColSolution() const
01265 {
01266 if(hasBeenOptimized_)
01267 return problem_->x_sol();
01268 else
01269 return problem_->x_init();
01270 }
01271
01273 const double *
01274 OsiTMINLPInterface::getRowPrice() const
01275 {
01276 if(hasBeenOptimized_)
01277 return problem_->duals_sol();
01278 else
01279 return problem_->duals_init();
01280 }
01281
01283 const double *
01284 OsiTMINLPInterface::getReducedCost() const
01285 {
01286 (*handler_)<<"WARNING : trying to access reduced cost in Ipopt always retrun 0"<<CoinMessageEol;
01287 if(reducedCosts_==NULL) {
01288 reducedCosts_ = new double [getNumCols()];
01289 CoinFillN(reducedCosts_,getNumCols(),0.);
01290 }
01291 return reducedCosts_;
01292 }
01293
01296 const double *
01297 OsiTMINLPInterface::getRowActivity() const
01298 {
01299 return problem_->g_sol();
01300 }
01301
01305 int
01306 OsiTMINLPInterface::getIterationCount() const
01307 {
01308 return app_->IterationCount();
01309 }
01310
01311
01314 void
01315 OsiTMINLPInterface::setColLower( int elementIndex, double elementValue )
01316 {
01317
01318 problem_->SetVariableLowerBound(elementIndex,elementValue);
01319 hasBeenOptimized_ = false;
01320 }
01321
01324 void
01325 OsiTMINLPInterface::setColUpper( int elementIndex, double elementValue )
01326 {
01327
01328 problem_->SetVariableUpperBound(elementIndex,elementValue);
01329 hasBeenOptimized_ = false;
01330 }
01331
01334 void
01335 OsiTMINLPInterface::setColLower( const double* array )
01336 {
01337 problem_->SetVariablesLowerBounds(problem_->num_variables(),
01338 array);
01339 hasBeenOptimized_ = false;
01340 }
01341
01344 void
01345 OsiTMINLPInterface::setColUpper( const double* array )
01346 {
01347 problem_->SetVariablesUpperBounds(problem_->num_variables(),
01348 array);
01349 hasBeenOptimized_ = false;
01350 }
01351
01354 void
01355 OsiTMINLPInterface::setRowLower( int elementIndex, double elementValue )
01356 {
01357 throw SimpleError("Not implemented yet but should be if necessary.",
01358 "setRowLower");
01359 hasBeenOptimized_ = false;
01360 }
01361
01364 void
01365 OsiTMINLPInterface::setRowUpper( int elementIndex, double elementValue )
01366 {
01367 throw SimpleError("Not implemented yet but should be if necessary.",
01368 "setRowUpper");
01369 hasBeenOptimized_ = false;
01370 }
01371
01373 void
01374 OsiTMINLPInterface::setRowType(int index, char sense, double rightHandSide,
01375 double range)
01376 {
01377 throw SimpleError("Not implemented yet but should be if necessary.",
01378 "setRowType");
01379 hasBeenOptimized_ = false;
01380 }
01381
01382
01385 void
01386 OsiTMINLPInterface::setObjSense(double s)
01387 {
01388 throw SimpleError("Can not change objective sense of an Ipopt problem.",
01389 "setObjSense");
01390 hasBeenOptimized_ = false;
01391 }
01392
01402 void
01403 OsiTMINLPInterface::setColSolution(const double *colsol)
01404 {
01405 if(colsol != NULL)
01406 problem_->setxInit(getNumCols(), colsol);
01407 else
01408 problem_->resetStartingPoint();
01409 hasBeenOptimized_ = false;
01410 }
01411
01422 void
01423 OsiTMINLPInterface::setRowPrice(const double * rowprice)
01424 {
01425 problem_->setDualsInit(getNumCols()*2 + getNumRows(), rowprice);
01426 hasBeenOptimized_ = false;
01427 }
01428
01435 CoinWarmStart *
01436 OsiTMINLPInterface::getEmptyWarmStart () const
01437 {return app_->getEmptyWarmStart();}
01438
01440 CoinWarmStart*
01441 OsiTMINLPInterface::getWarmStart() const
01442 {
01443 if (warmStartMode_ >= Optimum) {
01444 return internal_getWarmStart();;
01445 }
01446 else if(warmStartMode_ == FakeBasis){
01447 return build_fake_basis();
01448 }
01449 else {
01450 return getEmptyWarmStart();
01451 }
01452 }
01455 bool
01456 OsiTMINLPInterface::setWarmStart(const CoinWarmStart* ws)
01457 {
01458 if (warmStartMode_ >= Optimum) {
01459 return internal_setWarmStart(ws);
01460 }
01461 else {
01462 return true;
01463 }
01464 }
01465
01466
01467 #define EPSILON 1e-4
01468
01469 CoinWarmStart*
01470 OsiTMINLPInterface::build_fake_basis() const{
01471 CoinWarmStartBasis * r_val = new CoinWarmStartBasis();
01472 int n_con = problem_->num_constraints();
01473 r_val->setSize(problem_->num_variables(), n_con);
01474 const double * act = problem_->g_sol();
01475 const double * lb = problem_->g_l();
01476 const double * ub = problem_->g_u();
01477 for(int i = 0 ; i < n_con ; i++){
01478 if(lb[i] > ub[i] - EPSILON){
01479 r_val->setArtifStatus(i, CoinWarmStartBasis::isFree);
01480 }
01481 if(act[i] > ub[i] - EPSILON){
01482 r_val->setArtifStatus(i, CoinWarmStartBasis::atLowerBound);
01483 }
01484 else if(act[i] < lb[i] + EPSILON){
01485 r_val->setArtifStatus(i, CoinWarmStartBasis::atLowerBound);
01486 }
01487 else {
01488 r_val->setArtifStatus(i, CoinWarmStartBasis::basic);
01489 }
01490 }
01491 return r_val;
01492 }
01493
01495 CoinWarmStart*
01496 OsiTMINLPInterface::internal_getWarmStart() const
01497 {
01498 if (warmStartMode_ >= Optimum && warmstart_) {
01499 return warmstart_->clone();
01500 }
01501 else {
01502 return getEmptyWarmStart();
01503 }
01504 }
01507 bool
01508 OsiTMINLPInterface::internal_setWarmStart(const CoinWarmStart* ws)
01509 {
01510 delete warmstart_;
01511 warmstart_ = NULL;
01512 hasBeenOptimized_ = false;
01513 if (warmStartMode_ >= Optimum) {
01514 if (ws == NULL) {
01515 return true;
01516 }
01517 if(app_->warmStartIsValid(ws)) {
01518 warmstart_ = ws->clone();
01519 return true;
01520 }
01521
01522
01523 const CoinWarmStartPrimalDual* pdws =
01524 dynamic_cast<const CoinWarmStartPrimalDual*>(ws);
01525 if (pdws) {
01526
01527
01528 warmstart_ = new IpoptWarmStart(*pdws);
01529 return true;
01530 }
01531 return false;
01532 }
01533 else {
01534 return true;
01535 }
01536 }
01537
01539 void
01540 OsiTMINLPInterface::setContinuous(int index)
01541 {
01542 problem_->SetVariableType(index, TMINLP::CONTINUOUS);
01543 hasBeenOptimized_ = false;
01544 }
01546 void
01547 OsiTMINLPInterface::setInteger(int index)
01548 {
01549 problem_->SetVariableType(index, TMINLP::INTEGER);
01550 hasBeenOptimized_ = false;
01551 }
01552
01554 double
01555 OsiTMINLPInterface::getObjValue() const
01556 {
01557 return problem_->obj_value();
01558 }
01559
01560
01561
01562
01563
01564 bool
01565 OsiTMINLPInterface::setIntParam(OsiIntParam key, int value)
01566 {
01567
01568
01569 bool retval = false;
01570 switch (key) {
01571 case OsiMaxNumIteration:
01572 retval = false;
01573 break;
01574 case OsiMaxNumIterationHotStart:
01575 if( value >= 0 ) {
01576 retval = false;
01577 }
01578 else
01579 retval = false;
01580 break;
01581 case OsiLastIntParam:
01582 retval = false;
01583 break;
01584 default:
01585 retval = false;
01586 (*handler_)<< "Unhandled case in setIntParam\n"<<CoinMessageEol;
01587 break;
01588 }
01589 return retval;
01590 }
01591
01592
01593
01594 bool
01595 OsiTMINLPInterface::setDblParam(OsiDblParam key, double value)
01596 {
01597
01598
01599 bool retval = false;
01600 switch (key) {
01601 case OsiDualObjectiveLimit:
01602 OsiDualObjectiveLimit_ = value;
01603 retval = true;
01604 break;
01605 case OsiPrimalObjectiveLimit:
01606 (*handler_)<<"Can not set primal objective limit parameter"<<CoinMessageEol;
01607 retval = false;
01608 break;
01609 case OsiDualTolerance:
01610 (*handler_)<<"Can not set dual tolerance parameter"<<CoinMessageEol;
01611 retval = false;
01612 break;
01613 case OsiPrimalTolerance:
01614 (*handler_)<<"Can not set primal tolerance parameter"<<CoinMessageEol;
01615 retval = false;
01616 case OsiObjOffset:
01617 retval = OsiSolverInterface::setDblParam(key,value);
01618 break;
01619 case OsiLastDblParam:
01620 retval = false;
01621 break;
01622 default:
01623 retval = false;
01624 (*handler_) << "Unhandled case in setDblParam"<<CoinMessageEol;
01625 break;
01626 }
01627 return retval;
01628 }
01629
01630
01631
01632
01633 bool
01634 OsiTMINLPInterface::setStrParam(OsiStrParam key, const std::string & value)
01635 {
01636
01637
01638 bool retval=false;
01639 switch (key) {
01640 case OsiProbName:
01641 OsiSolverInterface::setStrParam(key,value);
01642 return retval = true;
01643 case OsiSolverName:
01644 return false;
01645 case OsiLastStrParam:
01646 return false;
01647 }
01648 return false;
01649 }
01650
01651
01652
01653 bool
01654 OsiTMINLPInterface::getIntParam(OsiIntParam key, int& value) const
01655 {
01656
01657
01658 value = -COIN_INT_MAX;
01659 bool retval = false;
01660 switch (key) {
01661 case OsiMaxNumIteration:
01662 retval = false;
01663 break;
01664 case OsiMaxNumIterationHotStart:
01665 retval = false;
01666 break;
01667 case OsiLastIntParam:
01668 retval = false;
01669 break;
01670 default:
01671 retval = false;
01672 (*handler_) << "Unhandled case in setIntParam"<<CoinMessageEol;
01673 }
01674 return retval;
01675 }
01676
01677
01678
01679 bool
01680 OsiTMINLPInterface::getDblParam(OsiDblParam key, double& value) const
01681 {
01682
01683
01684 bool retval = false;
01685 switch (key) {
01686 case OsiDualObjectiveLimit:
01687 value = OsiDualObjectiveLimit_;
01688 retval = true;
01689 break;
01690 case OsiPrimalObjectiveLimit:
01691 value = getInfinity();
01692 retval = true;
01693 break;
01694 case OsiDualTolerance:
01695 retval = false;
01696 break;
01697 case OsiPrimalTolerance:
01698 options()->GetNumericValue("tol", value,"");
01699 value = 1e-07;
01700 retval = true;
01701 break;
01702 case OsiObjOffset:
01703 retval = OsiSolverInterface::getDblParam(key, value);
01704 break;
01705 case OsiLastDblParam:
01706 retval = false;
01707 break;
01708 }
01709 return retval;
01710 }
01711
01712
01713
01714
01715 bool
01716 OsiTMINLPInterface::getStrParam(OsiStrParam key, std::string & value) const
01717 {
01718
01719
01720 switch (key) {
01721 case OsiProbName:
01722 OsiSolverInterface::getStrParam(key, value);
01723 break;
01724 case OsiSolverName:
01725 value = "Ipopt";
01726 break;
01727 case OsiLastStrParam:
01728 return false;
01729 }
01730
01731 return true;
01732 }
01733
01734 void
01735 OsiTMINLPInterface::set_linearizer(Ipopt::SmartPtr<TMINLP2OsiLP> linearizer)
01736 {
01737 linearizer_ = linearizer->clone();
01738 linearizer_->set_tols(tiny_, veryTiny_, rhsRelax_, infty_);
01739 linearizer_->set_model(GetRawPtr(problem_));
01740 }
01741
01742 Ipopt::SmartPtr<TMINLP2OsiLP>
01743 OsiTMINLPInterface::linearizer(){
01744 return linearizer_;
01745 }
01746
01747
01748 void
01749 OsiTMINLPInterface::randomStartingPoint()
01750 {
01751 int numcols = getNumCols();
01752 const double * colLower = getColLower();
01753 const double * colUpper = getColUpper();
01754 double * sol = new double[numcols];
01755 const Number * x_init = problem_->x_init_user();
01756 const double* perturb_radius = NULL;
01757 if (randomGenerationType_ == perturb_suffix) {
01758 const TMINLP::PerturbInfo* pertubinfo = tminlp_->perturbInfo();
01759 if (pertubinfo) {
01760 perturb_radius = pertubinfo->GetPerturbationArray();
01761 }
01762 if (!perturb_radius) {
01763 throw SimpleError("Can't use perturb_radius if no radii are given.",
01764 "randomStartingPoint");
01765 }
01766 }
01767 for(int i = 0 ; i < numcols ; i++) {
01768 if(randomGenerationType_ == uniform || x_init[i] < colLower[i] || x_init[i] > colUpper[i]) {
01769 double lower = std::min(-maxRandomRadius_,colUpper[i] - maxRandomRadius_);
01770 lower = std::max(colLower[i], lower);
01771 double upper = std::max(maxRandomRadius_,colLower[i] + maxRandomRadius_);
01772 upper = std::min(colUpper[i],upper);
01773 lower = std::min(upper,lower);
01774 upper = std::max(upper, lower);
01775 double interval = upper - lower;
01776 sol[i] = CoinDrand48()*(interval) + lower;
01777 }
01778 else if (randomGenerationType_ == perturb) {
01779 const double lower = std::max(x_init[i] - max_perturbation_, colLower[i]);
01780 const double upper = std::min(x_init[i] + max_perturbation_, colUpper[i]);
01781 const double interval = upper - lower;
01782 sol[i] = lower + CoinDrand48()*(interval);
01783 }
01784 else if (randomGenerationType_ == perturb_suffix) {
01785 const double radius = perturb_radius[i];
01786 const double lower = std::max(x_init[i] - radius*max_perturbation_, colLower[i]);
01787 const double upper = std::min(x_init[i] + radius*max_perturbation_, colUpper[i]);
01788 const double interval = upper - lower;
01789 sol[i] = lower + CoinDrand48()*(interval);
01790 }
01791 }
01792 app_->disableWarmStart();
01793 setColSolution(sol);
01794 delete [] sol;
01795 }
01796
01797
01798
01800 int OsiTMINLPInterface::initializeJacobianArrays()
01801 {
01802 Index n, m, nnz_h_lag;
01803 TNLP::IndexStyleEnum index_style;
01804 problem_to_optimize_->get_nlp_info( n, m, nnz_jac, nnz_h_lag, index_style);
01805
01806 if(jRow_ != NULL) delete [] jRow_;
01807 if(jCol_ != NULL) delete [] jCol_;
01808 if(jValues_ != NULL) delete [] jValues_;
01809
01810 jRow_ = new Index[nnz_jac];
01811 jCol_ = new Index[nnz_jac];
01812 jValues_ = new Number[nnz_jac];
01813 problem_to_optimize_->eval_jac_g(n, NULL, 0, m, nnz_jac, jRow_, jCol_, NULL);
01814 if(index_style == Ipopt::TNLP::FORTRAN_STYLE)
01815 {
01816 for(int i = 0 ; i < nnz_jac ; i++){
01817 jRow_[i]--;
01818 jCol_[i]--;
01819 }
01820 }
01821
01822 if(constTypes_ != NULL) delete [] constTypes_;
01823
01824
01825 constTypes_ = new TNLP::LinearityType[getNumRows()];
01826 problem_to_optimize_->get_constraints_linearity(getNumRows(), constTypes_);
01827
01828 for(int i = 0; i < getNumRows() ; i++) {
01829 if(constTypes_[i]==TNLP::NON_LINEAR) {
01830
01831 nNonLinear_++;
01832 }
01833 }
01834 return nnz_jac;
01835 }
01836
01837
01838 double
01839 OsiTMINLPInterface::getConstraintsViolation(const double *x, double &obj)
01840 {
01841 int numcols = getNumCols();
01842 int numrows = getNumRows();
01843 double * g = new double[numrows];
01844 tminlp_->eval_g(numcols, x, 1, numrows, g);
01845 const double * rowLower = getRowLower();
01846 const double * rowUpper = getRowUpper();
01847
01848
01849 double norm = 0;
01850 for(int i = 0; i< numrows ; i++) {
01851 if(!constTypes_ || constTypes_[i] == TNLP::NON_LINEAR) {
01852 double rowViolation = 0;
01853 if(rowLower[i] > -1e10)
01854 rowViolation = std::max(0.,rowLower[i] - g[i]);
01855
01856 if(rowUpper[i] < 1e10);
01857 rowViolation = std::max(rowViolation, g[i] - rowUpper[i]);
01858
01859 norm = rowViolation > norm ? rowViolation : norm;
01860 }
01861 }
01862 tminlp_->eval_f(numcols, x, 1, obj);
01863 delete [] g;
01864 return norm;
01865 }
01866
01868 double
01869 OsiTMINLPInterface::getNonLinearitiesViolation(const double *x, const double obj)
01870 {
01871 double f;
01872 double norm = getConstraintsViolation(x, f);
01873 assert((f - obj) > -1e-08);
01874 norm = (f - obj) > norm ? f - obj : norm;
01875 return norm;
01876 }
01877
01878
01879
01880
01881 static inline
01882 bool cleanNnz(double &value, double colLower, double colUpper,
01883 double rowLower, double rowUpper, double colsol,
01884 double & lb, double &ub, double tiny, double veryTiny,
01885 double infty)
01886 {
01887
01888 if(fabs(value)>= tiny) return 1;
01889
01890 if(fabs(value)<veryTiny) return 0;
01891
01892
01893 bool colUpBounded = colUpper < infty;
01894 bool colLoBounded = colLower > -infty;
01895 bool rowNotLoBounded = rowLower <= - infty;
01896 bool rowNotUpBounded = rowUpper >= infty;
01897 bool pos = value > 0;
01898 if(colUpBounded && pos && rowNotUpBounded) {
01899 lb += value * (colsol - colUpper);
01900 return 0;
01901 }
01902 else
01903 if(colUpBounded && !pos && rowNotLoBounded) {
01904 ub += value * (colsol - colUpper);
01905 return 0;
01906 }
01907 else
01908 if(colLoBounded && !pos && rowNotUpBounded) {
01909 lb += value * (colsol - colLower);
01910 return 0;
01911 }
01912 else
01913 if(colLoBounded && pos && rowNotLoBounded) {
01914 ub += value * (colsol - colLower);
01915 return 0;
01916 }
01917
01918 return 1;
01919 }
01920
01923 void
01924 OsiTMINLPInterface::getOuterApproximation(OsiCuts &cs, const double * x,
01925 int getObj, const double * x2,
01926 double theta, bool global)
01927 {
01928 if(IsValid(linearizer_) && x2 == NULL){
01929 linearizer_->get_oas(cs, x, getObj, global);
01930 return;
01931 }
01932 int n,m, nnz_jac_g, nnz_h_lag;
01933 TNLP::IndexStyleEnum index_style;
01934 problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
01935 if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
01936 initializeJacobianArrays();
01937 assert(jRow_ != NULL);
01938 assert(jCol_ != NULL);
01939 vector<double> g(m);
01940 problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
01941 problem_to_optimize_->eval_g(n,x,1,m,g());
01942
01943 vector<CoinPackedVector> cuts(nNonLinear_ + 1);
01944 vector<double> lb(nNonLinear_ + 1);
01945 vector<double> ub(nNonLinear_ + 1);
01946
01947 vector<int> row2cutIdx(m,-1);
01948 int numCuts = 0;
01949
01950 const double * rowLower = getRowLower();
01951 const double * rowUpper = getRowUpper();
01952 const double * colLower = getColLower();
01953 const double * colUpper = getColUpper();
01954 const double * duals = getRowPrice() + 2 * n;
01955 double infty = getInfinity();
01956
01957 for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
01958 if(constTypes_[rowIdx] == TNLP::NON_LINEAR) {
01959 row2cutIdx[rowIdx] = numCuts;
01960 if(rowLower[rowIdx] > - infty_)
01961 lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
01962 else
01963 lb[numCuts] = - infty;
01964 if(rowUpper[rowIdx] < infty_)
01965 ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
01966 else
01967 ub[numCuts] = infty;
01968 if(rowLower[rowIdx] > -infty && rowUpper[rowIdx] < infty)
01969 {
01970 if(duals[rowIdx] >= 0)
01971 lb[numCuts] = - infty;
01972 if(duals[rowIdx] <= 0)
01973 ub[numCuts] = infty;
01974 }
01975
01976 numCuts++;
01977 }
01978 }
01979
01980
01981 for(int i = 0 ; i < nnz_jac_g ; i++) {
01982 const int &rowIdx = jRow_[i];
01983 const int & cutIdx = row2cutIdx[ rowIdx ];
01984 if(cutIdx != -1) {
01985 const int &colIdx = jCol_[i];
01986
01987 if(cleanNnz(jValues_[i],colLower[colIdx], colUpper[colIdx],
01988 rowLower[rowIdx], rowUpper[rowIdx],
01989 x[colIdx],
01990 lb[cutIdx],
01991 ub[cutIdx], tiny_, veryTiny_, infty_)) {
01992 cuts[cutIdx].insert(colIdx,jValues_[i]);
01993 if(lb[cutIdx] > - infty)
01994 lb[cutIdx] += jValues_[i] * x[colIdx];
01995 if(ub[cutIdx] < infty)
01996 ub[cutIdx] += jValues_[i] * x[colIdx];
01997 }
01998 }
01999 }
02000
02001 vector<int> cut2rowIdx(0);
02002 if (IsValid(cutStrengthener_) || oaHandler_->logLevel() > 0) {
02003 cut2rowIdx.resize(numCuts);
02004 for(int rowIdx = 0 ; rowIdx < m ; rowIdx++){
02005 if(row2cutIdx[rowIdx] >= 0){
02006 cut2rowIdx[row2cutIdx[rowIdx]] = rowIdx;
02007 }
02008 }
02009 }
02010
02011 for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
02012
02013 if(x2 != NULL) {
02014 double rhs = cuts[cutIdx].dotProduct(x2);
02015 double violation = 0.;
02016 violation = std::max(violation, rhs - ub[cutIdx]);
02017 violation = std::max(violation, lb[cutIdx] - rhs);
02018 if(violation < theta && oaHandler_->logLevel() > 0) {
02019 oaHandler_->message(CUT_NOT_VIOLATED_ENOUGH, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
02020 continue;}
02021 if(oaHandler_->logLevel() > 0)
02022 oaHandler_->message(VIOLATED_OA_CUT_GENERATED, oaMessages_)<<cut2rowIdx[cutIdx]<<violation<<CoinMessageEol;
02023 }
02024 OsiRowCut newCut;
02025
02026
02027 if (IsValid(cutStrengthener_)) {
02028 const int& rowIdx = cut2rowIdx[cutIdx];
02029 bool retval =
02030 cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
02031 GetRawPtr(problem_), rowIdx,
02032 cuts[cutIdx], lb[cutIdx], ub[cutIdx], g[rowIdx],
02033 rowLower[rowIdx], rowUpper[rowIdx],
02034 n, x, infty);
02035 if (!retval) {
02036 (*messageHandler()) << "error in cutStrengthener_->ComputeCuts\n";
02037
02038 }
02039 }
02040 if(global) {
02041 newCut.setGloballyValidAsInteger(1);
02042 }
02043 if(lb[cutIdx] > infty) lb[cutIdx] -= rhsRelax_*std::max(fabs(lb[cutIdx]), 1.);
02044 if(ub[cutIdx] < infty) ub[cutIdx] += rhsRelax_*std::max(fabs(ub[cutIdx]), 1.);
02045 newCut.setLb(lb[cutIdx]);
02046 newCut.setUb(ub[cutIdx]);
02047 newCut.setRow(cuts[cutIdx]);
02048 if(oaHandler_->logLevel()>2){
02049 oaHandler_->print(newCut);}
02050 cs.insert(newCut);
02051 }
02052
02053 if(getObj == 2 || (getObj && !problem_->hasLinearObjective())) {
02054 vector<double> obj(n);
02055 problem_to_optimize_->eval_grad_f(n, x, 1,obj());
02056 double f;
02057 problem_to_optimize_->eval_f(n, x, 1, f);
02058
02059 CoinPackedVector v;
02060 v.reserve(n);
02061 lb[nNonLinear_] = -f;
02062 ub[nNonLinear_] = -f;
02063
02064 for(int i = 0; i<n ; i++) {
02065 if(cleanNnz(obj[i],colLower[i], colUpper[i],
02066 -getInfinity(), 0,
02067 x[i],
02068 lb[nNonLinear_],
02069 ub[nNonLinear_],tiny_, 1e-15, infty_)) {
02070 v.insert(i,obj[i]);
02071 lb[nNonLinear_] += obj[i] * x[i];
02072 ub[nNonLinear_] += obj[i] * x[i];
02073 }
02074 }
02075 v.insert(n,-1);
02076
02077 bool genCut = true;
02078 if(x2 != NULL) {
02079 double rhs = v.dotProduct(x2);
02080 double violation = std::max(0., rhs - ub[nNonLinear_]);
02081 if(violation < theta) genCut = false;
02082 }
02083 if(genCut) {
02084 if (IsValid(cutStrengthener_)) {
02085 lb[nNonLinear_] = -infty;
02086 bool retval =
02087 cutStrengthener_->ComputeCuts(cs, GetRawPtr(tminlp_),
02088 GetRawPtr(problem_), -1,
02089 v, lb[nNonLinear_], ub[nNonLinear_],
02090 ub[nNonLinear_], -infty, 0.,
02091 n, x, infty);
02092 if (!retval) {
02093 (*handler_)<< "error in cutStrengthener_->ComputeCuts"<<CoinMessageEol;
02094
02095 }
02096 }
02097 OsiRowCut newCut;
02098 if(global)
02099 newCut.setGloballyValidAsInteger(1);
02100
02101 newCut.setRow(v);
02102 newCut.setLb(-COIN_DBL_MAX);
02103 newCut.setUb(ub[nNonLinear_]);
02104 cs.insert(newCut);
02105 }
02106 }
02107 }
02108
02110 void
02111 OsiTMINLPInterface::getBendersCut(OsiCuts &cs,
02112 bool global){
02113 int n,m, nnz_jac_g, nnz_h_lag;
02114 TNLP::IndexStyleEnum index_style;
02115 problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
02116 if(jRow_ == NULL || jCol_ == NULL || jValues_ == NULL)
02117 initializeJacobianArrays();
02118 assert(jRow_ != NULL);
02119 assert(jCol_ != NULL);
02120 vector<double> g(m);
02121 const double * x = getColSolution();
02122 problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
02123 problem_to_optimize_->eval_g(n,x,1,m,g());
02124
02125 vector<double> cut(n+1,0.);
02126 vector<bool> keep(m+1,false);
02127 double lb = 0;
02128 double ub = 0;
02129
02130 const double * rowLower = getRowLower();
02131 const double * rowUpper = getRowUpper();
02132 const double * colLower = getColLower();
02133 const double * colUpper = getColUpper();
02134 const double * duals = getRowPrice() + 2 * n;
02135
02136
02137 for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
02138 if(constTypes_[rowIdx] == TNLP::NON_LINEAR && fabs(duals[rowIdx]) > 1e-06)
02139 {
02140 const double & lam = duals[rowIdx];
02141 keep[rowIdx] = true;
02142 assert(lam < 0 || rowUpper[rowIdx] < 1e10);
02143 assert(lam > 0 || rowLower[rowIdx] > -1e10);
02144 if(lam < 0){
02145 assert(rowLower[rowIdx] > -1e10);
02146 ub += lam*(rowLower[rowIdx] -g[rowIdx]);
02147 }
02148 else {
02149 assert(rowUpper[rowIdx] < 1e10);
02150 ub += lam*(rowUpper[rowIdx] -g[rowIdx]);
02151 }
02152 }
02153 }
02154
02155
02156 for(int i = 0 ; i < nnz_jac_g ; i++) {
02157 const int &rowIdx = jRow_[i];
02158 if (!keep[rowIdx]) continue;
02159 const int &colIdx = jCol_[i];
02160
02161 const double & lam = duals[rowIdx];
02162 double coeff = lam*jValues_[i];
02163 if(cleanNnz(coeff,colLower[colIdx], colUpper[colIdx],
02164 rowLower[rowIdx], rowUpper[rowIdx], x[colIdx], lb,
02165 ub, tiny_, veryTiny_, infty_)) {
02166 cut[colIdx] += coeff;
02167 ub += coeff * x[colIdx];
02168 }
02169 }
02170
02171 CoinPackedVector v;
02172 if(!problem_->hasLinearObjective() &&
02173 isProvenOptimal()) {
02174 vector<double> obj(n);
02175 problem_to_optimize_->eval_grad_f(n, x, 1,obj());
02176 double f;
02177 problem_to_optimize_->eval_f(n, x, 1, f);
02178
02179 ub = -f;
02180
02181 for(int i = 0; i<n ; i++) {
02182 if(cleanNnz(obj[i],colLower[i], colUpper[i], -getInfinity(), 0,
02183 x[i], lb, ub,tiny_, 1e-15, infty_)) {
02184 cut[i] += obj[i];
02185 ub += obj[i] * x[i];
02186 }
02187 }
02188 v.insert(n,-1);
02189 }
02190 for(int i = 0 ; i < n ; i++){
02191 if(fabs(cut[i])>1e-020){
02192 v.insert(i, cut[i]);
02193 }
02194 }
02195 OsiRowCut newCut;
02196 if(global)
02197 newCut.setGloballyValidAsInteger(1);
02198 newCut.setLb(-COIN_DBL_MAX);
02199 newCut.setUb(ub);
02200 newCut.setRow(v);
02201 cs.insert(newCut);
02202 }
02203
02206 void
02207 OsiTMINLPInterface::getConstraintOuterApproximation(OsiCuts &cs, int rowIdx,
02208 const double * x,
02209 const double * x2, bool global)
02210 {
02211 double g;
02212 int * indices = new int[getNumCols()];
02213 double * values = new double[getNumCols()];
02214 int nnz;
02215 problem_->eval_grad_gi(getNumCols(), x, 1, rowIdx, nnz, indices, values);
02216 problem_->eval_gi(getNumCols(),x,1, rowIdx, g);
02217
02218 CoinPackedVector cut;
02219 double lb;
02220 double ub;
02221
02222
02223 const double rowLower = getRowLower()[rowIdx];
02224 const double rowUpper = getRowUpper()[rowIdx];
02225 const double * colLower = getColLower();
02226 const double * colUpper = getColUpper();
02227 const double dual = (getRowPrice() + 2 * getNumCols())[rowIdx];
02228 double infty = getInfinity();
02229
02230 if(rowLower > - infty_)
02231 lb = rowLower - g;
02232 else
02233 lb = - infty;
02234 if(rowUpper < infty_)
02235 ub = rowUpper - g;
02236 else
02237 ub = infty;
02238 if(rowLower > -infty && rowUpper < infty)
02239 {
02240 if(dual >= 0)
02241 lb = - infty;
02242 if(dual <= 0)
02243 ub = infty;
02244 }
02245
02246 for(int i = 0 ; i < nnz; i++) {
02247 const int &colIdx = indices[i];
02248
02249 if(cleanNnz(values[i],colLower[colIdx], colUpper[colIdx],
02250 rowLower, rowUpper,
02251 x[colIdx],
02252 lb,
02253 ub, tiny_, veryTiny_, infty_)) {
02254 cut.insert(colIdx,values[i]);
02255 if(lb > - infty)
02256 lb += values[i] * x[colIdx];
02257 if(ub < infty)
02258 ub += values[i] * x[colIdx];
02259 }
02260 }
02261
02262 OsiRowCut newCut;
02263
02264 if(global) {
02265 newCut.setGloballyValidAsInteger(1);
02266 }
02267
02268 newCut.setLb(lb);
02269 newCut.setUb(ub);
02270 newCut.setRow(cut);
02271 cs.insert(newCut);
02272
02273 delete [] indices;
02274 delete [] values;
02275 }
02276
02277 void
02278 OsiTMINLPInterface::switchToFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
02279 double a, double s, int L){
02280 if(! IsValid(feasibilityProblem_)) {
02281 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02282 }
02283 feasibilityProblem_->set_use_feasibility_pump_objective(true);
02284 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
02285 feasibilityProblem_->setLambda(a);
02286 feasibilityProblem_->setSigma(s);
02287 feasibilityProblem_->setNorm(L);
02288 feasibilityProblem_->set_use_cutoff_constraint(false);
02289 feasibilityProblem_->set_use_local_branching_constraint(false);
02290 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
02291 feasibility_mode_ = true;
02292 }
02293
02294 void
02295 OsiTMINLPInterface::switchToFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
02296 double rhs_local_branching_constraint){
02297 if(! IsValid(feasibilityProblem_)) {
02298 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02299 }
02300 feasibilityProblem_->set_use_feasibility_pump_objective(false);
02301 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
02302 feasibilityProblem_->set_use_cutoff_constraint(false);
02303 feasibilityProblem_->set_use_local_branching_constraint(true);
02304 feasibilityProblem_->set_rhs_local_branching_constraint(rhs_local_branching_constraint);
02305 problem_to_optimize_ = GetRawPtr(feasibilityProblem_);
02306 feasibility_mode_ = true;
02307 }
02308
02309 void
02310 OsiTMINLPInterface::switchToOriginalProblem(){
02311 problem_to_optimize_ = GetRawPtr(problem_);
02312 feasibility_mode_ = false;
02313 }
02314
02315 double
02316 OsiTMINLPInterface::solveFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
02317 double a, double s, int L)
02318 {
02319 if(! IsValid(feasibilityProblem_)) {
02320 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02321 }
02322 feasibilityProblem_->set_use_feasibility_pump_objective(true);
02323 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
02324 feasibilityProblem_->setLambda(a);
02325 feasibilityProblem_->setSigma(s);
02326 feasibilityProblem_->setNorm(L);
02327 feasibilityProblem_->set_use_cutoff_constraint(false);
02328 feasibilityProblem_->set_use_local_branching_constraint(false);
02329 nCallOptimizeTNLP_++;
02330 totalNlpSolveTime_-=CoinCpuTime();
02331 SmartPtr<TNLPSolver> app2 = app_->clone();
02332 app2->options()->SetIntegerValue("print_level", (Index) 0);
02333 optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
02334 totalNlpSolveTime_+=CoinCpuTime();
02335 hasBeenOptimized_=true;
02336 return getObjValue();
02337 }
02338
02339 double
02340 OsiTMINLPInterface::solveFeasibilityProblem(size_t n,const double * x_bar,const int *inds,
02341 int L, double cutoff)
02342 {
02343 if(! IsValid(feasibilityProblem_)) {
02344 throw SimpleError("No feasibility problem","getFeasibilityOuterApproximation");
02345 }
02346 feasibilityProblem_->set_use_feasibility_pump_objective(true);
02347 feasibilityProblem_->set_dist_to_point_obj(n,(const Number *) x_bar,(const Index *) inds);
02348 feasibilityProblem_->setLambda(1.0);
02349 feasibilityProblem_->setSigma(0.0);
02350 feasibilityProblem_->setNorm(L);
02351 feasibilityProblem_->set_use_cutoff_constraint(true);
02352 feasibilityProblem_->set_cutoff(cutoff);
02353 feasibilityProblem_->set_use_local_branching_constraint(false);
02354 nCallOptimizeTNLP_++;
02355 totalNlpSolveTime_-=CoinCpuTime();
02356 SmartPtr<TNLPSolver> app2 = app_->clone();
02357 app2->options()->SetIntegerValue("print_level", (Index) 0);
02358 optimizationStatus_ = app2->OptimizeTNLP(GetRawPtr(feasibilityProblem_));
02359 totalNlpSolveTime_+=CoinCpuTime();
02360 hasBeenOptimized_=true;
02361 return getObjValue();
02362 }
02363
02364 #if 0
02365 double
02366 OsiTMINLPInterface::getFeasibilityOuterApproximation(int n,const double * x_bar,const int *inds, OsiCuts &cs, bool addOnlyViolated, bool global)
02367 {
02368 double ret_val = solveFeasibilityProblem(n, x_bar, inds, 1, 0, 2);
02369 getOuterApproximation(cs, getColSolution(), 0, (addOnlyViolated? x_bar:NULL)
02370 , global);
02371 return ret_val;
02372 }
02373 #endif
02374
02375 static bool WarnedForNonConvexOa=false;
02376
02377 void
02378 OsiTMINLPInterface::extractLinearRelaxation(OsiSolverInterface &si,
02379 const double * x, bool getObj
02380 )
02381 {
02382 if(IsValid(linearizer_)){
02383 linearizer_->extract(&si, x, getObj);
02384 return;
02385 }
02386 int n;
02387 int m;
02388 int nnz_jac_g;
02389 int nnz_h_lag;
02390 TNLP::IndexStyleEnum index_style;
02391
02392 problem_to_optimize_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
02393
02394
02395
02396 initializeJacobianArrays();
02397
02398
02399 problem_to_optimize_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, jValues_);
02400
02401
02402 vector<double> g(m);
02403 problem_to_optimize_->eval_g(n, x, 1, m, g());
02404
02405 vector<double> rowLow(m);
02406 vector<double> rowUp(m);
02407 vector<int> nonBindings(m);;
02408 int numNonBindings = 0;
02409 const double * rowLower = getRowLower();
02410 const double * rowUpper = getRowUpper();
02411 vector<double> colLower(n);
02412 vector<double> colUpper(n);
02413 std::copy(getColLower(), getColLower() + n, colLower.begin());
02414 std::copy(getColUpper(), getColUpper() + n, colUpper.begin());
02415 const double * duals = getRowPrice() + 2*n;
02416 assert(m==getNumRows());
02417 double infty = si.getInfinity();
02418 for(int i = 0 ; i < m ; i++) {
02419 if(constTypes_[i] == TNLP::NON_LINEAR) {
02420
02421 if(rowLower[i] > -infty_ && rowUpper[i] < infty_ && fabs(duals[i]) == 0.)
02422 {
02423 nonBindings[numNonBindings++] = i;
02424 continue;
02425 }
02426 else
02427 if(rowLower[i] > - infty_){
02428 rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
02429 if(! WarnedForNonConvexOa && rowUpper[i] < infty_){
02430 messageHandler()->message(WARNING_NON_CONVEX_OA, messages_)<<CoinMessageEol;
02431 WarnedForNonConvexOa = true;
02432 }
02433 }
02434 else
02435 rowLow[i] = - infty;
02436 if(rowUpper[i] < infty_)
02437 rowUp[i] = (rowUpper[i] - g[i]) + 1e-07;
02438 else
02439 rowUp[i] = infty;
02440
02441
02442 if(rowLower[i] > -infty_ && rowUpper[i] < infty_)
02443 {
02444 if(duals[i] >= 0.)
02445 rowLow[i] = - infty;
02446 if(duals[i] <= 0.)
02447 rowUp[i] = infty;
02448 }
02449 }
02450 else {
02451 if(rowLower[i] > -infty_){
02452 rowLow[i] = (rowLower[i]);
02453 }
02454 else
02455 rowLow[i] = - infty;
02456 if(rowUpper[i] < infty_){
02457 rowUp[i] = (rowUpper[i]);
02458 }
02459 else
02460 rowUp[i] = infty;
02461 }
02462 }
02463
02464
02465
02466
02467
02468 for(int i = 0 ; i < nnz_jac_g ; i++) {
02469 if(constTypes_[jRow_[i]] != TNLP::LINEAR){
02470 if(
02471 cleanNnz(jValues_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
02472 rowLower[jRow_[i]], rowUpper[jRow_[i]],
02473 x[jCol_[i]],
02474 rowLow[jRow_[i]],
02475 rowUp[jRow_[i]], tiny_, veryTiny_, infty_)) {
02476 if(rowLow[jRow_[i]] > - infty)
02477 rowLow[jRow_[i]] += jValues_[i] * x[jCol_ [i]];
02478 if(rowUp[jRow_[i]] < infty)
02479 rowUp[jRow_[i]] += jValues_[i] *x[jCol_[i]];
02480 }
02481 }
02482 else {
02483
02484
02485
02486 }
02487 }
02488
02489
02490 CoinPackedMatrix mat(true, jRow_, jCol_, jValues_, nnz_jac_g);
02491 mat.setDimensions(m,n);
02492
02493
02494 mat.deleteRows(numNonBindings, nonBindings());
02495
02496 int numcols=getNumCols();
02497 vector<double> obj(numcols);
02498 for(int i = 0 ; i < numcols ; i++){
02499 obj[i] = 0.;
02500 if(colLower[i] <= - infty_) colLower[i] = - infty;
02501 if(colUpper[i] >= infty_) colUpper[i] = infty;
02502 }
02503
02504
02505
02506 for(size_t i = 0 ; i < rowLow.size() ; i++){
02507 if(rowLow[i] > infty) rowLow[i] -= rhsRelax_*std::max(fabs(rowLow[i]), 1.);
02508 if(rowUp[i] < infty) rowUp[i] += rhsRelax_*std::max(fabs(rowUp[i]), 1.);
02509 }
02510
02511 si.loadProblem(mat, colLower(), colUpper(), obj(), rowLow(), rowUp());
02512 for(int i = 0 ; i < numcols ; i++) {
02513 if(isInteger(i))
02514 si.setInteger(i);
02515 }
02516 if(getObj) {
02517 bool addObjVar = false;
02518 if(problem_->hasLinearObjective()){
02519 double zero;
02520 vector<double> x0(n,0.);
02521 problem_to_optimize_->eval_f(n, x0(), 1, zero);
02522 si.setDblParam(OsiObjOffset, -zero);
02523
02524 problem_to_optimize_->eval_grad_f(n, x, 1,obj());
02525 si.setObjective(obj());
02526 }
02527 else {
02528 addObjVar = true;
02529 }
02530
02531 if(addObjVar){
02532 addObjectiveFunction(si, x);
02533 }
02534 }
02535
02536
02537 }
02538
02539 void
02540 OsiTMINLPInterface::addObjectiveFunction(OsiSolverInterface &si,
02541 const double *x){
02542 const double * colLower = getColLower();
02543 const double * colUpper = getColUpper();
02544 int numcols = getNumCols();
02545 assert(numcols == si.getNumCols() );
02546 vector<double> obj(numcols);
02547 problem_to_optimize_->eval_grad_f(numcols, x, 1,obj());
02548
02549
02550 CoinPackedVector a;
02551 si.addCol(a,-si.getInfinity(), si.getInfinity(), 1.);
02552
02553
02554
02555 double ub;
02556 problem_to_optimize_->eval_f(numcols, x, 1, ub);
02557 ub*=-1;
02558 double lb = -1e300;
02559 CoinPackedVector objCut;
02560 CoinPackedVector * v = &objCut;
02561 v->reserve(numcols+1);
02562 for(int i = 0; i<numcols ; i++) {
02563 if(si.getNumRows())
02564 {
02565 if(cleanNnz(obj[i],colLower[i], colUpper[i],
02566 -getInfinity(), 0,
02567 x[i],
02568 lb,
02569 ub, tiny_, veryTiny_, infty_)) {
02570 v->insert(i,obj[i]);
02571 lb += obj[i] * x[i];
02572 ub += obj[i] * x[i];
02573 }
02574 }
02575 else
02576 {
02577 if(cleanNnz(obj[i],colLower[i], colUpper[i],
02578 -getInfinity(), 0,
02579 x[i],
02580 lb,
02581 ub, 1e-03, 1e-08, infty_)) {
02582 v->insert(i,obj[i]);
02583 lb += obj[i] * x[i];
02584 ub += obj[i] * x[i];
02585 }
02586 }
02587 }
02588 v->insert(numcols,-1);
02589 si.addRow(objCut, lb, ub);
02590 }
02591
02593 void
02594 OsiTMINLPInterface::applyRowCuts(int numberCuts, const OsiRowCut * cuts)
02595 {
02596 if(numberCuts)
02597 freeCachedRowRim();
02598 const OsiRowCut ** cutsPtrs = new const OsiRowCut*[numberCuts];
02599 for(int i = 0 ; i < numberCuts ; i++)
02600 {
02601 cutsPtrs[i] = &cuts[i];
02602 }
02603 problem_->addCuts(numberCuts, cutsPtrs);
02604 delete [] cutsPtrs;
02605 }
02606
02607 void
02608 OsiTMINLPInterface::solveAndCheckErrors(bool warmStarted, bool throwOnFailure,
02609 const char * whereFrom)
02610 {
02611 if (BonminAbortAll == true) return;
02612 totalNlpSolveTime_-=CoinCpuTime();
02613 if(warmStarted)
02614 optimizationStatus_ = app_->ReOptimizeTNLP(GetRawPtr(problem_to_optimize_));
02615 else
02616 optimizationStatus_ = app_->OptimizeTNLP(GetRawPtr(problem_to_optimize_));
02617 totalNlpSolveTime_+=CoinCpuTime();
02618 nCallOptimizeTNLP_++;
02619 hasBeenOptimized_ = true;
02620
02621 if(getRowCutDebugger()){
02622
02623 if(! (isProvenOptimal() && getObjValue() <= getRowCutDebugger()->optimalValue())){
02624 std::string probName;
02625 getStrParam(OsiProbName, probName);
02626 fprintf(stderr, "Problem on optimal path is infeasible!\n");
02627 throw newUnsolvedError(app_->errorCode(), problem_, probName);
02628 }
02629 }
02630
02631
02632 if(!hasPrintedOptions) {
02633 hasPrintedOptions = 1;
02634
02635 app_->options()->SetStringValue("print_user_options","no", false, true);
02636 }
02637
02638 bool otherDisagree = false ;
02639 if(!app_->isRecoverable(optimizationStatus_))
02640 {
02641 std::string probName;
02642 getStrParam(OsiProbName, probName);
02643 throw newUnsolvedError(app_->errorCode(), problem_, probName);
02644 }
02645 else if(testOthers_ && !app_->isError(optimizationStatus_)){
02646 Ipopt::SmartPtr<TMINLP2TNLP> problem_copy = problem_->clone();
02647
02648 int f =1;
02649 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
02650 i != debug_apps_.end() ; i++){
02651 TNLPSolver::ReturnStatus otherStatus = (*i)->OptimizeTNLP(GetRawPtr(problem_copy));
02652 messageHandler()->message(LOG_LINE, messages_)
02653 <<'d'<<f++<<statusAsString(otherStatus)<<problem_copy->obj_value()
02654 <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
02655 if(!(*i)->isError(otherStatus)){
02656 CoinRelFltEq eq(1e-05);
02657 if(otherStatus != optimizationStatus_){
02658 otherDisagree = true;
02659 messageHandler()->message(SOLVER_DISAGREE_STATUS, messages_)
02660 <<app_->solverName()<<statusAsString()
02661 <<(*i)->solverName()<<statusAsString(otherStatus)<<CoinMessageEol;
02662 }
02663 else if(isProvenOptimal() && !eq(problem_->obj_value(),problem_copy->obj_value()))
02664 {
02665 otherDisagree = true;
02666 messageHandler()->message(SOLVER_DISAGREE_VALUE, messages_)
02667 <<app_->solverName()<<problem_->obj_value()
02668 <<(*i)->solverName()<<problem_copy->obj_value()<<CoinMessageEol;
02669 }
02670 }
02671 }
02672 }
02673 else if(app_->isError(optimizationStatus_) && ! debug_apps_.empty()){
02674 int f =1;
02675 for(std::list<Ipopt::SmartPtr<TNLPSolver> >::iterator i = debug_apps_.begin();
02676 i != debug_apps_.end() && app_->isError(optimizationStatus_) ; i++){
02677 optimizationStatus_ = (*i)->OptimizeTNLP(GetRawPtr(problem_));
02678 messageHandler()->message(LOG_LINE, messages_)
02679 <<'d'<<f++<<statusAsString(optimizationStatus_)<<problem_->obj_value()
02680 <<(*i)->IterationCount()<<(*i)->CPUTime()<<"retry with "+(*i)->solverName()<<CoinMessageEol;
02681 }
02682 }
02683 try{
02684 totalIterations_ += app_->IterationCount();
02685 }
02686 catch(SimpleError &E)
02687 {
02688 if (throwOnFailure)
02689 {
02690 throw SimpleError("No statistics available from Ipopt",whereFrom);
02691 }
02692 else {
02693 return;
02694 }
02695 }
02696 if(problem_->hasUpperBoundingObjective()){
02697 const double * sol = getColSolution();
02698 bool integerSol = true;
02699 double intTol = 1e-08;
02700 if(objects()){
02701 int nObjects = numberObjects();
02702 OsiObject ** object = objects();
02703 for(int i = 0 ; i< nObjects ; i++){
02704 int dummy;
02705 if(object[i]->infeasibility(this,dummy) > intTol)
02706 {
02707 integerSol=false;
02708 break;
02709 }
02710 }
02711 }
02712 else{
02713 int numcols = getNumCols();
02714 for(int i = 0 ; i < numcols ; i++){
02715 if(isInteger(i) || isBinary(i)){
02716 if(fabs(sol[i] - floor(sol[i]+0.5)) > intTol){
02717 integerSol = false;
02718 break;
02719 }
02720 }
02721 }
02722 }
02723 if(integerSol&&isProvenOptimal()){
02724 double help= problem_->evaluateUpperBoundingFunction(sol);
02725
02726
02727 OsiAuxInfo * auxInfo = getAuxiliaryInfo();
02728 Bonmin::AuxInfo * bonInfo = dynamic_cast<Bonmin::AuxInfo *>(auxInfo);
02729 if(bonInfo!=0)
02730 {
02731
02732 if(help<bonInfo->bestObj2())
02733 {
02734 bonInfo->setBestObj2(help);
02735 bonInfo->setBestSolution2(getNumCols(), const_cast<double *>(getColSolution()));
02736
02737 messageHandler()->message(ALTERNATE_OBJECTIVE, messages_)
02738 <<help<<CoinMessageEol;
02739 }
02740 }
02741 else {
02742 printf("\nWARNING: the algorithm selected does not consider the second objective function\n");
02743 }
02744 }
02745 }
02746
02747 messageHandler()->message(IPOPT_SUMMARY, messages_)
02748 <<whereFrom<<optimizationStatus_<<app_->IterationCount()<<app_->CPUTime()<<CoinMessageEol;
02749
02750 if((nCallOptimizeTNLP_ % 20) == 1)
02751 messageHandler()->message(LOG_HEAD, messages_)<<CoinMessageEol;
02752
02753
02754 if ( (numIterationSuspect_ >= 0 && (getIterationCount()>numIterationSuspect_ || isAbandoned())) ||
02755 ( otherDisagree )){
02756 messageHandler()->message(SUSPECT_PROBLEM,
02757 messages_)<<nCallOptimizeTNLP_<<CoinMessageEol;
02758 std::string subProbName;
02759 getStrParam(OsiProbName, subProbName);
02760 std::ostringstream os;
02761 os<<"_"<<nCallOptimizeTNLP_;
02762 subProbName+=os.str();
02763 problem_->outputDiffs(subProbName, NULL);
02764 }
02765
02766 }
02767
02769
02771 void OsiTMINLPInterface::initialSolve()
02772 {
02773 initialSolve("");
02774 }
02775
02777
02779 void OsiTMINLPInterface::initialSolve(const char * whereFrom)
02780 {
02781 assert(IsValid(app_));
02782 assert(IsValid(problem_));
02783
02784 if (BonminAbortAll == true) return;
02785
02786 delete warmstart_;
02787 warmstart_ = NULL;
02788
02789 if(!hasPrintedOptions) {
02790 int printOptions;
02791 app_->options()->GetEnumValue("print_user_options",printOptions,app_->prefix());
02792 if(printOptions)
02793 app_->options()->SetStringValue("print_user_options","yes",true,true);
02794 }
02795 if(warmStartMode_ >= Optimum)
02796 app_->disableWarmStart();
02797 solveAndCheckErrors(0,1,"initialSolve");
02798
02799
02800 if(!hasPrintedOptions) {
02801 hasPrintedOptions = 1;
02802 app_->options()->SetStringValue("print_user_options","no");
02803 app_->options()->SetIntegerValue("print_level",0);
02804 }
02805
02806 messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
02807 <<statusAsString()
02808 <<getObjValue()
02809 <<app_->IterationCount()
02810 <<app_->CPUTime()
02811 <<whereFrom<<CoinMessageEol;
02812
02813 if(BonminAbortAll){
02814 return;
02815 }
02816 int numRetry = firstSolve_ ? numRetryInitial_ : numRetryResolve_;
02817 if(isAbandoned() ||
02818 ( isProvenPrimalInfeasible() && getObjValue() < infeasibility_epsilon_)) {
02819 resolveForRobustness(numRetryUnsolved_);
02820 }
02821 else if(numRetry)
02822 {
02823 resolveForCost(numRetry, numRetryInitial_ > 0);
02825 numRetryInitial_ = 0;
02826 }
02827 firstSolve_ = false;
02828
02829
02830
02831 if (! warmstart_ && ! isAbandoned()) {
02832 if (warmStartMode_ >= Optimum) {
02833 warmstart_ = app_->getWarmStart(problem_);
02834 }
02835 }
02836 }
02837
02839 void
02840 OsiTMINLPInterface::resolve(){
02841 resolve("");
02842 }
02845 void
02846 OsiTMINLPInterface::resolve(const char * whereFrom)
02847 {
02848 assert(IsValid(app_));
02849 assert(IsValid(problem_));
02850 if (BonminAbortAll == true) return;
02851 int has_warmstart = warmstart_ == NULL ? 0 : 1;
02852 if(warmstart_ == NULL) has_warmstart = 0;
02853 else if(!app_->warmStartIsValid(warmstart_)) has_warmstart = 1;
02854 else has_warmstart = 2;
02855 if (has_warmstart < 2) {
02856
02857 initialSolve(whereFrom);
02858 return;
02859 }
02860 app_->setWarmStart(warmstart_, problem_);
02861 delete warmstart_;
02862 warmstart_ = NULL;
02863
02864 if (INT_BIAS > 0.) {
02865 app_->options()->SetStringValue("warm_start_same_structure", "yes");
02866 }
02867 else {
02868 app_->options()->SetStringValue("warm_start_same_structure", "no");
02869 }
02870
02871 if(problem_->duals_init() != NULL)
02872 app_->enableWarmStart();
02873 else app_->disableWarmStart();
02874 solveAndCheckErrors(1,1,"resolve");
02875
02876 messageHandler()->message(LOG_LINE, messages_)<<' '<<nCallOptimizeTNLP_
02877 <<statusAsString()
02878 <<getObjValue()
02879 <<app_->IterationCount()
02880 <<app_->CPUTime()
02881 <<whereFrom
02882 <<"totot"
02883 <<CoinMessageEol;
02884
02885 if(isAbandoned() ||
02886 ( (getObjValue() < 1e-06) && isProvenPrimalInfeasible())) {
02887 resolveForRobustness(numRetryUnsolved_);
02888 }
02889 else if(numRetryResolve_ ||
02890 (numRetryInfeasibles_ && isProvenPrimalInfeasible() ))
02891 resolveForCost(std::max(numRetryResolve_, numRetryInfeasibles_), 0);
02892
02893
02894
02895 if (! warmstart_ && ! isAbandoned()) {
02896 if (warmStartMode_ >= Optimum) {
02897 warmstart_ = app_->getWarmStart(problem_);
02898 }
02899 }
02900 }
02901
02902
02904
02907 bool OsiTMINLPInterface::isAbandoned() const
02908 {
02909 if (pretendSucceededNext_)
02910 return false;
02911 pretendSucceededNext_ = false;
02912 return (
02913 (optimizationStatus_==TNLPSolver::iterationLimit)||
02914 (optimizationStatus_==TNLPSolver::computationError)||
02915 (optimizationStatus_==TNLPSolver::illDefinedProblem)||
02916 (optimizationStatus_==TNLPSolver::illegalOption)||
02917 (optimizationStatus_==TNLPSolver::externalException)|
02918 (optimizationStatus_==TNLPSolver::exception)
02919 );
02920 }
02921
02923 bool OsiTMINLPInterface::isProvenOptimal() const
02924 {
02925 return (optimizationStatus_==TNLPSolver::solvedOptimal) ||
02926 (optimizationStatus_==TNLPSolver::solvedOptimalTol);
02927 }
02929 bool OsiTMINLPInterface::isProvenPrimalInfeasible() const
02930 {
02931 return (optimizationStatus_ == TNLPSolver::provenInfeasible);
02932 }
02934 bool OsiTMINLPInterface::isProvenDualInfeasible() const
02935 {
02936 return (optimizationStatus_ == TNLPSolver::unbounded);
02937 }
02939 bool OsiTMINLPInterface::isPrimalObjectiveLimitReached() const
02940 {
02941 (*handler_)<<"Warning : isPrimalObjectiveLimitReached not implemented yet"<<CoinMessageEol;
02942 return 0;
02943 }
02945 bool OsiTMINLPInterface::isDualObjectiveLimitReached() const
02946 {
02947
02948 return (optimizationStatus_==TNLPSolver::unbounded);
02949
02950 }
02952 bool OsiTMINLPInterface::isIterationLimitReached() const
02953 {
02954 return (optimizationStatus_==TNLPSolver::iterationLimit);
02955 }
02956
02957 void
02958 OsiTMINLPInterface::extractInterfaceParams()
02959 {
02960 if (IsValid(app_)) {
02961 int logLevel;
02962 app_->options()->GetIntegerValue("nlp_log_level", logLevel,app_->prefix());
02963 messageHandler()->setLogLevel(logLevel);
02964
02965 #ifdef COIN_HAS_FILTERSQP
02966 FilterSolver * filter = dynamic_cast<FilterSolver *>(GetRawPtr(app_));
02967
02968 bool is_given =
02969 #endif
02970 app_->options()->GetNumericValue("max_random_point_radius",maxRandomRadius_,app_->prefix());
02971
02972 #ifdef COIN_HAS_FILTERSQP
02973 if(filter && !is_given){
02974
02975 maxRandomRadius_ = 10.;
02976 }
02977 #endif
02978
02979 int oaCgLogLevel = 0;
02980 app_->options()->GetIntegerValue("oa_cuts_log_level", oaCgLogLevel,app_->prefix());
02981 oaHandler_->setLogLevel(oaCgLogLevel);
02982
02983 int buffy;
02984 app_->options()->GetEnumValue("warm_start", buffy, app_->prefix());
02985 warmStartMode_ = (WarmStartModes) buffy;
02986
02987 app_->options()->GetIntegerValue("num_retry_unsolved_random_point", numRetryUnsolved_,app_->prefix());
02988 app_->options()->GetIntegerValue("num_resolve_at_root", numRetryInitial_,app_->prefix());
02989 app_->options()->GetIntegerValue("num_resolve_at_node", numRetryResolve_,app_->prefix());
02990 app_->options()->GetIntegerValue("num_resolve_at_infeasibles", numRetryInfeasibles_,app_->prefix());
02991 app_->options()->GetIntegerValue("num_iterations_suspect", numIterationSuspect_,app_->prefix());
02992 app_->options()->GetEnumValue("nlp_failure_behavior",pretendFailIsInfeasible_,app_->prefix());
02993 app_->options()->GetNumericValue
02994 ("warm_start_bound_frac" ,pushValue_,app_->prefix());
02995 app_->options()->GetNumericValue("tiny_element",tiny_,app_->prefix());
02996 app_->options()->GetNumericValue("very_tiny_element",veryTiny_,app_->prefix());
02997 app_->options()->GetNumericValue("oa_rhs_relax",rhsRelax_,app_->prefix());
02998 app_->options()->GetNumericValue("random_point_perturbation_interval",max_perturbation_,app_->prefix());
02999 app_->options()->GetEnumValue("random_point_type",randomGenerationType_,app_->prefix());
03000 int cut_strengthening_type;
03001 app_->options()->GetEnumValue("cut_strengthening_type", cut_strengthening_type,app_->prefix());
03002 double lo_inf, up_inf;
03003 app_->options()->GetNumericValue("nlp_lower_bound_inf",lo_inf, app_->prefix());
03004 app_->options()->GetNumericValue("nlp_upper_bound_inf",up_inf, app_->prefix());
03005 infty_ = std::min(fabs(lo_inf), fabs(up_inf));
03006
03007 #ifdef COIN_HAS_FILTERSQP
03008 is_given =
03009 #endif
03010 app_->options()->GetNumericValue("resolve_on_small_infeasibility", infeasibility_epsilon_, app_->prefix());
03011
03012 #ifdef COIN_HAS_FILTERSQP
03013 if(filter && !is_given){
03014
03015 infeasibility_epsilon_ = 1e-06;
03016 std::string o_name = app_->prefix();
03017 o_name += "resolve_on_small_infeasibility";
03018 app_->options()->SetNumericValue(o_name.c_str(), infeasibility_epsilon_, true, true);
03019 }
03020 #endif
03021 if (cut_strengthening_type != CS_None) {
03022
03023 cutStrengthener_ = new CutStrengthener(app_->clone(), app_->options());
03024 }
03025 }
03026 }
03027
03028 void
03029 OsiTMINLPInterface::SetStrongBrachingSolver(SmartPtr<StrongBranchingSolver> strong_branching_solver)
03030 {
03031 strong_branching_solver_ = strong_branching_solver;
03032 }
03033
03034
03035 #ifdef STRONG_COMPARE
03036 static double objorig;
03037 #endif
03038
03039 void
03040 OsiTMINLPInterface::markHotStart()
03041 {
03042 if (IsValid(strong_branching_solver_)) {
03043 #ifdef STRONG_COMPARE
03044
03045 OsiSolverInterface::markHotStart();
03046 objorig = getObjValue();
03047 #endif
03048 optimizationStatusBeforeHotStart_ = optimizationStatus_;
03049 strong_branching_solver_->markHotStart(this);
03050 }
03051 else {
03052
03053 OsiSolverInterface::markHotStart();
03054 }
03055 }
03056
03057 void
03058 OsiTMINLPInterface::solveFromHotStart()
03059 {
03060 if (IsValid(strong_branching_solver_)) {
03061 #ifdef STRONG_COMPARE
03062
03063 OsiSolverInterface::solveFromHotStart();
03064 double obj_nlp = getObjValue() - objorig;
03065 #endif
03066 optimizationStatus_ = strong_branching_solver_->solveFromHotStart(this);
03067 hasBeenOptimized_ = true;
03068 #ifdef STRONG_COMPARE
03069 double obj_other = getObjValue() - objorig;
03070 printf("AWDEBUG: Strong Branching results: NLP = %15.8e Other = %15.8e\n",
03071 obj_nlp, obj_other);
03072 #endif
03073 }
03074 else {
03075
03076 OsiSolverInterface::solveFromHotStart();
03077 }
03078 }
03079
03080 void
03081 OsiTMINLPInterface::unmarkHotStart()
03082 {
03083 if (IsValid(strong_branching_solver_)) {
03084 #ifdef STRONG_COMPARE
03085
03086 OsiSolverInterface::unmarkHotStart();
03087 #endif
03088 strong_branching_solver_->unmarkHotStart(this);
03089 optimizationStatus_ = optimizationStatusBeforeHotStart_;
03090 }
03091 else {
03092
03093 OsiSolverInterface::unmarkHotStart();
03094 }
03095 }
03096
03097 const double * OsiTMINLPInterface::getObjCoefficients() const
03098 {
03099 const int n = getNumCols();
03100 delete [] obj_;
03101 obj_ = NULL;
03102 obj_ = new double[n];
03103
03104 bool new_x = true;
03105 const double* x_sol = problem_->x_sol();
03106 bool retval = problem_->eval_grad_f(n, x_sol, new_x, obj_);
03107
03108 if (!retval) {
03109
03110 fprintf(stderr, "ERROR WHILE EVALUATING GRAD_F in OsiTMINLPInterface::getObjCoefficients()\n");
03111 delete [] obj_;
03112 obj_ = NULL;
03113 }
03114
03115 return obj_;
03116 }
03117
03119 void
03120 OsiTMINLPInterface::use(Ipopt::SmartPtr<TMINLP2TNLP> tminlp2tnlp){
03121 problem_ = tminlp2tnlp;
03122 problem_to_optimize_ = GetRawPtr(problem_);
03123 feasibilityProblem_->use(GetRawPtr(tminlp2tnlp));}
03124
03125 }