source: trunk/Couenne/src/problem/checkNLP.cpp @ 115

Last change on this file since 115 was 115, checked in by pbelotti, 11 years ago

fix right value of objective computed in checkNLP -- numerics can make a bound interval infeasible. Some cleanup

File size: 5.5 KB
Line 
1/*
2 * Name:    checkNLP.cpp
3 * Author:  Pietro Belotti
4 * Purpose: check NLP feasibility of incumbent integer solution
5 *
6 * (C) Carnegie-Mellon University, 2006-09.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#include "CoinHelperFunctions.hpp"
11#include "CouenneProblem.hpp"
12
13// check if solution is MINLP feasible
14bool CouenneProblem::checkNLP (const double *solution, double &obj, bool recompute) const {
15
16  const int infeasible = 1;
17  const int wrong_obj  = 2;
18
19  /*printf ("checking solution: [%.12e] ", obj);
20  for (int i=0; i<nOrigVars_; i++)
21    printf ("%.12e ", solution [i]);
22    printf ("\n");*/
23
24  // pre-check on original variables
25  for (int i=0; i < nOrigVars_; i++) {
26
27    CouNumber val = solution [i];
28
29    // check (original and auxiliary) variables' integrality
30
31    if ((variables_ [i] -> isInteger ()) &&
32        (variables_ [i] -> Type () == VAR) &&
33        (variables_ [i] -> Multiplicity () > 0) &&
34        (fabs (val - COUENNE_round (val)) > feas_tolerance_)) {
35
36      Jnlst()->Printf(Ipopt::J_ITERSUMMARY, J_PROBLEM,
37                      "checkNLP: integrality %d violated: %.6f [%g,%g]\n", 
38                      i, val, domain_.lb (i), domain_.ub (i));
39
40      return false;
41    }
42  }
43
44  CouNumber *sol = new CouNumber [nVars ()];
45
46  // copy solution, evaluate the corresponding aux, and then replace
47  // the original variables again for checking
48  CoinCopyN (solution, nOrigVars_, sol);
49  getAuxs (sol);
50  CoinCopyN (solution, nOrigVars_, sol);
51
52  // install NL solution candidate in evaluation structure
53  domain_.push (nVars (), sol, domain_.lb (), domain_.ub (), false);
54
55  /*printf ("checknlp: %d vars -------------------\n", domain_.current () -> Dimension ());
56  for (int i=0; i<domain_.current () -> Dimension (); i++)
57    printf ("%4d %.12e [%.12e %.12e]\n",
58    i, domain_.x (i), domain_.lb (i), domain_.ub (i));*/
59
60  expression *objBody = Obj (0) -> Body ();
61
62  // BUG: if Ipopt solution violates bounds of original variables and
63  // objective depends on originals, we may have a "computed object"
64  // out of bounds
65
66  //CouNumber realobj = (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
67  CouNumber realobj = obj;
68
69  if (objBody) 
70    realobj = 
71      (objBody -> Index () >= 0) ?
72      sol [objBody -> Index ()] : 
73      (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
74
75  /*printf ("%.12e %.12e %.12e ------------------------------\n",
76          realobj, sol [objBody -> Index ()],
77          (*(objBody -> Image () ? objBody -> Image () : objBody)) ());*/
78
79  bool retval = true;
80
81  try {
82
83    // check if objective corresponds
84
85    if (fabs (realobj - obj) / (1. + fabs (realobj)) > feas_tolerance_) {
86
87      Jnlst()->Printf(Ipopt::J_ITERSUMMARY, J_PROBLEM,
88                      "checkNLP: false objective. %g != %g (diff. %g)\n", 
89                      realobj, obj, realobj - obj);
90
91      if (!recompute)
92        throw wrong_obj;
93    }
94
95    if (recompute)
96      obj = realobj;
97
98    //printf ("recomputed: %.12e\n", obj);
99
100    for (int i=0; i < nOrigVars_; i++) {
101
102      if (variables_ [i] -> Multiplicity () <= 0) 
103        continue;
104
105      CouNumber val = domain_.x (i);
106
107      // check bounds
108
109      if ((val > domain_.ub (i) + feas_tolerance_) ||
110          (val < domain_.lb (i) - feas_tolerance_)) {
111
112        Jnlst()->Printf(Ipopt::J_ITERSUMMARY, J_PROBLEM,
113                        "checkNLP: variable %d out of bounds: %.6f [%g,%g] (diff %g)\n", 
114                        i, val, domain_.lb (i), domain_.ub (i),
115                        CoinMax (fabs (val - domain_.lb (i)), 
116                                 fabs (val - domain_.ub (i))));
117        throw infeasible;
118      }
119
120      // check (original and auxiliary) variables' integrality
121
122      if (variables_ [i] -> isInteger () &&
123          (fabs (val - COUENNE_round (val)) > feas_tolerance_)) {
124
125        Jnlst()->Printf(Ipopt::J_ITERSUMMARY, J_PROBLEM,
126                        "checkNLP: integrality %d violated: %.6f [%g,%g]\n", 
127                        i, val, domain_.lb (i), domain_.ub (i));
128
129        throw infeasible;
130      }
131
132      /*if (variables_ [i] -> Type () == AUX) {
133        printf ("checking aux ");
134        variables_ [i] -> print (); printf (" := ");
135        variables_ [i] -> Image () -> print ();
136        printf (" --- %.12e = %.12e [%.12e]; {",
137                (*(variables_ [i])) (),
138                (*(variables_ [i] -> Image ())) (),
139                (*(variables_ [i])) () -
140                (*(variables_ [i] -> Image ())) ());
141        for (int j=0; j<nVars (); j++)
142          printf ("%.12e ", (*(variables_ [j])) ());
143        printf ("}\n");
144        }*/
145
146      CouNumber delta;
147
148      // check if auxiliary has zero infeasibility
149
150      if ((variables_ [i] -> Type () == AUX) && 
151          (fabs (delta = (*(variables_ [i])) () - 
152                 (*(variables_ [i] -> Image ())) ()) > feas_tolerance_)) {
153
154        Jnlst()->Printf(Ipopt::J_ITERSUMMARY, J_PROBLEM,
155                        "checkNLP: auxiliarized %d violated (%g)\n", i, delta);
156
157        throw infeasible;
158      }
159    }
160
161    // check constraints
162
163    for (int i=0; i < nCons (); i++) {
164
165      CouenneConstraint *c = Con (i);
166
167      CouNumber
168        body = (*(c -> Body ())) (),
169        lhs  = (*(c -> Lb   ())) (),
170        rhs  = (*(c -> Ub   ())) ();
171
172      if ((rhs < COUENNE_INFINITY) &&
173           (body > rhs + feas_tolerance_ * (1 + CoinMax (fabs (body), fabs (rhs)))) || 
174          (lhs > -COUENNE_INFINITY) &&
175          (body < lhs - feas_tolerance_ * (1 + CoinMax (fabs (body), fabs (lhs))))) {
176
177        if (Jnlst()->ProduceOutput(Ipopt::J_ITERSUMMARY, J_PROBLEM)) {
178
179          Jnlst()->Printf
180            (Ipopt::J_ITERSUMMARY, J_PROBLEM,
181             "checkNLP: constraint %d violated (lhs=%+e body=%+e rhs=%+e, violation %g): ",
182             i, lhs, body, rhs, CoinMax (lhs-body, body-rhs));
183
184          c -> print ();
185        }
186
187        throw infeasible;
188      }
189    }
190  }
191
192  catch (int exception) {
193
194    switch (exception) {
195
196    case wrong_obj:
197      retval = false;
198      break;
199
200    case infeasible:
201    default:
202      retval = false;
203      break;
204    }
205  }
206
207  delete [] sol;
208  domain_.pop ();
209
210  return retval;
211}
Note: See TracBrowser for help on using the repository browser.