production.c

#line 2 "production.lw"
#include<scil/scil.h>
#include<scil/constraints/flow.h>
//#include<scil/constraints/cp_knapsack.h>

using namespace SCIL;
using namespace LEDA;

#define NI 5
#define NF 10
#define NC 5

class production_heuristic : public primal_heuristic {
  public:

  int heuristic(subproblem& S) {
    cout<<"hier\n";
    return 0;
  }
};

int main() {

  // Define Input

  int W[NI]={ 7,8,9,10,12 };
  int Cap[NF]={ 40, 45, 49, 53, 57, 63, 65, 71, 71, 73 };
  int D[NI][NC] = { { 1, 0, 0, 2, 1 },
                    { 2, 2, 2, 0, 0 },
                    { 0, 1, 1, 2, 1 },
                    { 0, 1, 0, 0, 1 },
                    { 1, 1, 0, 0, 2 } };
  double OC[NF];

  node CN[NC];
  node FN[NF];

  double x;
  random_source RS;
  RS.set_seed(99);
  for(int i=0; i<NF; i++) {
    RS>>x;
    OC[i]=100*x;
  }

  graph G;
  for(int i=0; i<NF; i++) FN[i]=G.new_node();
  for(int i=0; i<NC; i++) CN[i]=G.new_node();

  for(int i=0; i<NF; i++) for(int j=0; j<NC; j++) {
    G.new_edge(FN[i], CN[j]);
  }

  edge e;
  edge_array<double> C(G);
  forall_edges(e,G) {
    RS>>x;
    C[e]=1*x;
  }

  // Define Model

  ILP_Problem IP(Optsense_Min);

  array<row_map<node> > NVM(NI); // TODO init
  row_map<node> y;
  var v;

  node u;

  for(int i=0; i<NF; i++) {
    u=FN[i];
    for(int j=0; j<NI; j++) {
      v=IP.add_variable(0, 0, 10, Vartype_Integer);
      NVM[j][u]=v;
    }
    v=IP.add_variable(OC[i], 0, 1, Vartype_Integer);
    y[u]=v;

    row r;
    for(int j=0; j<NI; j++) {
      r+=W[j]*NVM[j][u]; };
    IP.add_basic_constraint(r <= Cap[i]*y[u]);
    //IP.add_sym_constraint(new cp_knapsack(r <= Cap[i]*y[u]));
  }

  for(int i=0; i<NC; i++) {
    for(int j=0; j<NI; j++) {
      NVM[j][CN[i]]=D[i][j];
    }
  }

  row_map<edge> EVM[NI]; // TODO init
  for(int i=0; i<NI; i++) {
    forall_edges(e, G) {
      v=IP.add_variable(C[e], 0, 100, Vartype_Float);
      EVM[i][e]=v;
    }
    IP.add_sym_constraint(new FLOW(G, NVM[i], EVM[i]));
  }

  // Solve

  IP.set_primal_heuristic(new production_heuristic());

  IP.optimize();

  // Display Solution

  for(int i=0; i<NF; i++) {
    cout<<IP.get_solution(y[FN[i]])<<" solution of facility "<<i<<endl;
    if(IP.get_solution(y[FN[i]])>0.5) {
      cout<<"It produces\n";
      for(int j=0; j<NI; j++) {
        cout<<IP.get_solution(NVM[j][FN[i]])<<" Items of item "<<j<<"\n";
      }
    }
  }

  return 0;
}

---