computational models bulk commit

Specialization/LP_Problems/Maintenance_Problem.cpp → Specialization/LP_Problems/Maintenance_Problem/Maintenance_Problem.cpp

 #include "Maintenance_Problem.h"
 
-Maintenance_Problem::Maintenance_Problem() : Linear_Program(true), _g({{"Flow", Attribute(max, 0)},{"Upper", Attribute(fix, 1)}, {"Selected", Attribute(fix, 0)}}) {
+Maintenance_Problem::Maintenance_Problem() : Linear_Program(true), _number_of_epochs(0), _g({{"Flow", Attribute(max, 0)},{"Upper", Attribute(fix, 1)}, {"Selected", Attribute(fix, 0)}}) {
 
 }
 
-Maintenance_Problem::Maintenance_Problem(Maintenance_Problem& mp) : Linear_Program(mp), _g(mp._g) {
+Maintenance_Problem::Maintenance_Problem(Maintenance_Problem& mp) : Linear_Program(mp), _number_of_epochs(mp._number_of_epochs), _g(mp._g) {
   mp._g = Graph();
 }
 
-Maintenance_Problem::Maintenance_Problem(Maintenance_Problem&& mp) : Linear_Program(std::move(mp)), _g(std::move(mp._g)) {
+Maintenance_Problem::Maintenance_Problem(Maintenance_Problem&& mp) : Linear_Program(std::move(mp)), _number_of_epochs(std::move(mp._number_of_epochs)), _g(std::move(mp._g)) {
   mp._g = Graph();
 }
 
 Maintenance_Problem::Maintenance_Problem(
   Graph& g, Node* source,
   Node* target,
-  size_t intervals
+  size_t epochs
 )
-  : Linear_Program(true), _g(g)
+  : Linear_Program(true), _number_of_epochs(epochs), _g(g)
 {
   std::set<Edge*> critical = g.gather_edges({"Selected"});
 
@@ -45,7 +45,7 @@ Maintenance_Problem::Maintenance_Problem(
   /*
     flow conservation constraints for all epochs
   */
-  for(size_t epoch = 0; epoch < intervals; ++epoch){
+  for(size_t epoch = 0; epoch < epochs; ++epoch){
     std::stringstream name_appendix;
     name_appendix << "epoch_" << epoch;
     lookup.push_back(lp_generator::grow_from_graph(
@@ -86,7 +86,7 @@ Maintenance_Problem::Maintenance_Problem(
     critical edges processed
   */
   for(Edge* e : critical){
-    std::map<Variable*, Coefficient> lhs;
+    std::unordered_map<Variable*, Coefficient> lhs;
     for(auto lookup_during_epoch : lookup){
       lhs.insert({lookup_during_epoch.second.find({e, "Selected"})->second.first, {1}});
     }
@@ -105,11 +105,11 @@ Maintenance_Problem::Maintenance_Problem(
   /*
     value bounded by flow
   */
-  std::pair<Variable*, size_t> target_var = this->polyeder().add_variable(Variable("target_variable", Integral));
+  std::pair<Variable*, size_t> target_var = this->polyeder().add_variable(Variable("target_variable", Continuous));
 
-  for(size_t epoch = 0; epoch < intervals; ++epoch){
+  for(size_t epoch = 0; epoch < epochs; ++epoch){
     auto lookup_during_epoch = lookup[epoch];
-    std::map<Variable*, Coefficient> lhs;
+    std::unordered_map<Variable*, Coefficient> lhs;
     for(Edge* e : source->incident()){
       lhs.insert({lookup_during_epoch.second.find({e, "Flow"})->second.first, {-1}});
     }
@@ -134,14 +134,50 @@ Graph& Maintenance_Problem::network(){
   return this->_g;
 }
 
+size_t Maintenance_Problem::number_of_epochs(){
+  return this->_number_of_epochs;
+}
+
+std::vector<SCIP*> Maintenance_Problem::all_computational_models(){
+  // basic model:
+  SCIP* scip_basic = Linear_Program::computational_model();
+
+  // semiassignment_branching
+  SCIP* scip_semiassign = Linear_Program::computational_model();
+  SCIP_CONS** conss = SCIPgetOrigConss(scip_semiassign);
+  std::vector<SCIP_CONS*> process_cons;
+  for(int cons_index = 0; cons_index < SCIPgetNOrigConss(scip_semiassign); cons_index++){
+    if(std::string(conss[cons_index]->name).find("processed") != std::string::npos){
+      process_cons.push_back(conss[cons_index]);
+    }
+  }
+  SCIPincludeConshdlrSemiassign(scip_semiassign);
+
+  SCIP_CONS** process_cons_array = new SCIP_CONS*[process_cons.size()];
+  std::copy(process_cons.begin(), process_cons.end(), process_cons_array);
+  SCIPincludeBranchruleSemiassign(scip_semiassign, process_cons.size(), process_cons_array);
+
+  // basic benders feasibility cuts:
+  SCIP* scip_benders_feas = this->nmp_generate_benders();
+
+  // benders feas with basic cut inequalities
+  SCIP* scip_benders_feas_basic_cut = this->nmp_generate_benders();
+  SCIP_BENDERS* benders_basic_cut = SCIPfindBenders(scip_benders_feas_basic_cut, "default");
+  SCIP_CALL_ABORT( SCIPincludeBenderscutBasicSepa(scip_benders_feas_basic_cut, benders_basic_cut));
+
+  return {scip_basic, scip_semiassign, scip_benders_feas, scip_benders_feas_basic_cut};
+}
+
 void Maintenance_Problem::operator=(Maintenance_Problem& mp){
   this->Linear_Program::operator=(mp);
+  this->_number_of_epochs = mp._number_of_epochs;
   this->_g = mp._g;
   mp._g = Graph();
 }
 
 void Maintenance_Problem::operator=(Maintenance_Problem&& mp){
   this->Linear_Program::operator=(mp);
+  this->_number_of_epochs = std::move(mp._number_of_epochs);
   this->_g = std::move(mp._g);
   mp._g = Graph();
 }

Specialization/LP_Problems/Maintenance_Problem.h → Specialization/LP_Problems/Maintenance_Problem/Maintenance_Problem.h

@@ -5,12 +5,27 @@
 #include <vector>
 #include <sstream>
 
-#include "../../Graphtheory/Graphtheory.h"
-#include "../../Linear_Programming/Linear_Program.h"
-#include "../../Linear_Programming/lp_generator.h"
+#include "../../../Graphtheory/Graphtheory.h"
+#include "../../../Linear_Programming/Linear_Program.h"
+#include "../../../Linear_Programming/lp_generator.h"
+
+// computational
+#include "computational_types/semiassignment_hdlr.h"
+#include "computational_types/semiassignment_branch.h"
+#include "computational_types/bulksepa/nmp_bulksepa.h"
+
+#include <scip/scip.h>
+#include <scip/scipshell.h>
+#include <scip/scipdefplugins.h>
+#include <scip/cons_benders.h>
+#include <scip/cons_benderslp.h>
+#include <scip/benders_default.h>
 
 class Maintenance_Problem : public Linear_Program {
+  size_t _number_of_epochs;
   Graph _g;
+
+  SCIP* nmp_generate_benders();
 public:
   Maintenance_Problem();
   Maintenance_Problem(Maintenance_Problem& mp);
@@ -18,6 +33,9 @@ public:
   Maintenance_Problem(Graph& g, Node* source, Node* target, size_t intervals);
 
   Graph& network();
+  size_t number_of_epochs();
+
+  std::vector<SCIP*> all_computational_models();
 
   void operator=(Maintenance_Problem& mp);
   void operator=(Maintenance_Problem&& mp);

Specialization/LP_Problems/Maintenance_Problem/computational_types/bulksepa/nmp_bulksepa.cpp

+#include "nmp_bulksepa.h"
+
+struct nmp_bulksepa_container{
+  SCIP* bulksepa_scip;
+  SCIP_VAR**  coefficient_vars;
+  int nedges;
+  int nepochs;
+};
+
+
+struct SCIP_VarData{
+  int nepochs;
+  SCIP_VAR** vars_of_repr_edge;
+};
+
+struct SCIP_BenderscutData
+{
+  SCIP* master;
+  SCIP_BENDERS* benders;
+  std::vector<SCIP_CONS*>* capture_inequalities;
+};
+
+SCIP_RETCODE create_bulksepa_scip(SCIP* master, SCIP_SOL* solution, SCIP_Real* min_cut_capacities, SCIP_VAR** mp_cut_vars, SCIP_VAR* target_var, int nepochs, int nedges, nmp_bulksepa_container* container){
+  SCIP* nmp_bulksepa;
+  SCIPcreate(&nmp_bulksepa);
+  char name[SCIP_MAXSTRLEN];
+  SCIP_Real bound = SCIPgetDualbound(master);
+
+  SCIP_CALL( SCIPincludeDefaultPlugins(nmp_bulksepa));
+  SCIP_CALL( SCIPcreateProbBasic(nmp_bulksepa, "ephemeral_bulksepa"));
+  SCIP_CALL( SCIPsetObjsense(nmp_bulksepa, SCIP_OBJSENSE_MINIMIZE));
+  SCIP_CALL( SCIPsetMessagehdlr(nmp_bulksepa, NULL));
+
+  // create vars
+  SCIP_VAR**  coefficient_vars;
+  SCIP_CALL( SCIPallocBlockMemoryArray(nmp_bulksepa, &coefficient_vars, nepochs*nedges));
+  for(int epoch = 0; epoch < nepochs; epoch++){
+    for(int edge = 0; edge < nedges; edge++){
+      SCIPsnprintf(name, SCIP_MAXSTRLEN, "coefficient_of_%s", SCIPvarGetName(mp_cut_vars[edge*nepochs+epoch]));
+      SCIP_CALL( SCIPcreateVarBasic(nmp_bulksepa, &coefficient_vars[edge*nepochs+epoch], name, 0, bound, SCIPgetSolVal(master, solution, mp_cut_vars[edge*nepochs+epoch]), SCIP_VARTYPE_CONTINUOUS));
+      SCIP_CALL( SCIPaddVar(nmp_bulksepa, coefficient_vars[edge*nepochs+epoch]));
+    }
+  }
+
+  SCIP_CALL( SCIPsetPresolving(nmp_bulksepa, SCIP_PARAMSETTING_OFF, TRUE));
+
+  // add cut_sepa_sepa_conshdlr
+  SCIP_CALL( SCIPincludeConshdlrBulkSepa(nmp_bulksepa, coefficient_vars, min_cut_capacities, nedges, nepochs, bound, FALSE));
+  container->bulksepa_scip = nmp_bulksepa;
+  container->coefficient_vars = coefficient_vars;
+  container->nedges = nedges;
+  container->nepochs = nepochs;
+
+  // TODO free
+  return SCIP_OKAY;
+}
+
+static
+SCIP_RETCODE exec_sepa(SCIP* scip, SCIP_BENDERS* benders, SCIP_BENDERSCUT* benderscut, SCIP_SOL* sol, int probnumber, SCIP_BENDERSENFOTYPE type, SCIP_RESULT* result)
+{
+  assert(benders != NULL);
+  SCIP* subproblem = SCIPbendersSubproblem(benders, probnumber);
+
+  // finding target_var
+  SCIP_VAR* target_var_subproblem = SCIPfindVar(subproblem, "target_variable");
+  SCIP_VAR* target_var;
+  SCIPgetBendersMasterVar(scip, benders, target_var_subproblem, &target_var);
+  assert(target_var != NULL);
+
+  if(SCIPisZero(scip, SCIPgetVarSol(scip, target_var)) || !SCIPisLPDualReliable(scip)){
+    *result = SCIP_DIDNOTRUN;
+    return SCIP_OKAY;
+  }
+
+
+  // gathering cons with dual val "!= 0" which form a min cut
+  SCIP_CONS** conss = SCIPgetConss(subproblem);
+  std::vector<SCIP_Real> min_cut_capacities;
+  std::vector<SCIP_CONS*> min_cut_conss;
+  for(int cons_index = 0; cons_index < SCIPgetNConss(subproblem); cons_index++){
+    if(std::string(SCIPconsGetName(conss[cons_index])).find("capacity") != std::string::npos){
+      SCIP_Real dual_val;
+      SCIP_Bool result;
+      SCIPconsGetDualfarkas(subproblem, conss[cons_index], &dual_val, &result);
+      if(!SCIPisFeasZero(subproblem, dual_val)){
+        SCIP_Bool rhs_result = FALSE;
+        min_cut_capacities.push_back(SCIPconsGetRhs(subproblem, conss[cons_index], &rhs_result));
+        assert(rhs_result);
+        min_cut_conss.push_back(conss[cons_index]);
+      }
+    }
+  }
+
+  // extracting the corresponding variables
+  std::vector<SCIP_VAR*> subproblem_cut_vars;
+  subproblem_cut_vars.reserve(min_cut_conss.size());
+  for(SCIP_CONS* cons : min_cut_conss){
+    SCIP_VAR* cons_vars[2];
+    SCIP_Bool result = FALSE;
+    SCIP_CALL( SCIPgetConsVars(subproblem, cons, cons_vars, 2, &result));
+    if(result == false){
+      continue;
+    }
+
+    for(int i = 0; i < 2; i++){
+      if(std::string(SCIPvarGetName(cons_vars[i])).find("Selected") != std::string::npos){
+        subproblem_cut_vars.push_back(cons_vars[i]);
+        continue;
+      }
+    }
+  }
+  int nedges = min_cut_capacities.size();
+  int nepochs = SCIPbendersGetNSubproblems(benders);
+
+  // finding mp decision variables of all epochs
+  SCIP_VAR** mp_cut_vars;
+  SCIP_CALL( SCIPallocBufferArray(scip, &mp_cut_vars, nedges*nepochs));
+  for(int edge = 0; edge < nedges; edge++){
+    SCIP_VAR* mp_var;
+    SCIP_CALL( SCIPgetBendersMasterVar(scip, benders, subproblem_cut_vars[edge], &mp_var));
+    assert(mp_var != NULL);
+    SCIP_VarData* vardata = SCIPvarGetData(mp_var);
+    assert(vardata != NULL);
+
+    for(int epoch = 0; epoch < nepochs; epoch++){
+      mp_cut_vars[edge*nepochs+epoch] = vardata->vars_of_repr_edge[epoch]; // TODO: bulkcpy possible?
+    }
+  }
+
+  nmp_bulksepa_container container;
+  SCIP_CALL( create_bulksepa_scip(scip, sol, min_cut_capacities.data(), mp_cut_vars, target_var, nepochs, nedges, &container));
+
+  SCIPsolve(container.bulksepa_scip);
+  SCIP_SOL* opt_sol = SCIPgetBestSol(container.bulksepa_scip);
+
+
+  // checking if problem separated solution; if yes add contraint
+  //std::cout << SCIPgetVarSol(scip, target_var) << ">" << SCIPgetSolOrigObj(container.bulksepa_scip, opt_sol) << "? if so add capture-ineq" << std::endl;
+  if(SCIPisGT(scip, SCIPgetVarSol(scip, target_var), SCIPgetSolOrigObj(container.bulksepa_scip, opt_sol))){
+    std::cout << "capture inequality found! " << SCIPgetNConss(container.bulksepa_scip) << " cons had to be separated." << std::endl;
+    char name[SCIP_MAXSTRLEN];
+
+    // extracting optimal solution from container.bulksepa_scip
+    SCIP_Real* best_coeffs;
+    SCIP_CALL( SCIPallocBufferArray(scip, &best_coeffs, nedges*nepochs));
+    for(int edge = 0; edge < nedges; edge++){
+      for(int epoch = 0; epoch < nepochs; epoch++){
+        best_coeffs[edge*nepochs+epoch] = SCIPgetSolVal(container.bulksepa_scip, opt_sol, container.coefficient_vars[edge*nepochs+epoch]);
+      }
+    }
+
+    SCIP_CONS* cons;
+    SCIP_BENDERSCUTDATA* benderscut_data = SCIPbenderscutGetData(benderscut);
+    assert(benderscut_data != NULL);
+    SCIPsnprintf(name, SCIP_MAXSTRLEN, "capture_inequality_%d", benderscut_data->capture_inequalities->size());
+    SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, nedges*nepochs, mp_cut_vars, best_coeffs, 0, SCIPinfinity(scip))); // TODO: cut? SCIPstoreBenderscutCons function does not exist?
+    SCIP_CALL( SCIPaddCoefLinear(scip, cons, target_var, -1));
+    SCIP_CALL( SCIPaddCons(scip, cons));
+    SCIP_CALL( SCIPreleaseCons(scip, &cons));
+
+    benderscut_data->capture_inequalities->push_back(cons);
+
+    // release buffer
+    SCIPfreeBufferArray(scip, &best_coeffs);
+    SCIPfreeBufferArray(scip, &mp_cut_vars);
+
+    *result = SCIP_CONSADDED;
+  }else{
+    //std::cout << "no admissible capture inequality found" << std::endl;
+    *result = SCIP_DIDNOTFIND;
+  }
+
+
+  for(int index = 0; index < nepochs*nedges; index++){
+    SCIPreleaseVar(container.bulksepa_scip, &container.coefficient_vars[index]);
+  }
+  SCIPfreeBlockMemoryArray(container.bulksepa_scip, &container.coefficient_vars, nedges*nepochs);
+  //SCIPfree(&container.bulksepa_scip);
+
+  return SCIP_OKAY;
+}
+
+SCIP_RETCODE free_benderscut(SCIP* scip, SCIP_BENDERSCUT* benderscut){
+  SCIP_BENDERSCUTDATA* data;
+  data = SCIPbenderscutGetData(benderscut);
+  assert(data != NULL);
+
+  delete data->capture_inequalities;
+
+  SCIPfreeBlockMemory(scip, &data);
+  SCIPbenderscutSetData(benderscut, NULL);
+  std::cin.ignore();
+
+  return SCIP_OKAY;
+}
+
+SCIP_RETCODE SCIPincludeBenderscutBasicSepa(
+   SCIP*                 scip,
+   SCIP_BENDERS*         benders
+   )
+{
+  SCIP_BENDERSCUTDATA* benderscutdata;
+  SCIP_BENDERSCUT* benderscut;
+
+  assert(benders != NULL);
+
+  SCIP_CALL( SCIPallocBlockMemory(scip, &benderscutdata) );
+  benderscutdata->master = scip;
+  benderscutdata->benders = benders;
+  benderscutdata->capture_inequalities = new std::vector<SCIP_CONS*>();
+
+  SCIP_CALL( SCIPincludeBenderscutBasic(scip, benders, &benderscut, "cut capture", "uses min cut from subproblem to separate facets of the corresponding maintenance problem contrained to cut network",
+       -100000, TRUE, exec_sepa, benderscutdata) );
+
+  SCIP_CALL( SCIPsetBenderscutFree(scip, benderscut, &free_benderscut));
+
+  assert(benderscut != NULL);
+  return SCIP_OKAY;
+}

Specialization/LP_Problems/Maintenance_Problem/computational_types/bulksepa/nmp_bulksepa.h

+#ifndef MAINTENANCE_PROBLEM_BASIC_CUT_SEPA_H
+
+#include <iostream>
+#include <vector>
+
+#include "scip/scip.h"
+#include <scip/scipdefplugins.h>
+
+#include "nmp_bulksepa_conshdlr.h"
+
+extern "C"{
+  SCIP_RETCODE SCIPincludeBenderscutBasicSepa(
+    SCIP*                 scip,
+    SCIP_BENDERS*         benders
+  );
+}
+
+#endif

Specialization/LP_Problems/Maintenance_Problem/computational_types/bulksepa/nmp_bulksepa_conshdlr.cpp

+#include "nmp_bulksepa_conshdlr.h"
+
+//#define LOG
+
+struct SCIP_ConsData
+{
+};
+
+
+struct cons_separator_container
+{
+  SCIP* cons_separator;
+  SCIP_VAR** vars; // [edge*nedges+epoch] for decision [nedges*nepochs] for target variable
+  SCIP_CONS** maintenance_select_conss;
+};
+
+struct SCIP_ConshdlrData
+{
+  SCIP* separator;
+  SCIP_VAR** separator_vars;
+  SCIP_Real* capacities;
+  int nedges;
+  int nepochs;
+  SCIP_Real target_upper_bound;
+  int number_of_added_cons;
+  SCIP_Bool use_reopt;
+  cons_separator_container* cons_separator_container_for_reopt;
+  SCIP_RETCODE (*create_cons_separator_function)(SCIP_ConshdlrData*, SCIP_SOL*, cons_separator_container*);
+};
+
+// helper functions
+SCIP_RETCODE create_cons_separator_benders(SCIP_ConshdlrData* data, SCIP_SOL* solution, cons_separator_container* container){
+  SCIP* separator = data->separator;
+  SCIP_VAR** separator_vars = data->separator_vars;
+  SCIP_Bool use_reopt = data->use_reopt;
+  int nedges = data->nedges;
+  int nepochs = data->nepochs;
+  SCIP_Real target_upper_bound = data->target_upper_bound;
+  SCIP_Real* capacities = data->capacities;
+
+  char name[SCIP_MAXSTRLEN];
+
+  SCIP* master;
+  SCIP_CALL( SCIPcreate(&master));
+  SCIP_CALL( SCIPincludeDefaultPlugins(master));
+  SCIP_CALL( SCIPcreateProbBasic(master, "mp"));
+  SCIP_CALL( SCIPsetObjsense(master, SCIP_OBJSENSE_MINIMIZE));
+  SCIP_CALL( SCIPenableReoptimization(master, use_reopt));
+  SCIP_CALL( SCIPsetMessagehdlr(master, NULL));
+
+
+  SCIP_CALL( SCIPsetIntParam(master, "presolving/maxrestarts",0) );
+  SCIP_CALL( SCIPsetIntParam(master, "heuristics/trysol/freq", 1) );
+
+  SCIP_VAR** mp_vars;
+  SCIP_CALL( SCIPallocBlockMemoryArray(master, &mp_vars, nedges*nepochs+1));
+  for(int edge = 0; edge < nedges; edge++){
+    for(int epoch = 0; epoch < nepochs; epoch++){
+      SCIPsnprintf(name, SCIP_MAXSTRLEN, "Selected_%d_%d", edge, epoch);
+      SCIP_CALL( SCIPcreateVarBasic(master, &mp_vars[edge*nepochs+epoch], name, 0, 1, SCIPgetSolVal(separator, solution, separator_vars[edge*nepochs+epoch]), SCIP_VARTYPE_BINARY));
+      SCIP_CALL( SCIPaddVar(master, mp_vars[edge*nepochs+epoch]));
+    }
+  }
+  SCIP_CALL( SCIPcreateVarBasic(master, &mp_vars[nedges*nepochs], "target_variable", 0, target_upper_bound, 1, SCIP_VARTYPE_CONTINUOUS));
+  SCIP_CALL( SCIPaddVar(master, mp_vars[nedges*nepochs]));
+
+
+  // const
+  SCIP_CONS** maintenance_select_conss;
+  SCIP_CALL( SCIPallocBlockMemoryArray(master, &maintenance_select_conss, nedges));
+  for(int edge = 0; edge < nedges; edge++){
+    SCIPsnprintf(name, SCIP_MAXSTRLEN, "maintenance_epoch_select_%d", edge);
+    SCIP_CALL( SCIPcreateConsBasicLinear(master, &maintenance_select_conss[edge], name, 0, NULL, NULL, 1, 1));
+    for( int epoch = 0; epoch < nepochs; epoch++){
+      SCIP_CALL( SCIPaddCoefLinear( master, maintenance_select_conss[edge], mp_vars[edge*nepochs+epoch], 1));
+    }
+    SCIP_CALL( SCIPaddCons(master, maintenance_select_conss[edge]));
+  }
+
+  // subproblems
+  SCIP** subproblems;
+  SCIP_CALL( SCIPallocBufferArray(master, &subproblems, nepochs));
+
+  SCIP_VAR* flow_var;
+  SCIP_VAR* decision_var;
+  SCIP_CONS* flow_cons;
+  SCIP_CONS* capacity_cons;
+  for(int epoch = 0; epoch < nepochs; epoch++){
+    SCIP_CALL( SCIPcreate(&subproblems[epoch]));
+    SCIP_CALL( SCIPincludeDefaultPlugins(subproblems[epoch]));
+    SCIPsnprintf( name, SCIP_MAXSTRLEN, "sub_%d", epoch);
+    SCIP_CALL( SCIPcreateProbBasic(subproblems[epoch], name));
+    SCIP_CALL( SCIPsetObjsense(subproblems[epoch], SCIP_OBJSENSE_MINIMIZE));
+
+
+    SCIPsnprintf( name, SCIP_MAXSTRLEN, "aggregate_flow_val_cons");
+
+    SCIP_CALL( SCIPcreateConsBasicLinear(subproblems[epoch], &flow_cons, name, 0, NULL, NULL, -SCIPinfinity(subproblems[epoch]), 0));
+
+    //flow vars
+    for(int edge = 0; edge < nedges; edge++){
+      SCIPsnprintf( name, SCIP_MAXSTRLEN, "flow_edge_%d", edge);
+
+      SCIP_CALL( SCIPcreateVarBasic(subproblems[epoch], &flow_var, name, 0, capacities[edge], 0, SCIP_VARTYPE_CONTINUOUS));
+      SCIP_CALL( SCIPaddVar(subproblems[epoch], flow_var));
+
+      SCIP_CALL( SCIPaddCoefLinear(subproblems[epoch], flow_cons, flow_var, -1));
+
+      SCIPsnprintf( name, SCIP_MAXSTRLEN, "Selected_%d_%d", edge, epoch);
+
+      SCIP_CALL( SCIPcreateVarBasic(subproblems[epoch], &decision_var, name, 0, 1, 0, SCIP_VARTYPE_CONTINUOUS));
+      SCIP_CALL( SCIPaddVar(subproblems[epoch], decision_var));
+
+      // edge capacity constraints
+      SCIPsnprintf( name, SCIP_MAXSTRLEN, "capacity_cons_%d", edge);
+
+      SCIP_CALL( SCIPcreateConsBasicLinear(subproblems[epoch], &capacity_cons, name, 0, NULL, NULL, -SCIPinfinity(subproblems[epoch]), capacities[edge]));
+      SCIP_CALL( SCIPaddCoefLinear(subproblems[epoch], capacity_cons, flow_var, 1));
+      SCIP_CALL( SCIPaddCoefLinear(subproblems[epoch], capacity_cons, decision_var, capacities[edge]));
+      SCIP_CALL( SCIPaddCons(subproblems[epoch], capacity_cons));
+      SCIP_CALL( SCIPreleaseCons(subproblems[epoch], &capacity_cons));
+
+      SCIP_CALL( SCIPreleaseVar(subproblems[epoch], &flow_var));
+      SCIP_CALL( SCIPreleaseVar(subproblems[epoch], &decision_var));
+    }
+
+    // abusing flow var to store target var pointer
+    SCIP_CALL( SCIPcreateVarBasic(subproblems[epoch], &flow_var, "target_variable", 0, target_upper_bound, 0, SCIP_VARTYPE_CONTINUOUS));
+    SCIP_CALL( SCIPaddVar(subproblems[epoch], flow_var));
+    SCIP_CALL( SCIPaddCoefLinear(subproblems[epoch], flow_cons, flow_var, 1));
+    SCIP_CALL( SCIPaddCons(subproblems[epoch], flow_cons));
+    SCIP_CALL( SCIPreleaseVar(subproblems[epoch], &flow_var));
+    SCIP_CALL( SCIPreleaseCons(subproblems[epoch], &flow_cons));
+  }
+
+  SCIP_CALL( SCIPcreateBendersDefault(master, subproblems, nepochs));
+  SCIP_CALL( SCIPsetBoolParam(master, "constraints/benders/active", TRUE) );
+  SCIP_CALL( SCIPsetBoolParam(master, "constraints/benderslp/active", TRUE) );
+  SCIP_CALL( SCIPsetIntParam(master, "constraints/benders/maxprerounds", 1) );
+  SCIP_CALL( SCIPsetIntParam(master, "presolving/maxrounds", 1) );
+
+  container->cons_separator = master;
+  container->vars = mp_vars;
+  container->maintenance_select_conss = maintenance_select_conss;
+
+  return SCIP_OKAY;
+}
+
+SCIP_RETCODE create_cons_separator_scip(SCIP_ConshdlrData* data, SCIP_SOL* solution, cons_separator_container* container){
+  SCIP* separator = data->separator;
+  SCIP_VAR** separator_vars = data->separator_vars;
+  SCIP_Real* capacities = data->capacities;
+  int nedges = data->nedges;
+  int nepochs = data->nepochs;
+  SCIP_Real target_upper_bound = data->target_upper_bound;
+  SCIP_Bool use_reopt = data->use_reopt;