source: branches/devel/Bonmin/src/IpoptInterface/TMINLP2TNLP.cpp @ 44

Last change on this file since 44 was 44, checked in by pbonami, 13 years ago

Can add linear cut to nonlinear formulation in Ipopt

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 16.3 KB
Line 
1// (C) Copyright International Business Machines Corporation and Carnegie Mellon University 2004
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// Carl D. Laird, Carnegie Mellon University,
7// Andreas Waechter, International Business Machines Corporation
8// Pierre Bonami, Carnegie Mellon University,
9//
10// Date : 12/01/2004
11
12
13#include "TMINLP2TNLP.hpp"
14#include "IpBlas.hpp"
15#include "IpAlgTypes.hpp"
16#include <string>
17#include <fstream>
18#include <sstream>
19namespace Ipopt
20{
21  TMINLP2TNLP::TMINLP2TNLP(const SmartPtr<TMINLP> tminlp,
22      const OptionsList& options)
23      :
24      tminlp_(tminlp),
25      n_(0),
26      m_(0),
27      var_types_(NULL),
28      x_l_(NULL),
29      x_u_(NULL),
30      orig_x_l_(NULL),
31      orig_x_u_(NULL),
32      g_l_(NULL),
33      g_u_(NULL),
34      x_init_(NULL),
35      x_init_user_(NULL),
36      x_sol_(NULL),
37      g_sol_(NULL),
38      duals_sol_(NULL),
39//      return_status_(NOT_SOLVED),
40      obj_value_(1e100),
41      curr_warm_starter_(NULL),
42      need_new_warm_starter_(true)
43  {
44    // read the nlp size and bounds information from
45    // the TMINLP and keep an internal copy. This way the
46    // caller can modify the bounds that are sent to Ipopt;
47    DBG_ASSERT(IsValid(tminlp_));
48
49    Index nnz_jac_g;
50    tminlp_->get_nlp_info(n_, m_, nnz_jac_g, nnz_h_lag_, index_style_);
51
52    // Allocate space for the variable types vector
53    var_types_ = new TMINLP::VariableType[n_];
54
55    // retrieve the variable types
56    tminlp_->get_var_types(n_, var_types_);
57
58    // Allocate space for the internal copy of the variable bounds
59    x_l_ = new Number[n_];
60    x_u_ = new Number[n_];
61    orig_x_l_ = new Number[n_];
62    orig_x_u_ = new Number[n_];
63
64    g_l_ = new Number[m_];
65    g_u_ = new Number[m_];
66
67    // retrieve the variable bounds
68    tminlp_->get_bounds_info(n_, x_l_, x_u_, m_, g_l_, g_u_);
69    IpBlasDcopy(n_, x_l_, 1, orig_x_l_, 1);
70    IpBlasDcopy(n_, x_u_, 1, orig_x_u_, 1);
71
72    //    // Check that the bounds make sense compared with the variable type
73    //    for (int i=0; i<n_; i++) {
74    //      throw_exception_on_bad_variable_bound(i);
75    //    }
76
77    // Allocate space for the initial point
78    x_init_ = new Number[3*n_ + m_];
79    tminlp_->get_starting_point(n_, true, x_init_, false, NULL, NULL,
80        m_, false, NULL);
81    x_init_user_ = new Number[n_];
82    IpBlasDcopy(n_, x_init_, 1, x_init_user_, 1);
83    duals_sol_ = NULL;
84    duals_init_ = NULL;
85    // Get values for parameters
86    options.GetNumericValue("nlp_lower_bound_inf", nlp_lower_bound_inf_, "");
87    options.GetNumericValue("nlp_upper_bound_inf", nlp_upper_bound_inf_, "");
88    options.GetBoolValue("warm_start_entire_iterate",
89        warm_start_entire_iterate_, "");
90  }
91
92  TMINLP2TNLP::~TMINLP2TNLP()
93  {
94    delete [] var_types_;
95    var_types_ = NULL;
96    delete [] x_l_;
97    x_l_ = NULL;
98    delete [] x_u_;
99    x_u_ = NULL;
100    delete [] orig_x_l_;
101    orig_x_l_ = NULL;
102    delete [] orig_x_u_;
103    orig_x_u_ = NULL;
104    delete [] g_l_;
105    g_l_ = NULL;
106    delete [] g_u_;
107    g_u_ = NULL;
108    delete [] x_init_;
109    x_init_ = NULL;
110    delete [] x_init_user_;
111    x_init_user_ = NULL;
112    delete [] x_sol_;
113    x_sol_ = NULL;
114    delete [] g_sol_;
115    g_sol_ = NULL;
116    delete [] duals_sol_;
117    duals_sol_ = NULL;
118  }
119
120  /** Copies the modification made to TNLP by the user (modifications
121      such as changing bound changing starting point,...).
122      this and other should be two instances of the same problem
123      I am trying to mimic a copy construction for Cbc
124      use with great care not safe.
125  */
126  void
127  TMINLP2TNLP::copyUserModification(TMINLP2TNLP& other)
128  {
129    DBG_ASSERT(x_l_);
130    DBG_ASSERT(x_u_);
131    DBG_ASSERT(other.x_l_);
132    DBG_ASSERT(other.x_u_);
133    DBG_ASSERT(n_ = other.n_);
134    DBG_ASSERT(m_ = other.m_);
135
136    IpBlasDcopy(n_, other.x_l_, 1, x_l_, 1);
137    IpBlasDcopy(n_, other.x_u_, 1, x_u_, 1);
138
139    if(other.duals_init_) {
140      duals_init_ = &x_init_[n_];
141      IpBlasDcopy(3*n_+ m_, other.x_init_, 1, x_init_, 1);
142    }
143    else
144      IpBlasDcopy(n_, other.x_init_, 1, x_init_, 1);
145    return_status_ = other.return_status_;
146
147    if(other.x_sol_ !=NULL) {
148      //DBG_ASSERT(return_status_ != NOT_SOLVED);
149      Set_x_sol(n_,other.x_sol_);
150    }
151
152    if(other.g_sol_!=NULL) {
153//      DBG_ASSERT(return_status_ != NOT_SOLVED);
154      g_sol_ = new Number [m_];
155      IpBlasDcopy(m_, other.g_sol_, 1, g_sol_, 1);
156    }
157
158    if(other.duals_sol_!=NULL) {
159//      DBG_ASSERT(return_status_ != NOT_SOLVED);
160      duals_sol_ = new Number[m_ + 2*n_];
161      IpBlasDcopy(2*n_+ m_, other.duals_sol_, 1, duals_sol_, 1);
162    }
163
164    obj_value_ = other.obj_value_;
165
166    curr_warm_starter_ = other.curr_warm_starter_;
167
168    nlp_lower_bound_inf_ = other.nlp_lower_bound_inf_;
169
170    nlp_upper_bound_inf_ = other.nlp_upper_bound_inf_;
171
172    need_new_warm_starter_ = other.need_new_warm_starter_;
173  }
174
175  void TMINLP2TNLP::SetVariableBounds(Index var_no, Number x_l, Number x_u)
176  {
177    DBG_ASSERT(var_no >= 0 && var_no < n_);
178    x_l_[var_no] = x_l;
179    x_u_[var_no] = x_u;
180    //    throw_exception_on_bad_variable_bound(var_no);
181  }
182
183  void TMINLP2TNLP::SetVariableLowerBound(Index var_no, Number x_l)
184  {
185    DBG_ASSERT(var_no >= 0 && var_no < n_);
186    x_l_[var_no] = x_l;
187    //    throw_exception_on_bad_variable_bound(var_no);
188  }
189
190  void TMINLP2TNLP::SetVariableUpperBound(Index var_no, Number x_u)
191  {
192    DBG_ASSERT(var_no >= 0 && var_no < n_);
193    x_u_[var_no] = x_u;
194    //    throw_exception_on_bad_variable_bound(var_no);
195  }
196
197  void TMINLP2TNLP::SetStartingPoint(Index n,const Number* x_init)
198  {
199    DBG_ASSERT(n == n_);
200    IpBlasDcopy(n_, x_init, 1, x_init_, 1);
201  }
202
203  void TMINLP2TNLP::resetStartingPoint()
204  {
205    curr_warm_starter_ = NULL;
206    IpBlasDcopy(n_, x_init_user_, 1, x_init_, 1);
207  }
208
209  void TMINLP2TNLP::setxInit(Index ind, const Number val)
210  {
211    x_init_[ind] = val;
212  }
213
214  void TMINLP2TNLP::setxInit(Index n,const Number* x_init)
215  {
216    DBG_ASSERT(n == n_);
217    IpBlasDcopy(n_, x_init, 1, x_init_, 1);
218  }
219
220  void TMINLP2TNLP::setDualInit(Index ind, const Number val)
221  {
222    if(!duals_init_)
223      duals_init_ = &x_init_[n_];
224    duals_init_[ind] = val;
225  }
226
227  void TMINLP2TNLP::setDualsInit(Index n, const Number* duals_init)
228  {
229    DBG_ASSERT(n == m_ );
230
231    if(!duals_init_)
232      duals_init_ = &x_init_[n_];
233
234    IpBlasDcopy(n, duals_init, 1, duals_init_, 1);
235
236  }
237
238  /** Set the contiuous solution */
239  void TMINLP2TNLP::Set_x_sol(Index n, const Number* x_sol)
240  {
241    DBG_ASSERT(n == n_);
242    if (!x_sol_) {
243      x_sol_ = new Number[n];
244    }
245    IpBlasDcopy(n_, x_sol, 1, x_sol_, 1);
246  }
247
248  /** Change the type of the variable */
249  void TMINLP2TNLP::SetVariableType(Index n, TMINLP::VariableType type)
250  {
251    DBG_ASSERT(n >= 0 && n < n_);
252    var_types_[n] = type;
253  }
254
255  bool TMINLP2TNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
256      Index& nnz_h_lag, TNLP::IndexStyleEnum& index_style)
257  {
258    bool return_code = tminlp_->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
259    m += tminlp_->nLinearCuts_;
260    nnz_jac_g += tminlp_->linearCutsNnz_;
261    return return_code;
262  }
263
264  bool TMINLP2TNLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
265      Index m, Number* g_l, Number* g_u)
266  {
267    IpBlasDcopy(n, x_l_, 1, x_l, 1);
268    IpBlasDcopy(n, x_u_, 1, x_u, 1);
269    IpBlasDcopy(m_, g_l_, 1, g_l, 1);
270    IpBlasDcopy(m_, g_u_, 1, g_u, 1);
271    IpBlasDcopy(tminlp_->nLinearCuts_, tminlp_->lower_, 1, &g_l[m_],1); 
272    IpBlasDcopy(tminlp_->nLinearCuts_, tminlp_->upper_, 1, &g_u[m_],1); 
273    return true;
274  }
275
276  bool TMINLP2TNLP::get_starting_point(Index n, bool init_x, Number* x,
277      bool init_z, Number* z_L, Number* z_U,
278      Index m, bool init_lambda,
279      Number* lambda)
280  {
281    assert(m==m_ + tminlp_->nLinearCuts_);
282    assert(n==n_);
283    if (init_x == true) {
284      if(x_init_==NULL)
285        return false;
286      IpBlasDcopy(n, x_init_, 1, x, 1);
287    }
288    if (init_z == true) {
289      if(duals_init_ == NULL)
290        return false;
291      IpBlasDcopy(n, &duals_init_[m], 1, z_L, 1);
292      IpBlasDcopy(n, &duals_init_[m + n], 1, z_U, 1);
293
294    }
295    if(init_lambda == true) {
296      if(duals_init_ == NULL)
297        return false;
298      IpBlasDcopy(m_, duals_init_, 1, lambda, 1);
299      for(int i = 0 ; i < tminlp_->nLinearCuts_; i++)
300      {
301        lambda [i + m_ ] = 0.;
302      }
303    }
304
305    need_new_warm_starter_ = true;
306    return true;
307  }
308
309  bool TMINLP2TNLP::get_warm_start_iterate(IteratesVector& warm_start_iterate)
310  {
311    if (IsNull(curr_warm_starter_)) {
312      return false;
313    }
314
315    bool retval = curr_warm_starter_->WarmStartIterate(n_, x_l_, x_u_,
316        warm_start_iterate);
317
318    need_new_warm_starter_ = true;
319    return retval;
320  }
321
322  bool TMINLP2TNLP::eval_f(Index n, const Number* x, bool new_x,
323      Number& obj_value)
324  {
325    return tminlp_->eval_f(n, x, new_x, obj_value);
326  }
327
328  bool TMINLP2TNLP::eval_grad_f(Index n, const Number* x, bool new_x,
329      Number* grad_f)
330  {
331    return tminlp_->eval_grad_f(n, x, new_x, grad_f);
332  }
333
334  bool TMINLP2TNLP::eval_g(Index n, const Number* x, bool new_x,
335      Index m, Number* g)
336  {
337    int return_code = tminlp_->eval_g(n, x, new_x, m_, g);
338    // Add the linear cuts
339    int nnz = 0;
340    for(int i = 0 ; i < tminlp_->nLinearCuts_ ; i++)
341    {
342      int iplusm = m_ + i;
343      g[iplusm] = 0;
344      while(tminlp_->iRow_[nnz]==i) { g[iplusm] += tminlp_->elems_[nnz] * x[tminlp_->jCol_[nnz]]; nnz++;}
345    }
346    return return_code;
347  }
348
349  bool TMINLP2TNLP::eval_jac_g(Index n, const Number* x, bool new_x,
350      Index m, Index nele_jac, Index* iRow,
351      Index *jCol, Number* values)
352  {
353    int return_code = tminlp_->eval_jac_g(n, x, new_x, m_ , nele_jac - tminlp_->nLinearCuts_, 
354                                          iRow, jCol, values);
355    bool isFortran = index_style_ == TNLP::FORTRAN_STYLE;
356    if(iRow != NULL)
357    {
358      DBG_ASSERT(iCol != NULL);
359      DBG_ASSERT(values == NULL)
360     
361      int nnz = nele_jac - tminlp_->linearCutsNnz_ ;
362      for(int i = 0; i < tminlp_->linearCutsNnz_ ; i++ , nnz++)
363      {
364         iRow[nnz] = tminlp_->iRow_[i] + m_ + isFortran; 
365         jCol[nnz] = tminlp_->jCol_[i] + isFortran;
366      }
367    }
368    else
369    {
370      DBG_ASSERT(iCol == NULL);
371      DBG_ASSERT(values != NULL);
372      IpBlasDcopy(tminlp_->linearCutsNnz_ , tminlp_->elems_, 1, &values[nele_jac - tminlp_->linearCutsNnz_],1);
373    }
374  return return_code;
375  }
376
377  bool TMINLP2TNLP::eval_h(Index n, const Number* x, bool new_x,
378      Number obj_factor, Index m, const Number* lambda,
379      bool new_lambda, Index nele_hess,
380      Index* iRow, Index* jCol, Number* values)
381  {
382    return tminlp_->eval_h(n, x, new_x, obj_factor, m, lambda,
383        new_lambda, nele_hess,
384        iRow, jCol, values);
385  }
386
387  void TMINLP2TNLP::finalize_solution(SolverReturn status,
388      Index n, const Number* x, const Number* z_L, const Number* z_U,
389      Index m, const Number* g, const Number* lambda,
390      Number obj_value)
391  {
392    assert(n==n_);
393    assert(m==m_ + tminlp_->nLinearCuts_);
394    if (!x_sol_) {
395      x_sol_ = new Number[n];
396    }
397    IpBlasDcopy(n, x, 1, x_sol_, 1);
398
399    if(!g_sol_) {
400      g_sol_ = new Number [m];
401    }
402    IpBlasDcopy(m, g, 1, g_sol_, 1);
403    if (!duals_sol_) {
404      duals_sol_ = new Number[m + 2*n];
405    }
406    IpBlasDcopy(m, lambda, 1, duals_sol_, 1);
407
408    IpBlasDcopy(n, z_L, 1 , &duals_sol_[m], 1);
409    IpBlasDcopy(n, z_U, 1 , &duals_sol_[m + n], 1);
410
411    return_status_ = status;
412    obj_value_ = obj_value;
413
414    if (IsValid(curr_warm_starter_)) {
415      curr_warm_starter_->Finalize();
416    }
417  }
418
419  void TMINLP2TNLP::throw_exception_on_bad_variable_bound(Index i)
420  {
421    DBG_ASSERT(i >= 0 && i < n_);
422
423    if (var_types_[i] == TMINLP::BINARY) {
424      ASSERT_EXCEPTION(x_l_[i] == 0 && x_u_[i] == 1,
425          TMINLP_INVALID_VARIABLE_BOUNDS,
426          "Invalid variable bounds in TMINLP. All binaries must have 0,1 for variable bounds."
427                      );
428    }
429    else if (var_types_[i] == TMINLP::INTEGER) {
430      if(fabs(x_l_[i])<INT_MAX && fabs(x_u_[i]) < INT_MAX)//round to closest valid integer
431      {
432        int x_l = (int)ceil(x_l_[i]);
433        int x_u = (int)floor(x_u_[i]);
434        std::cerr<<"Inconsistent bounds on an integer"<<std::endl;
435        ASSERT_EXCEPTION(x_l_[i] == (Number)x_l && x_u_[i] == (Number)x_u,
436            TMINLP_INVALID_VARIABLE_BOUNDS,
437            "Invalid variable bounds in TMINLP. All integer variables must have integer bounds."
438                        );
439      }
440    }
441  }
442  bool TMINLP2TNLP::checkZeroDimension(ApplicationReturnStatus
443      &optimization_status)
444  {
445    for(int i = 0 ; i < n_ ; i++) {
446      if(x_u_[i] - x_l_[i] > 1e-5)
447        return 0;
448    }
449
450    //Problem has no variables just check if the unique solution given by the bounds is
451    // feasible or not.
452    eval_f(n_, x_l_, true, obj_value_);
453    if (!x_sol_) {
454      x_sol_ = new Number[n_];
455    }
456    IpBlasDcopy(n_, x_l_, 1, x_sol_, 1);
457
458    if(!g_sol_) {
459      g_sol_ = new Number [m_];
460    }
461    eval_g(n_, x_l_, true, m_, g_sol_);
462    optimization_status = Solve_Succeeded;
463    for(int i = 0 ; i < m_ ; i++) {
464      if(g_sol_[i] - g_l_[i] <  - 1e-07 || g_sol_[i] - g_u_[i] > 1e-07) {
465        optimization_status = Infeasible_Problem_Detected;
466        return 1;
467      }
468    }
469    return 1;
470  }
471
472  bool TMINLP2TNLP::intermediate_callback(AlgorithmMode mode,
473      Index iter, Number obj_value,
474      Number inf_pr, Number inf_du,
475      Number mu, Number d_norm,
476      Number regularization_size,
477      Number alpha_du, Number alpha_pr,
478      Index ls_trials,
479      const IpoptData* ip_data,
480      IpoptCalculatedQuantities* ip_cq)
481  {
482    // If we don't have this swtiched on, we assume that also the
483    // "warm_start" option for bonmin is set not to refer to the
484    // interior warm start object
485    if (!warm_start_entire_iterate_) {
486      return true;
487    }
488
489    if (need_new_warm_starter_) {
490      // Create a new object for later warm start information
491      curr_warm_starter_ = new IpoptInteriorWarmStarter(n_, x_l_, x_u_,
492          nlp_lower_bound_inf_,
493          nlp_upper_bound_inf_,
494          warm_start_entire_iterate_);
495      need_new_warm_starter_ = false;
496    }
497
498    return curr_warm_starter_->UpdateStoredIterates(mode, *ip_data, *ip_cq);
499  }
500
501
502  /** Procedure to ouptut relevant informations to reproduce a sub-problem.
503  Compare the current problem to the problem to solve
504  and writes files with bounds which have changed and current starting point.
505
506  */
507  void
508  TMINLP2TNLP::outputDiffs(const std::string& probName, const std::string * varNames)
509  {
510    const int &numcols = n_;
511    const int &numrows = m_;
512
513    const double * currentLower = x_l();
514    const double * currentUpper = x_u();
515
516    const double * originalLower = orig_x_l();
517    const double * originalUpper = orig_x_u();
518    CoinRelFltEq eq;
519    std::string fBoundsName = probName;
520    std::ostringstream os;
521    fBoundsName+=".bounds";
522    std::string fModName = probName;
523    fModName+= ".mod";
524    std::ofstream fBounds;
525    std::ofstream fMod;
526    bool hasVarNames = 0;
527
528    if(varNames!=NULL )
529      hasVarNames=1;
530    if(hasVarNames)
531      fMod.open(fModName.c_str());
532    fBounds.open(fBoundsName.c_str());
533
534    for(int i = 0 ; i < numcols ; i++) {
535      if(!eq(currentLower[i],originalLower[i])) {
536        if(hasVarNames)
537          fMod<<"bounds"<<i<<": "
538          <<varNames[i]<<" >= "
539          <<currentLower[i]<<";\n";
540
541
542        fBounds<<"LO"<<"\t"<<i<<"\t"<<currentLower[i]<<std::endl;
543      }
544      if(!eq(currentUpper[i],originalUpper[i])) {
545        if(hasVarNames)
546          fMod<<"bounds"<<i<<": "
547          <<varNames[i]<<" <= "
548          <<currentUpper[i]<<";\n";
549
550        fBounds<<"UP"<<"\t"<<i<<"\t"<<currentUpper[i]<<std::endl;
551      }
552    }
553
554    //write a file with starting point
555    std::string fStartPointName=probName;
556    fStartPointName+=".start";
557
558    std::ofstream fStartPoint(fStartPointName.c_str());
559    const double * primals = x_init();
560    const double * duals = duals_init();
561    fStartPoint.precision(17);
562    fStartPoint<<numcols<<"\t"<<2*numcols+numrows<<std::endl;
563    for(int i = 0 ; i < numcols ; i++)
564      fStartPoint<<primals[i]<<std::endl;
565    int end = 2*numcols + numrows;
566    if(duals) {
567      for(int i = 0 ; i < end; i++)
568        fStartPoint<<duals[i]<<std::endl;
569    }
570
571  }
572
573  /** force solution to be fractionnal.*/
574  void
575  TMINLP2TNLP::force_fractionnal_sol()
576  {
577    for(int i=0 ; i < n_ ; i++) {
578      if( ( var_types_[i] == TMINLP::INTEGER ||
579          var_types_[i] == TMINLP::BINARY )&&
580          x_l_[i] < x_u_[i] + 0.5)//not fixed
581      {
582        x_sol_[i] = ceil(x_l_[i]) + 0.5;//make it integer infeasible
583      }
584    }
585  }
586
587}
588// namespace Ipopt
589
Note: See TracBrowser for help on using the repository browser.