MeshKit  1.0
IAMilp.cpp
Go to the documentation of this file.
00001 // IAMilp.cpp
00002 // Interval Assignment for Meshkit
00003 //
00004 #include "meshkit/IAMilp.hpp"
00005 #include "meshkit/IAData.hpp"
00006 #include "meshkit/IASolution.hpp"
00007 #include <assert.h>
00008 
00009 #include <iostream>
00010 #include <sstream>
00011 #include <string>
00012 
00013 //usr/local/include/glpk.h
00014 #include "meshkit/glpk.h"
00015 #include <stdio.h>
00016 #include <math.h>
00017 
00018 class GlpRepresentation
00019 {
00020   public: 
00021   glp_prob *glp;
00022 };
00023 
00024 
00025 // to do: design description of sum-even constraints, represent those in the MILP
00026 
00027 // constructor
00028 IAMilp::IAMilp(const IAData *data_ptr, IASolution *solution_ptr): 
00029 data(data_ptr), solution(solution_ptr), lp(NULL), naturalConstraintsSet(false), deltasSet(false),
00030 weightsSet(false), maxesSet(false), solved(false),
00031 xStart(1), deltaStart(1), constraint_tolerance(1.0e-4), integrality_tolerance(1.0e-4),
00032 debugging(true) //zzyk true
00033 {
00034   assert(data_ptr);
00035   assert(solution_ptr);
00036   printf("\nIA MILP Problem size:\n");
00037   printf("  number of variables: %lu\n", data->I.size());
00038   printf("  number of constraints: %lu\n\n", data->constraints.size());
00039 }
00040 
00041 
00042 IAMilp::~IAMilp() {data = NULL;}
00043 
00044 // index of first delta variable
00045 // organized by x_i - I_i = dp_i - dn_i
00046 // where col(x_i) = i+1, col(dp_i) = i + I_i + 1, col(dn_i) = i + 2*I_i +1
00047 
00048 inline 
00049 int IAMilp::x_i(int i)
00050 {
00051   return i+xStart;
00052 }
00053 
00054 inline 
00055 int IAMilp::delta_plus_i(int i) 
00056 { 
00057   return i + deltaStart;
00058 }
00059 
00060 inline 
00061 int IAMilp::delta_minus_i(int i) 
00062 { 
00063   return i + deltaStart + data->num_variables(); 
00064 }
00065 
00066 inline 
00067 int IAMilp::delta_j(int j) 
00068 { 
00069   return j + deltaConstraintStart;
00070 }
00071 
00072 
00073 bool IAMilp::create_problem()
00074 {
00075   lp = new GlpRepresentation();
00076   lp->glp = glp_create_prob();
00077   if (!lp || !lp->glp) 
00078     return false;
00079   if (debugging) {
00080     glp_set_prob_name(lp->glp, "IA MILP");
00081     glp_set_obj_name(lp->glp, "interval deviations");
00082   }
00083   glp_set_obj_dir(lp->glp, GLP_MIN); // minimize objective function
00084   return true;
00085 }
00086 
00087 bool IAMilp::destroy_problem()
00088 {
00089   naturalConstraintsSet = false;
00090   deltasSet = false;
00091   weightsSet = false;
00092   maxesSet = false;
00093   solved = false;
00094   glp_delete_prob(lp->glp);
00095   delete lp;
00096   lp = NULL;
00097   return true;
00098 }
00099 
00100 int glpk_bound_type(double low, double high)
00101 {
00102   if (low < high) 
00103   {
00104     if (low == MESHKIT_IA_lowerUnbound)
00105       if (high == MESHKIT_IA_upperUnbound )
00106         return GLP_FR;  // -inf, inf
00107       else
00108         return GLP_UP;  // -inf, h
00109       else
00110         if (high == MESHKIT_IA_upperUnbound )
00111           return GLP_LO;  // l, inf
00112         else
00113           return  GLP_DB; // l, h
00114   }
00115   return GLP_FX; // equality l = h
00116  }
00117 
00118 bool IAMilp::set_natural_constraints()
00119 {
00120   // don't call this twice
00121   if (lp==NULL)
00122     return false;
00123   if (naturalConstraintsSet)
00124     return true;
00125   
00126   
00127   // allocate columns for variables (curves)
00128   // and set upper and lower bounds on them
00129   // cols must be allocated before the rows reference their indices
00130   int nv = data->num_variables(); 
00131   xStart = glp_add_cols(lp->glp, nv); // set here, used by x_i(i)
00132   for (int i=0; i<data->num_variables(); ++i)
00133   {
00134     if (debugging)
00135     {
00136       std::stringstream ss;
00137       ss << "x" << x_i(i);
00138       glp_set_col_name(lp->glp, x_i(i), ss.str().c_str());
00139     }
00140     glp_set_col_bnds(lp->glp, x_i(i), GLP_LO, 1., MESHKIT_IA_upperUnbound); // 1, inf    
00141     // enforce integer solution
00142     // glp_set_col_kind sets (changes) the kind of j-th column (structural variable) as specified by the parameter kind:
00143     // GLP_CV = continuous variable; GLP_IV = integer variable; GLP_BV = binary variable.
00144     glp_set_col_kind(lp->glp, x_i(i), GLP_IV);
00145   }
00146 
00147   // allocate rows for constraints
00148   const int num_rows = (int) data->constraints.size();
00149   const int next_free_row = glp_add_rows(lp->glp, num_rows);
00150   assert(next_free_row==1);
00151   
00152   // upper and lower bounds on constraints
00153   std::vector<int> ind;
00154   std::vector<double> val;
00155   for (int j=0; j<data->constraints.size(); ++j)
00156   {
00157     if (debugging)
00158     {
00159       std::stringstream ss;
00160       ss << "constraint " << j+1;
00161       glp_set_row_name(lp->glp, j+1, ss.str().c_str());
00162     }
00163     const int glp_j = next_free_row+j;
00164     const double low = data->constraints[j].lowerBound;
00165     const double high = data->constraints[j].upperBound;
00166     const int bound_type = glpk_bound_type( low, high );
00167     glp_set_row_bnds(lp->glp, glp_j, bound_type, low, high);
00168     
00169     // coefficients
00170     const int nnz = (int) data->constraints[j].M.size();
00171     ind.reserve(nnz);
00172     val.reserve(nnz);
00173     for (int k = 0; k<nnz; ++k)
00174     {
00175       ind[k] = x_i(data->constraints[j].M[k].col);
00176       val[k] = data->constraints[j].M[k].val;
00177     }
00178     // -1 since glp indexes from 1...
00179     glp_set_mat_row(lp->glp, glp_j, nnz, ind.data()-1, val.data()-1);    
00180   }
00181     
00182   naturalConstraintsSet = true;
00183   return true;
00184 }
00185 
00186 void IAMilp::get_x_bounds(const int i, double &lobound, double &hibound)
00187 {
00188   const double xbound = solution->x_solution[i]; 
00189   const double Ibound = data->I[i];
00190   if (xbound < Ibound)
00191   {
00192     lobound = xbound;
00193     hibound = Ibound;
00194   }
00195   else
00196   {
00197     lobound = Ibound;
00198     hibound = xbound;
00199   }
00200 }
00201 
00202 bool IAMilp::set_deltas()
00203 {
00204   if (!naturalConstraintsSet)
00205     return false;
00206   if (deltasSet)
00207     return true;
00208   
00209   // allocate columns for deltas 
00210   // and set upper and lower bounds on them
00211   int nv = 2*data->num_variables(); 
00212   deltaStart = glp_add_cols(lp->glp, nv);
00213   for (int i=0; i<data->num_variables(); ++i)
00214   {
00215     if (debugging)
00216     {
00217       std::stringstream ssp;
00218       ssp << "delta_plus_x" << x_i(i);
00219       glp_set_col_name(lp->glp, delta_plus_i(i), ssp.str().c_str());
00220 
00221       std::stringstream ssm;
00222       ssm << "delta_minus_x" << x_i(i);
00223       glp_set_col_name(lp->glp, delta_minus_i(i), ssm.str().c_str());
00224     }
00225     glp_set_col_bnds(lp->glp, delta_plus_i(i), GLP_LO, 0., MESHKIT_IA_upperUnbound); // 0, inf
00226     glp_set_col_bnds(lp->glp, delta_minus_i(i), GLP_LO, 0., MESHKIT_IA_upperUnbound); // 0, inf
00227   }
00228   
00229   // allocate rows for delta-setting constraints
00230   const int num_rows = data->num_variables();
00231   deltaConstraintStart = glp_add_rows(lp->glp, num_rows);
00232   for (int j=0; j<num_rows; ++j)
00233   {
00234     if (debugging)
00235     {
00236       std::stringstream ss;
00237       // index is correct, contraint for the variable == column of x_j
00238       ss << "row" << delta_j(j) << "_delta_constraints_for_x" << x_i(j); 
00239       glp_set_row_name(lp->glp, delta_j(j), ss.str().c_str());
00240     }
00241     // old way:
00242     // x_i - I_i = p_i - m_i  <=> x_i - p_i + m_i = I_i
00243     // rhs bound
00244     // const double bound = data->I[j]; 
00245 
00246     // new way: 
00247     // base delta on passed in nlp x_solution *and* original problem goals
00248     // we trust and use the nlp solution value
00249     // lower <= x - delta_plus + delta_minus <= upper
00250     // and we have non-zero deltas only when the milp solution is outside the interval [x_nlp,I]
00251     
00252     // to do: not sure if the above is sufficient, I may need to have a small weight to push variables back towards the goals, 
00253     // to decide who gets to reap the improvement when it is possible to reduce several back 
00254     // towards the goal.
00255     double lobound, hibound;
00256     get_x_bounds(j,lobound,hibound);
00257     glp_set_row_bnds(lp->glp, delta_j(j), GLP_DB, lobound, hibound);
00258     
00259     // coefficients
00260     const int ind[3] = {x_i(j), delta_plus_i(j), delta_minus_i(j)};
00261     const double val[3] = {1., -1., 1.};
00262     // -1 since indexes from 1
00263     glp_set_mat_row(lp->glp, delta_j(j), 3, ind-1, val-1);
00264   }
00265   
00266   deltasSet = true;
00267   return true;
00268 }
00269 
00270 bool IAMilp::set_sum_even_constratins()
00271 {
00272   return true;
00273 }
00274 
00275 // old paper started with an integer solution, but not us, so the bounds here will not be identical.
00276 bool IAMilp::set_bounds_1() 
00277 {
00278   // from old paper:
00279   // Curve intervals can increase by at most one
00280   if (debugging) printf("MILP bounds 1: [floor x, ceil x + 1]\n");
00281   
00282   for (int i=0; i<data->num_variables(); ++i)
00283   {
00284     const double x_nlp = solution->x_solution[i];  // nlp solution
00285     // const double Ibound = data->I[i];
00286     double lo = floor(x_nlp);
00287     if (lo<1.) lo = 1.;
00288     const double hi = ceil(x_nlp + 0.99);
00289     glp_set_col_bnds(lp->glp, x_i(i), GLP_DB, lo, hi);
00290   }
00291   
00292   return true;
00293 }
00294 
00295 bool IAMilp::set_bounds_2() 
00296 {
00297   // from old paper:
00298   // Curve intervals can increase or decrease by one
00299   if (debugging) printf("MILP bounds 2: [floor x-1, ceil x+1]\n");
00300   
00301   for (int i=0; i<data->num_variables(); ++i)
00302   {
00303     const double x_nlp = solution->x_solution[i];  // nlp solution
00304     // const double Ibound = data->I[i];
00305     double lo = floor(x_nlp - 0.99);
00306     if (lo<1.) lo = 1.;
00307     const double hi = ceil(x_nlp + 0.99);
00308     glp_set_col_bnds(lp->glp, x_i(i), GLP_DB, lo, hi);
00309   }
00310   
00311   return true;
00312 }
00313 
00314 bool IAMilp::set_bounds_3() 
00315 {
00316   // from old paper:
00317   // Curve intervals can double, but can’t decrease.
00318   if (debugging) printf("MILP bounds 3: [floor x, ceil 2x]\n");
00319   
00320   for (int i=0; i<data->num_variables(); ++i)
00321   {
00322     const double x_nlp = solution->x_solution[i];  // nlp solution
00323     // const double Ibound = data->I[i];
00324     double lo = floor(x_nlp);
00325     if (lo<1.) lo = 1.;
00326     const double hi = ceil(x_nlp * 2);
00327     glp_set_col_bnds(lp->glp, x_i(i), GLP_DB, lo, hi);
00328   }
00329   
00330   return true;
00331 }
00332 
00333 bool IAMilp::set_bounds_4() 
00334 {
00335   // from old paper:
00336   // Curve intervals can double, and can decrease by at most one. Interval-sum variables are bounded between
00337   // the floor and twice the ceiling of the pseudo-relaxed solution.
00338   if (debugging) printf("MILP bounds 4: [floor x - 1, ceil 2x]\n");
00339 
00340   
00341   for (int i=0; i<data->num_variables(); ++i)
00342   {
00343     const double x_nlp = solution->x_solution[i];  // nlp solution
00344     // const double Ibound = data->I[i];
00345     double lo = floor(x_nlp - 0.99);
00346     if (lo<1.) lo = 1.;
00347     const double hi = ceil(x_nlp * 2);
00348     glp_set_col_bnds(lp->glp, x_i(i), GLP_DB, lo, hi);
00349   }
00350   
00351   return true;
00352 }
00353 
00354 bool IAMilp::set_bounds_A() 
00355 {
00356   // new idea: bound the change based on the number of curves involved in the constraint.  
00357   if (debugging) printf("MILP bounds A: [floor x - c, ceil x +c], where c is max number of curves in a constraint involving x.\n");
00358   
00359   // count the number of edges each curves constraint list
00360   std::vector<int> num_same_side( data->num_variables(), 0 );
00361   std::vector<int> num_opposite_side( data->num_variables(), 0 );
00362   // allocate rows for constraints
00363   const int num_constraints = (int) data->constraints.size();
00364   for (int j=0; j<num_constraints; ++j)
00365   {
00366     // number of non-zeros in the constraint
00367     const int nnz = (int) data->constraints[j].M.size();
00368     // number of constraints on each side
00369     int num_plus(0), num_minus(0);
00370     for (int k = 0; k<nnz; ++k)
00371     {
00372       if (data->constraints[j].M[k].val > 0)
00373         ++num_plus;
00374       else
00375         ++num_minus;
00376     }
00377     for (int k = 0; k<nnz; ++k)
00378     {
00379       int i = data->constraints[j].M[k].col;
00380       int ss, os;
00381       if (data->constraints[j].M[k].val > 0)
00382       {
00383         ss = num_plus; 
00384         os = num_minus;
00385       }
00386       else
00387       {
00388         os = num_plus; 
00389         ss = num_minus;
00390       }
00391       if (num_same_side[i] < ss)
00392         num_same_side[i] = ss;
00393       if (num_opposite_side[i] < os)
00394         num_opposite_side[i] = os;
00395     }
00396   }  
00397   for (int i=0; i<data->num_variables(); ++i)
00398   {
00399     const double x_nlp = solution->x_solution[i];  // nlp solution
00400     double lo = floor(x_nlp - num_opposite_side[i] - num_same_side[i]); // to do: work out some reasonable bounds here
00401     if (lo<1.) lo = 1.;
00402     double hi = ceil(x_nlp + num_opposite_side[i] + num_same_side[i]);  // to do: work out some reasonable bounds here
00403     glp_set_col_bnds(lp->glp, x_i(i), GLP_DB, lo, hi);
00404   }  
00405   return true;
00406 }
00407 
00408 
00409 bool IAMilp::set_bounds_B() 
00410 {
00411   // unbounded
00412   if (debugging) printf("MILP bounds B: [1, infinity]\n");
00413 
00414   for (int i=0; i<data->num_variables(); ++i)
00415   {
00416     glp_set_col_bnds(lp->glp, x_i(i), GLP_LO, 1., MESHKIT_IA_upperUnbound);
00417   }  
00418   return true;
00419 }
00420 
00421 bool IAMilp::weight_deltas_1()
00422 {
00423   if (!deltasSet)
00424     return false;
00425 
00426   // attach weights based on passed in nlp x_solution *and* original problem goals
00427   weights_minus.resize(data->num_variables());
00428   weights_plus.resize(data->num_variables());
00429   
00430   for (int i=0; i<data->num_variables(); ++i)
00431   {
00432     double lobound, hibound;
00433     get_x_bounds(i,lobound,hibound);
00434     assert( lobound + 1e-2 > 1.0 );
00435     assert( lobound <= hibound );
00436     weights_minus[i] = 1.3  / lobound;
00437     weights_plus[i] = 1.0 / hibound; // or should this be 1.0 / lobound?
00438   }
00439 
00440   weightsSet = true;
00441   return true;
00442 }
00443 
00444 bool IAMilp::set_maxes()
00445 {
00446   if (!weightsSet)
00447     return false;
00448   
00449   // Mp = max ( delta_plus )
00450   // Mm = max ( delta_minus )
00451   // M = max( Mp, Mm )
00452   
00453   // Mwp = max ( w_p * delta_plus )
00454   // Mwm = max ( w_m * delta_minus )
00455   // Mw = max( Mwp, Mwm )
00456 
00457   // allocate rows for max-delta-setting constraints
00458   const int num_rows = 2*data->num_variables()+2;
00459   mwp_j = glp_add_rows(lp->glp, num_rows);
00460   mwm_j = mwp_j + data->num_variables();
00461   mw_j = mwp_j + 2 * data->num_variables();
00462   
00463   mwp_i = glp_add_cols(lp->glp, 3);
00464   mwm_i = mwp_i + 1;
00465   mw_i = mwp_i + 2;
00466   
00467   glp_set_col_bnds( lp->glp, mwp_i, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00468   glp_set_col_bnds( lp->glp, mwm_i, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00469   glp_set_col_bnds( lp->glp, mw_i, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00470   
00471   if (debugging)
00472   {
00473     std::stringstream ss;
00474     ss << "var max weighted deltas plus";
00475     glp_set_col_name(lp->glp, mwp_i, ss.str().c_str());
00476     ss.clear(); ss.str("");
00477     ss << "var max weighted deltas minus";
00478     glp_set_col_name(lp->glp, mwm_i, ss.str().c_str());
00479     ss.clear(); ss.str("");
00480     ss << "var max weighted deltas plus and minus";
00481     glp_set_col_name(lp->glp, mw_i, ss.str().c_str());
00482 
00483     for (int i = 0; i<data->num_variables(); ++i)
00484     {
00485       ss.clear(); ss.str("");
00486       int r = mwp_j + i;
00487       ss << "row" << r << " max weighted delta_plus" << x_i(i) << " constraint";
00488       glp_set_row_name(lp->glp, r, ss.str().c_str());
00489       ss.clear(); ss.str("");
00490       r = mwm_j + i;
00491       ss << "row" << r << " max weighted delta_minus" << x_i(i) << " constraint";
00492       glp_set_row_name(lp->glp, r, ss.str().c_str());
00493     }
00494     ss.clear(); ss.str("");
00495     int r = mw_j; 
00496     ss << "row" << r << " weighted_max_delta_ge_deltaplus";
00497     glp_set_row_name(lp->glp, r, ss.str().c_str());
00498     ss.clear(); ss.str(""); 
00499     ++r; 
00500     ss << "row" << r << " weighted_max_delta_ge_deltaminus";
00501     glp_set_row_name(lp->glp, r, ss.str().c_str());
00502   }
00503 
00504   for (int i = 0; i<data->num_variables(); ++i)
00505   {
00506     // mwp - dp_i >= 0
00507     glp_set_row_bnds(lp->glp, mwp_j + i, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00508 
00509     // mwn - dp_i >= 0
00510     glp_set_row_bnds(lp->glp, mwm_j + i, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00511 
00512     // coefficients
00513     const int indp[2] = {mwp_i, delta_plus_i(i)};
00514     const double valp[2] = {1., -weights_plus[i] };
00515     // -1 since indexes from 1
00516     glp_set_mat_row(lp->glp, mwp_j + i, 2, indp-1, valp-1);
00517 
00518     // coefficients
00519     const int indm[2] = {mwm_i, delta_minus_i(i)};
00520     const double valm[2] = {1., -weights_minus[i]};
00521     // -1 since indexes from 1
00522     glp_set_mat_row(lp->glp, mwm_j + i, 2, indm-1, valm-1);
00523   }
00524 
00525   // mw > mwp, mwm
00526   glp_set_row_bnds(lp->glp, mw_j, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00527   const double val[2] = {1., -1.};
00528   const int indpp[2] = {mw_i, mwp_i};
00529   glp_set_mat_row(lp->glp, mw_j, 2, indpp-1, val-1);
00530   glp_set_row_bnds(lp->glp, mw_j+1, GLP_LO, 0., MESHKIT_IA_upperUnbound);
00531   const int indmm[2] = {mw_i, mwm_i};
00532   glp_set_mat_row(lp->glp, mw_j+1, 2, indmm-1, val-1);
00533 
00534   maxesSet = true;
00535   return true;
00536   
00537 }
00538 
00539 bool IAMilp::set_objectives_1()
00540 {
00541   // to do: set objective function
00542   // to do: set bounds on x_i for some strategy
00543   
00544   /*
00545   // sum of deltas for testing
00546   for (int i = 0; i<data->num_variables(); ++i)
00547   {
00548     glp_set_obj_coef(lp->glp, delta_plus_i(i),  1.0);    
00549     glp_set_obj_coef(lp->glp, delta_minus_i(i), 1.0);    
00550   }
00551   */
00552   
00553   // obj = (num_variables + 1) * max_weighted_detlas + sum weighted_deltas
00554   glp_set_obj_coef(lp->glp, mw_i, 1 + data->num_variables() );
00555   // sum : this is the part that makes the MILP slow
00556   // to do: experiment with lex min max, as that might still be faster than branch and cut
00557   for (int i = 0; i<data->num_variables(); ++i)
00558   {
00559 //    glp_set_obj_coef(lp->glp, delta_plus_i(i),  weights_plus[i]);    
00560 //    glp_set_obj_coef(lp->glp, delta_minus_i(i), weights_minus[i]);    
00561   }
00562   
00563   return true;
00564 }
00565 
00566 
00567 bool IAMilp::glpk_solve(bool &optimal)
00568 {
00569   assert(!solved); // already solved?
00570   optimal = false;
00571   
00572   // print problem we are solving
00573   if (debugging)
00574     glp_write_lp(lp->glp, 0, "zzykoutput"); // only works in command line, not within xcode
00575 
00576   
00577   // to do: identify independent sub-problems and solve them separately
00578   // to do: sum-even constraints
00579   // to do: figure out how to tell it about the feasible relaxed solution to the nlp, or to limit its time searching for the simplex solution
00580   
00581   // MILP branch and cut: limit of 40 variables or so is practical for time
00582   // time grows super-linearly
00583   // this is for the lp time limit I think
00584   //void lpx_set_real_parm(LPX *lp, int parm, double val); val is milliseconds?
00585   // see tm_lim instead
00586   lpx_set_real_parm( lp->glp, LPX_K_TMLIM, 10. );  // this doesn't seem to do anything
00587   // simplex for relaxed solution
00588 //  bool relaxed_success = glp_simplex(lp->glp, NULL) == 0;
00589 /*   {
00590    solved = true;
00591    return true;
00592    }
00593    */
00594   
00595   
00596   // find integer solution using branch and cut
00597   glp_iocp parm;
00598   glp_init_iocp(&parm);
00599   if (!debugging)
00600     parm.msg_lev = GLP_MSG_OFF;
00601   // parm.pp_tech = GLP_PP_ALL; 
00602 //  parm.presolve = GLP_OFF; // not needed if do explicit solve above
00603   parm.presolve = GLP_ON; // solve relaxed solution first, remove redundant constraints, fixed variables, ...
00604   // parm.mip_gap = 0; // experiment with larger values to get good-enough integer solutions sooner
00605   parm.tm_lim = 500.; // time limit, in milliseconds. 1000 = 1 second. to do: modify this for which set of bounds
00606   int status = glp_intopt(lp->glp, &parm);
00607   if (status == GLP_ETMLIM) // timed out
00608   {
00609     // OK if it is non-optimal, but not if it is non-integer
00610     solved = solution_is_integer() && solution_satisfies_constraints();
00611     return solved;
00612   }
00613   else if (status == 0 || status == GLP_EMIPGAP)
00614   {
00615     optimal = true;
00616     solved = true;
00617     return true;
00618   }
00619   
00620   return false;
00621 }
00622 
00623 
00624 bool IAMilp::solution_satisfies_constraints()
00625 {
00626   bool unsatisfied_found(false);
00627   for (int j = 0; j < data->constraints.size(); ++j)
00628   {
00629     double slack=0.;
00630     const IAData::constraintRow & c = data->constraints[j];
00631     for (std::vector<IAData::sparseEntry>::const_iterator i = c.M.begin(); i < c.M.end(); ++i)
00632     {
00633       const double xv = glp_mip_col_val( lp->glp, x_i(i->col));
00634       slack += xv * i->val;
00635     }
00636     if (data->constraints.front().upperBound == data->constraints.front().lowerBound)
00637     {
00638       if ( fabs(slack - data->constraints.front().upperBound) > constraint_tolerance)
00639       {
00640         unsatisfied_found = true;
00641         if (debugging) 
00642           print_constraint(j);
00643         else
00644           return false;
00645       }
00646     }
00647     else
00648     {
00649       printf(" in [%1.1f,%1.1f]", data->constraints.front().upperBound, data->constraints.front().lowerBound );
00650       if ( ( (data->constraints.front().upperBound - slack) < -constraint_tolerance ) ||
00651            ( (data->constraints.front().lowerBound - slack) > constraint_tolerance ) )
00652       {
00653         unsatisfied_found = true;
00654         if (debugging) 
00655         {
00656           printf("Unsatisfied constraint: ");
00657           print_constraint(j);
00658         }
00659         else
00660           return false;
00661       }
00662     }
00663   }
00664   return !unsatisfied_found;
00665 }
00666 
00667 void IAMilp::print_constraint(int j)
00668 {
00669   //printf("constraint %d: ", j);
00670   double slack=0.;
00671   const IAData::constraintRow & c = data->constraints[j];
00672   for (std::vector<IAData::sparseEntry>::const_iterator i = c.M.begin(); i < c.M.end(); ++i)
00673   {
00674     // nlp solution const double xv = solution->x_solution[i->col];
00675     // const double xv = glp_get_col_prim( lp->glp, x_i(i->col) ); relaxed solution only
00676     const double xv = glp_mip_col_val( lp->glp, x_i(i->col));
00677     slack += xv * i->val;
00678     printf(" %1.0f x%d (%1.3f) ", 
00679            i->val, i->col, xv );
00680   }
00681   if (data->constraints.front().upperBound == data->constraints.front().lowerBound)
00682   {
00683     printf(" = %1.1f ", data->constraints.front().upperBound );  
00684     if ( fabs(slack - data->constraints.front().upperBound) > constraint_tolerance )
00685       printf(" UNSATISFIED ");
00686   }
00687   else
00688   {
00689     printf(" in [%1.1f,%1.1f]", data->constraints.front().upperBound, data->constraints.front().lowerBound );
00690     if ( (data->constraints.front().upperBound - slack) < -constraint_tolerance )
00691       printf(" UNSATISFIED BEYOND UPPERBOUND ");
00692     if ( (data->constraints.front().lowerBound - slack) > constraint_tolerance )
00693       printf(" UNSATISFIED BELOW LOWERBOUND ");
00694     
00695   }
00696   printf(" (%1.1f)\n", slack );
00697 }
00698 
00699 bool IAMilp::solution_is_integer()
00700 {
00701   bool non_integer_found(false);
00702   for (int i=0; i<data->num_variables(); ++i)
00703   {
00704     const double x = glp_mip_col_val(lp->glp, x_i(i));
00705     if ( x < 0. || fabs(x - round(x)) > integrality_tolerance)
00706     {
00707       if (debugging)
00708       {
00709         printf("Noninteger variable: ");
00710         print_solution(i);
00711         non_integer_found = true;
00712       }
00713       else
00714         return false;
00715     }
00716   }
00717   return !non_integer_found;
00718 }
00719 
00720 void IAMilp::print_solution(int i)
00721 {
00722   const double x = glp_mip_col_val(lp->glp, x_i(i));
00723   const double xp = glp_mip_col_val( lp->glp, delta_plus_i(i) );
00724   const double xm = glp_mip_col_val( lp->glp, delta_minus_i(i) );
00725   printf("%d: goal (%1.1f) x_nlp (%1.1f): x (%1.1f) plus (%1.1f) minus (%1.1f)\n", 
00726          x_i(i), data->I[i], solution->x_solution[i], x, xp, xm);
00727 }
00728 
00729 
00730 bool IAMilp::get_solution()
00731 {
00732   // print solution
00733   if (debugging)
00734   {
00735     printf("MILP solution:\n");  
00736     // print delta values
00737     printf("x* and deltas\n");
00738     for (int i=0; i<data->num_variables(); ++i)
00739     {
00740       print_solution(i);
00741     }
00742     printf("constraints\n");
00743     for (int j = 0; j < data->constraints.size(); ++j)
00744     {
00745       print_constraint(j);
00746     }
00747     // todo: print objective value
00748   }
00749   
00750   // collect the solution into the vector
00751   for (int i = 0; i<data->num_variables(); ++i)
00752   {
00753     solution->x_solution[i] = glp_mip_col_val( lp->glp, x_i(i) );
00754   }
00755 
00756   // objective function values
00757   // z = glp_get_obj_val(lp);
00758   return true;
00759 }
00760 
00761 // return the solution
00762 bool IAMilp::solve()
00763 {
00764   bool success;
00765   success = create_problem();
00766   assert(success);
00767 
00768   success = set_natural_constraints();
00769   assert(success);
00770   
00771   success = set_sum_even_constratins();
00772   assert(success);
00773   
00774   success = set_deltas();
00775   assert(success);
00776   
00777   success = weight_deltas_1();
00778   assert(success);
00779   
00780   success = set_maxes();
00781   assert(success);
00782 
00783   success = set_objectives_1();
00784   assert(success); 
00785 
00786   // cycle through bounds
00787   // accept any integer solution, even if suboptimal
00788   for (int c=0; c<6; ++c)
00789   {
00790     // preference order: 2, A, 1, 3, 4, B
00791     switch (c)
00792     {
00793       case 0: 
00794         success = set_bounds_1(); 
00795         break;
00796       case 1:
00797         success = set_bounds_2(); 
00798         break;
00799       case 2:
00800         success = set_bounds_A();
00801         break;
00802       case 3: 
00803         success = set_bounds_3(); 
00804         break;
00805       case 4: 
00806         success = set_bounds_4(); 
00807         break;
00808       case 5: 
00809         success = set_bounds_B(); 
00810         break;
00811     }
00812     assert(success);
00813     
00814     // to do: cycle through solving sum of deltas, then add max of deltas or other heuristics to improve the solution
00815     // to do: grant longer time limit for set_bounds_B
00816   
00817     bool optimal(false);
00818     success = glpk_solve(optimal);
00819     if (success && !optimal)
00820     {
00821       if (debugging)
00822         printf("Stopping at sub-optimal but integer MILP solution.\n");
00823     }
00824     if (success)
00825       break;
00826   }
00827   
00828   success = get_solution();
00829   assert(success);
00830   
00831   destroy_problem();
00832   return true;
00833 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines