Hierarchical Index

00001 #define USE_MONOMIALS
00002 #define USE_REFORMULATION
00003 #define CYCLIC
00004 
00005 #include<scil/scil.h>
00006 #include<scil/constraints/tour.h>
00007 
00008 #include<scil/core/monomial.h>
00009 #include<scil/core/polynomial.h>
00010 #include<vector>
00011 
00012 #include<boost/graph/iteration_macros.hpp>
00013 
00014 #include <stdlib.h>
00015 
00016 using namespace SCIL;
00017 
00018    typedef boost::adjacency_list<vecS, vecS, undirectedS> Graph;
00019    typedef boost::graph_traits<Graph> GraphTraits;
00020    typedef GraphTraits::vertex_descriptor vertex_descriptor;
00021    typedef GraphTraits::edge_descriptor edge_descriptor;
00022 
00023 class reformulation_constraint : public sym_constraint{
00024   private:
00025      std::list<cons_obj*> lc;
00026   public:
00027      reformulation_constraint(ILP_Problem &IP, Graph &G, map<vertex_descriptor, list<int> > &needed_tools, var_map<edge_descriptor> &tour_vars, map<vertex_descriptor, var_map<int> > &tool_vars, int capacity, bool cyclic, vertex_descriptor &home, int num_tools){
00028         row r,q;
00029         polynomial p;
00030         int i;
00031         vertex_descriptor v;
00032         edge_descriptor e;
00033         double rhs;
00034         BGL_FORALL_VERTICES(v,G, Graph) if(cyclic || v!=home) {
00035            r=0;
00036            q=0;
00037            BOOST_FOREACH(i, needed_tools[v]) q+=tool_vars[v][i];
00038            for(i=0; i<num_tools; i++) r+=tool_vars[v][i];
00039            r-=q;
00040            rhs = (capacity-(needed_tools[v].size()));
00041            BGL_FORALL_EDGES(e,G, Graph){
00042               if(source(e, G) == v || target(e, G) == v){
00043                  p = (r - rhs) * tour_vars[e];
00044                  std::list<monomial> mons = p.get_p();
00045                  row pl;
00046                  std::list<monomial>::const_iterator mo;
00047                  foreach(mo, mons){
00048                     pl += IP.add_polynomial(*mo, false);
00049                  }
00050                  lc.push_back(pl<=0);
00051               }
00052            }
00053         }
00054      }
00055 
00056      ~reformulation_constraint(){
00057         for(std::list<cons_obj*>::iterator it = lc.begin(); it != lc.end(); it++)
00058            delete *it;
00059      }
00060 
00061      virtual void init(subproblem& S){
00062         return;
00063      }
00064 
00065      status feasible(solution& S){
00066         return feasible_solution;
00067      }
00068 
00069      status standard_separation(subproblem& S) {
00070         std::list<cons_obj*>::iterator c;
00071         cons_obj* cons;
00072         int num = 0;
00073         foreach(c, lc) {
00074            if( (*c)->violation(S) > 1.0e-4 ) {
00075               row r;
00076               (*c)->non_zero_entries(r);
00077               cons = r<=0;
00078               cons->set_qrStatus(NONE);
00079               S.add_basic_constraint(cons);
00080               num++;
00081            }
00082         }
00083         if( num ) {
00084            return constraint_found;
00085         }
00086 
00087         return no_constraint_found;
00088      }
00089 };
00090 
00091 int main(int argc, char** argv) {
00092 
00093   edge_descriptor e;
00094   vertex_descriptor v;
00095   int i,j;
00096   row r,q;
00097   int num_jobs,num_tools,capacity;
00098   int f;
00099   char line[500];
00100   double rhs;
00101   polynomial p,p1;
00102 
00103   /* read instance */
00104 
00105   scanf("n=%d\n",&num_jobs);
00106   scanf("m=%d\n",&num_tools);
00107   scanf("min=%d\n",&f);
00108   scanf("max=%d\n",&f);
00109   scanf("c=%d\n",&capacity);
00110   for( int i = 1; i <= atoi(argv[1]); i++ ) {
00111   do scanf("%500s %d:\n",line,&f);
00112   while(strcmp(line,"problem"));
00113   }
00114   printf("======================================================================\n");
00115   printf("Solving problem %d (%d jobs, %d tools, capacity %d)\n",
00116          f,num_jobs,num_tools,capacity);
00117   printf("======================================================================\n");
00118   scanf("%500s\n",line);
00119   scanf("\n");
00120 
00121   Graph G;
00122   vector<vertex_descriptor> ref(num_jobs+1);
00123   for(i=0;i<num_jobs;i++) ref[i]=add_vertex(G);
00124   for(i=0;i<num_jobs;i++) for(j=i+1;j<num_jobs;j++)
00125     add_edge(ref[i],ref[j], G);
00126   vertex_descriptor home = ref[0];
00127   bool cyclic = true;
00128 
00129   map<vertex_descriptor, list<int> > needed_tools;
00130   for(i=0;i<num_tools;i++) for(j=0;j<num_jobs;j++) {
00131     scanf("%d",&f);
00132     if(f) needed_tools[ref[j]].push_back(i);
00133   };
00134 
00135 #ifndef CYCLIC
00136   cyclic = false;
00137   ref[num_jobs++]=home=add_vertex(G);
00138   BGL_FORALL_VERTICES(v,G, Graph) if(v!=home) add_edge(home,v, G);
00139 #endif
00140 
00141   /* set up ilp */
00142 
00143   ILP_Problem IP(Optsense_Min);
00144 
00145   // tour variables x_ij
00146   var_map<edge_descriptor> tour_vars;
00147   BGL_FORALL_EDGES(e,G, Graph) {
00148     tour_vars[e]=IP.add_binary_variable(0);
00149   }
00150 
00151   // tool variables y_it
00152   map<vertex_descriptor, var_map<int> > tool_vars;
00153   BGL_FORALL_VERTICES(v,G, Graph) if(cyclic || v!=home) for(i=0; i<num_tools; i++) {
00154     tool_vars[v][i]=IP.add_binary_variable(0);
00155   };
00156 
00157   // tour constraints
00158   IP.add_sym_constraint(new TOUR<Graph>(G,tour_vars));
00159 
00160   // capacity constraints
00161   std::list<cons_obj*> lc;
00162   BGL_FORALL_VERTICES(v,G, Graph) if(cyclic || v!=home) {
00163      r=0;
00164      q=0;
00165      BOOST_FOREACH(i, needed_tools[v]) q+=tool_vars[v][i];
00166      for(i=0; i<num_tools; i++) r+=tool_vars[v][i];
00167      r-=q;
00168      rhs = (capacity-(needed_tools[v].size()));
00169      IP.add_basic_constraint(r<=rhs,Static);
00170   }
00171 
00172 #ifdef USE_REFORMULATION
00173   //adds a symbolic constraint for reformulation of capacity constraints
00174   IP.add_sym_constraint(new reformulation_constraint(IP,G,needed_tools,tour_vars,tool_vars,capacity,cyclic,home,num_tools));
00175 #endif
00176 
00177   // tool constraints
00178   BGL_FORALL_VERTICES(v,G, Graph) if(cyclic || v!=home) BOOST_FOREACH(i,needed_tools[v]) {
00179     r=tool_vars[v][i];
00180     IP.add_basic_constraint(r==1,Static);
00181   }
00182 
00183 #ifdef USE_MONOMIALS
00184 
00185   // monomials in objective function
00186   BGL_FORALL_EDGES(e,G, Graph) {
00187     for(i=0; i<num_tools; i++) {
00188       if(cyclic || (source(e, G)!=home && target(e, G)!=home))
00189         IP.add_polynomial(-1 * tour_vars[e] * tool_vars[source(e, G)][i] * tool_vars[target(e, G)][i]);
00190       if(cyclic || source(e, G)!=home)
00191         IP.add_polynomial(0.5 * tour_vars[e] * tool_vars[source(e, G)][i]);
00192       if(cyclic || target(e, G)!=home)
00193         IP.add_polynomial(0.5 * tour_vars[e] * tool_vars[target(e, G)][i]);
00194     }
00195   };
00196 
00197 #else
00198 
00199   // linearization variables z_ijt
00200   map<edge_descriptor, var_map<int> > lin_vars;
00201   BGL_FORALL_EDGES(e,G, Graph) for(i=0; i<num_tools; i++) {
00202     lin_vars[e][i]=IP.add_binary_variable(0.5);
00203   }
00204 
00205 
00206   // linearization constraints
00207   BGL_FORALL_EDGES(e,G, Graph) {
00208     for(i=0; i<num_tools; i++) {
00209       r=lin_vars[e][i];
00210       r-=tour_vars[e];
00211       if(cyclic || source(e, G)!=home) r+=tool_vars[source(e, G)][i];
00212       if(cyclic || target(e, G)!=home) r-=tool_vars[target(e, G)][i];
00213       IP.add_basic_constraint(r>=-1.0,Static);
00214       r=lin_vars[e][i];
00215       r-=tour_vars[e];
00216       if(cyclic || source(e, G)!=home) r-=tool_vars[source(e, G)][i];
00217       if(cyclic || target(e, G)!=home) r+=tool_vars[target(e, G)][i];
00218       IP.add_basic_constraint(r>=-1.0,Static);
00219     }
00220   }
00221 
00222 #endif
00223 
00224   /* solve ilp */
00225 
00226   // optimize
00227   IP.optimize();
00228 
00229   // determine order
00230   map<vertex_descriptor, int> inv_ref;
00231   for(i=0;i<num_jobs;i++) inv_ref[ref[i]]=i;
00232   int *s1=new int[num_jobs];
00233   int *s2=new int[num_jobs];
00234   for(i=0;i<num_jobs;i++) s1[i]=s2[i]=-1;
00235   BGL_FORALL_EDGES(e,G, Graph) if(IP.get_solution(tour_vars[e])) {
00236     if(s1[inv_ref[source(e, G)]]==-1) s1[inv_ref[source(e, G)]] = inv_ref[target(e, G)];
00237     else s2[inv_ref[source(e, G)]] = inv_ref[target(e, G)];
00238     if(s1[inv_ref[target(e, G)]]==-1) s1[inv_ref[target(e, G)]] = inv_ref[source(e, G)];
00239     else s2[inv_ref[target(e, G)]] = inv_ref[source(e, G)];
00240   }
00241   list<int> order;
00242   j=num_jobs-1;
00243   order.push_back(j);
00244   while(order.size()<num_jobs) {
00245     i=j;
00246     if(s1[i]!=-1) j=s1[i];
00247     else j=s2[i];
00248     if(s1[j]==i) s1[j]=-1;
00249     else s2[j]=-1;
00250     order.push_back(j);
00251   }
00252   delete[] s1;
00253   delete[] s2;
00254 
00255   // print solution
00256   printf("======================================================================\n");
00257   printf("optimal solution:\n");
00258   int last=order.back();
00259   int c1,c2;
00260   c2=0;
00261   BOOST_FOREACH(i,order) {
00262     if(cyclic || ref[i]!=home) {
00263       printf("  job %2d, tools",i+1);
00264       for(j=0;j<num_tools;j++)
00265         if(IP.get_solution(tool_vars[ref[i]][j]))
00266           printf(" %2d",j+1);
00267       c1=0;
00268       if(cyclic || ref[last]!=home) {
00269         for(j=0;j<num_tools;j++)
00270           if(IP.get_solution(tool_vars[ref[i]][j])
00271              && !IP.get_solution(tool_vars[ref[last]][j]))
00272             c1++;
00273         printf(", %2d setup(s)",c1);
00274       }
00275       else {
00276         for(j=0;j<num_tools;j++)
00277           if(IP.get_solution(tool_vars[ref[i]][j]))
00278             c1++;
00279         printf(", %2d setup(s)",c1);
00280       }
00281       c2+=c1;
00282       printf("\n");
00283     }
00284     last=i;
00285   }
00286   printf("total number of setups: %2d\n",c2);
00287   printf("======================================================================\n");
00288 
00289   return 0;
00290 }
00291
ToolSwitching.cc
