diff --git a/CMakeLists.txt b/CMakeLists.txt
index e0c311c52..3265ecb4f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.22)
-project(idol VERSION 0.10.5)
+project(idol VERSION 0.10.6)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${idol_SOURCE_DIR}/cmake")
diff --git a/bin/method-managers/robust/BBBB.cpp b/bin/method-managers/robust/BBBB.cpp
index 59ad5397b..7caf0d917 100644
--- a/bin/method-managers/robust/BBBB.cpp
+++ b/bin/method-managers/robust/BBBB.cpp
@@ -3,6 +3,7 @@
 //
 #include "BBBB.h"
 #include "../milp/MILPMethodManager.h"
+#include "idol/bilevel/optimizers/wrappers/MibS/MibS.h"
 #include "idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h"
 #include "idol/mixed-integer/modeling/expressions/operations/operators.h"
 #include "idol/mixed-integer/optimizers/branch-and-bound/BranchAndBound.h"
@@ -11,7 +12,82 @@
 #include "idol/mixed-integer/optimizers/callbacks/heuristics/IntegerMaster.h"
 #include "idol/mixed-integer/optimizers/callbacks/heuristics/RENS.h"
 #include "idol/mixed-integer/optimizers/callbacks/heuristics/SimpleRounding.h"
+#include "idol/mixed-integer/optimizers/wrappers/Gurobi/Gurobi.h"
+#include "idol/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/ColumnAndConstraintGeneration.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/separation/OptimalitySeparation.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.h"
 #include "idol/robust/optimizers/critical-value-column-and-constraint-generation/CriticalValueColumnAndConstraintGeneration.h"
+#include "idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h"
+
+template<class NodeInfoT = idol::DefaultNodeInfo>
+class EvaluateCallback : public idol::BranchAndBoundCallbackFactory<NodeInfoT> {
+public:
+    class Strategy : public idol::BranchAndBoundCallback<NodeInfoT> {
+        const unsigned int m_max_evaluations = 10;
+        const unsigned int m_node_frequency = 20;
+    protected:
+        void operator()(idol::CallbackEvent t_event) override;
+    };
+
+    idol::BranchAndBoundCallback<NodeInfoT>* operator()() override { return new Strategy(); }
+    [[nodiscard]] idol::BranchAndBoundCallbackFactory<NodeInfoT>* clone() const override { return new EvaluateCallback<NodeInfoT>(*this); }
+};
+
+template <class NodeInfoT>
+void EvaluateCallback<NodeInfoT>::Strategy::operator()(idol::CallbackEvent t_event) {
+
+    if (t_event != idol::InvalidSolution) {
+        return;
+    }
+
+    if (this->node_count() % m_node_frequency != 0) {
+        return;
+    }
+
+    const auto& max_min_relaxation = this->relaxation().optimizer().template as<idol::Optimizers::Robust::MaxMinRelaxation>();
+    const auto& cvccg = max_min_relaxation.get_formulation().model.optimizer().template as<idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration>();
+    const auto& branching_candidates = cvccg.branching_candidates();
+
+    auto relaxation_fixed = this->relaxation().copy();
+    relaxation_fixed.optimizer().set_param_time_limit(1e-2);
+    relaxation_fixed.optimize();
+
+    unsigned int n_evaluations = 0;
+    for (const auto& uncertainty : cvccg.get_formulation().uncertainties()) {
+
+        std::vector<idol::CVCCG::Formulation::CurrentlyPresentCut> cuts;
+        std::copy(uncertainty.currently_present_cuts().begin(), uncertainty.currently_present_cuts().end(), std::back_inserter(cuts));
+        std::sort(cuts.begin(), cuts.end(), [](const auto& t_a, const auto& t_b) {
+            return t_a.scenario->scenario.objective_value() < t_b.scenario->scenario.objective_value();
+        });
+
+        for (const auto& cut : cuts) {
+
+            for (const auto& var : branching_candidates) {
+                const double val = cut.scenario->scenario.get(var);
+                relaxation_fixed.set_var_lb(var, val);
+                relaxation_fixed.set_var_ub(var, val);
+            }
+            const double remaining_time = this->original_model().optimizer().get_remaining_time();
+            relaxation_fixed.optimizer().set_param_time_limit(remaining_time);
+            relaxation_fixed.optimize();
+
+            const auto status = relaxation_fixed.get_status();
+            if (status == idol::Optimal || status == idol::Feasible) {
+                auto* info = new idol::DefaultNodeInfo();
+                info->set_primal_solution(idol::save_primal(relaxation_fixed));
+                this->submit_heuristic_solution(info);
+            }
+
+            n_evaluations++;
+            if (n_evaluations >= m_max_evaluations) {
+                return;
+            }
+        }
+    }
+
+}
 
 std::string RobustMethods::BBBB::description() const {
     return "Bilevel-based branch-and-bound.";
@@ -51,6 +127,34 @@ void RobustMethods::BBBB::set_optimizer(idol::Model& t_model, const RobustMethod
 
     const auto sub_milp_optimizer = MILPMethodManager::get_sub_milp_optimizer(args);
 
+    for (const auto& var : t_model.vars()) {
+        if (t_model.get_var_type(var) == idol::Continuous) {
+            t_model.set_var_type(var, idol::Integer);
+        }
+    }
+
+    std::list<idol::PrimalPoint> initial_scenarios;
+
+    if (false) {
+
+        auto ccg = idol::Robust::ColumnAndConstraintGeneration(robust_description, bilevel_description);
+        ccg.with_master_optimizer(idol::Gurobi());
+        ccg.add_separation(idol::Robust::CCG::OptimalitySeparation().with_bilevel_optimizer(idol::Bilevel::MibS()));
+        ccg.with_logs(true);
+        ccg.with_iteration_limit(2);
+
+        t_model.use(ccg);
+        t_model.optimize();
+
+        const auto& optimizer = t_model.optimizer().as<idol::Optimizers::Robust::ColumnAndConstraintGeneration>();
+        for (const auto& scenario : optimizer.get_formulation().generated_scenarios()) {
+            initial_scenarios.emplace_back(scenario);
+        }
+
+        std::cout << "Initializing with " << initial_scenarios.size() << " scenarios." << std::endl;
+
+    }
+
     // Branching Candidates
     std::list<idol::Var> branching_candidates;
     for (const auto& var : t_model.vars()) {
@@ -65,10 +169,16 @@ void RobustMethods::BBBB::set_optimizer(idol::Model& t_model, const RobustMethod
     branch_and_bound.with_logs(true);
     branch_and_bound.with_logger(idol::Logs::BranchAndBound::Info().with_frequency_in_seconds(0).with_node_logs(false));
 
+    branch_and_bound.add_callback(EvaluateCallback());
+
     auto max_min_relaxation = idol::Robust::MaxMinRelaxation(robust_description, bilevel_description);
     max_min_relaxation.with_master_optimizer(*sub_milp_optimizer);
     max_min_relaxation.with_deterministic_optimizer(*sub_milp_optimizer);
     max_min_relaxation.with_indicator(false);
 
+    for (const auto& scenario : initial_scenarios) {
+        max_min_relaxation.add_initial_scenario(scenario);
+    }
+
     t_model.use(branch_and_bound + max_min_relaxation);
 }
diff --git a/lib/include/idol/general/utils/SparseVector.h b/lib/include/idol/general/utils/SparseVector.h
index f7c6ea88b..aba91b76c 100644
--- a/lib/include/idol/general/utils/SparseVector.h
+++ b/lib/include/idol/general/utils/SparseVector.h
@@ -120,6 +120,8 @@ class idol::SparseVector {
 
     void clear() { m_map.clear(); }
 
+    std::pair<double, double> range() const;
+
     void reserve(unsigned int t_capacity) {
 #ifdef IDOL_USE_TSL
         m_map.reserve(t_capacity);
@@ -234,6 +236,17 @@ bool idol::SparseVector<IndexT, ValueT>::is_zero(double t_tolerance) const {
     return true;
 }
 
+template <class IndexT, class ValueT>
+std::pair<double, double> idol::SparseVector<IndexT, ValueT>::range() const {
+    double min = Inf, max = -Inf;
+    for (const auto& [var, val] : m_map) {
+        const double abs_val = std::abs(val);
+        min = std::min(min, abs_val);
+        max = std::max(max, abs_val);
+    }
+    return std::make_pair(min, max);
+}
+
 template<class IndexT, class ValueT>
 idol::SparseVector<IndexT, ValueT>&
 idol::SparseVector<IndexT, ValueT>::operator+=(const SparseVector &t_vector) {
diff --git a/lib/include/idol/mixed-integer/modeling/models/Model.h b/lib/include/idol/mixed-integer/modeling/models/Model.h
index 9ff822fe3..74b83520f 100644
--- a/lib/include/idol/mixed-integer/modeling/models/Model.h
+++ b/lib/include/idol/mixed-integer/modeling/models/Model.h
@@ -303,6 +303,8 @@ class idol::Model {
     void reset_minor_representation();
 
     static Model read_from_file(Env& t_env, const std::string& t_filename);
+
+    void print_statistics(std::ostream& t_os = std::cout) const;
 };
 
 template<class T, unsigned int N>
diff --git a/lib/include/idol/mixed-integer/optimizers/branch-and-bound/Optimizers_BranchAndBound.h b/lib/include/idol/mixed-integer/optimizers/branch-and-bound/Optimizers_BranchAndBound.h
index ac07f5034..3b537f6f0 100644
--- a/lib/include/idol/mixed-integer/optimizers/branch-and-bound/Optimizers_BranchAndBound.h
+++ b/lib/include/idol/mixed-integer/optimizers/branch-and-bound/Optimizers_BranchAndBound.h
@@ -466,7 +466,7 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::submit_heuristic_solution(Node
     }
 
     set_as_incumbent(t_node);
-    log_node_after_solve(t_node);
+    //log_node_after_solve(t_node);
 
     //if (m_branching_rule->is_valid(t_node)) {
         // New incumbent by submission
@@ -585,6 +585,10 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_before_optimize() {
 template<class NodeInfoT>
 void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_optimize() {
 
+    if (get_param_logs()) {
+        parent().print_statistics(std::cout);
+    }
+
     if (!m_presolve.empty()) {
         m_presolved_model.reset(working_model().clone());
         m_presolve.execute(*m_presolved_model);
diff --git a/lib/include/idol/mixed-integer/optimizers/branch-and-bound/callbacks/BranchAndBoundCallback.h b/lib/include/idol/mixed-integer/optimizers/branch-and-bound/callbacks/BranchAndBoundCallback.h
index 6845d80a6..da5be03cb 100644
--- a/lib/include/idol/mixed-integer/optimizers/branch-and-bound/callbacks/BranchAndBoundCallback.h
+++ b/lib/include/idol/mixed-integer/optimizers/branch-and-bound/callbacks/BranchAndBoundCallback.h
@@ -324,21 +324,27 @@ idol::BranchAndBoundCallbackI<NodeInfoT>::operator()(Optimizers::BranchAndBound<
                                                      Model *t_relaxation) {
     SideEffectRegistry result;
 
-    m_parent = t_parent;
-    m_node = t_current_node;
-    m_relaxation = t_relaxation;
-    m_registry = &result;
+    const bool is_nested = m_parent != nullptr;
+
+    if (!is_nested) {
+        m_parent = t_parent;
+        m_node = t_current_node;
+        m_relaxation = t_relaxation;
+        m_registry = &result;
+    }
 
     for (auto &cb: m_callbacks) {
-        cb->m_interface = this;
+        if (!is_nested) { cb->m_interface = this; }
         cb->operator()(t_event);
-        cb->m_interface = nullptr;
+        if (!is_nested) { cb->m_interface = nullptr; }
     }
 
-    m_parent = nullptr;
-    m_node.reset();
-    m_relaxation = nullptr;
-    m_registry = nullptr;
+    if (!is_nested) {
+        m_parent = nullptr;
+        m_node.reset();
+        m_relaxation = nullptr;
+        m_registry = nullptr;
+    }
 
     return result;
 }
diff --git a/lib/include/idol/mixed-integer/optimizers/branch-and-bound/logs/Info.h b/lib/include/idol/mixed-integer/optimizers/branch-and-bound/logs/Info.h
index 1b2f1848f..f6527021f 100644
--- a/lib/include/idol/mixed-integer/optimizers/branch-and-bound/logs/Info.h
+++ b/lib/include/idol/mixed-integer/optimizers/branch-and-bound/logs/Info.h
@@ -147,8 +147,10 @@ void idol::Logs::BranchAndBound::Info<NodeInfoT>::Strategy::log_node_after_solve
 
     const double total_time = parent.time().count();
 
-    if (!m_root_node_has_been_printed) {
-        log_root_node(t_node);
+    if (!m_root_node_has_been_printed && t_node.id() == 0) {
+        if (t_node.id() == 0) {
+            log_root_node(t_node);
+        }
         return;
     }
 
diff --git a/lib/include/idol/mixed-integer/optimizers/callbacks/cutting-planes/CglCutCallback.h b/lib/include/idol/mixed-integer/optimizers/callbacks/cutting-planes/CglCutCallback.h
index 358c59d2b..760cbde73 100644
--- a/lib/include/idol/mixed-integer/optimizers/callbacks/cutting-planes/CglCutCallback.h
+++ b/lib/include/idol/mixed-integer/optimizers/callbacks/cutting-planes/CglCutCallback.h
@@ -65,7 +65,8 @@ class idol::CglCutCallback : public BranchAndBoundCallbackFactory<NodeInfoT> {
 
         std::list<TempCtr> to_idol_cuts(OsiCuts& t_cuts);
         TempCtr to_idol_cut(OsiRowCut& t_cut);
-        std::vector<std::pair<TempCtr, double>> sort_cuts_by_effectiveness(const std::list<TempCtr>& t_cuts);
+        std::vector<std::pair<TempCtr, double>> sort_and_filter_cuts_by_effectiveness(const std::list<TempCtr>& t_cuts);
+        void standard_scaling(std::list<TempCtr>& t_cuts);
     protected:
         NodeCutContext& get_cut_context();
         const NodeCutContext& get_cut_context() const { return const_cast<Strategy*>(this)->get_cut_context(); }
@@ -184,7 +185,8 @@ void idol::CglCutCallback<NodeInfoT>::Strategy::operator()(CallbackEvent t_event
 
         auto osi_cuts = cut_family->generate(*m_osi_solver, be_aggressive ? 100 : 0);
         auto idol_cuts = to_idol_cuts(osi_cuts);
-        auto sorted_cuts = sort_cuts_by_effectiveness(idol_cuts);
+        standard_scaling(idol_cuts);
+        auto sorted_cuts = sort_and_filter_cuts_by_effectiveness(idol_cuts);
 
         for (auto& [cut, effectiveness] : sorted_cuts) {
 
@@ -290,7 +292,7 @@ idol::TempCtr idol::CglCutCallback<NodeInfoT>::Strategy::to_idol_cut(OsiRowCut&
 }
 
 template <class NodeInfoT>
-std::vector<std::pair<idol::TempCtr, double>> idol::CglCutCallback<NodeInfoT>::Strategy::sort_cuts_by_effectiveness(
+std::vector<std::pair<idol::TempCtr, double>> idol::CglCutCallback<NodeInfoT>::Strategy::sort_and_filter_cuts_by_effectiveness(
     const std::list<TempCtr>& t_cuts) {
     std::vector<std::pair<TempCtr, double>> result;
     const auto primal_solution = this->node().info().primal_solution();
@@ -320,6 +322,34 @@ std::vector<std::pair<idol::TempCtr, double>> idol::CglCutCallback<NodeInfoT>::S
     return result;
 }
 
+template <class NodeInfoT>
+void idol::CglCutCallback<NodeInfoT>::Strategy::standard_scaling(std::list<TempCtr>& t_cuts) {
+
+    for (auto& cut : t_cuts) {
+
+        auto& row = cut.lhs();
+        double& rhs = cut.rhs();
+        double infinity_norm = 0;
+        for (const auto& [var, coeff] : row) {
+            infinity_norm = std::max(std::abs(coeff), infinity_norm);
+        }
+        infinity_norm = std::max(infinity_norm, std::abs(rhs));
+
+        if (is_zero(infinity_norm, Tolerance::Sparsity)) {
+            continue;
+        }
+
+        int e = 0;
+        std::frexp(infinity_norm, &e);
+        double closest_power_of_2 = std::ldexp(1.0, e - 1); // scale = 2^(e-1) or 2^e depending on your normalization choice
+        if (closest_power_of_2 != 1.) {
+            row /= closest_power_of_2;
+            rhs /= closest_power_of_2;
+        }
+
+    }
+}
+
 template <class NodeInfoT>
 idol::CglCutCallback<NodeInfoT>::Strategy::NodeCutContext& idol::CglCutCallback<
     NodeInfoT>::Strategy::get_cut_context() {
diff --git a/lib/include/idol/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.h b/lib/include/idol/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.h
index e5dc5a618..27bbb1993 100644
--- a/lib/include/idol/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.h
+++ b/lib/include/idol/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.h
@@ -8,18 +8,20 @@
 #include "Optimizers_DantzigWolfeDecomposition.h"
 
 class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
+
+    enum NumericalPolicy { Default, ColumnPoolCleanUp, NoDualSmoothing, Failure };
+
     DantzigWolfeDecomposition& m_parent;
     double m_best_bound_stop;
 
-    enum NumericalPolicy { Default, ColumnPoolCleanUp, NoDualSmoothing, Failure };
-    const unsigned int m_max_n_iterations_without_generating_column = 1000;
+    const unsigned int m_max_n_iterations_without_generating_column = 200;
     unsigned int m_n_iterations_without_generating_column = 0;
     NumericalPolicy m_numerical_policy = Default;
 
     SolutionStatus m_status = Loaded;
     SolutionReason m_reason = NotSpecified;
     std::optional<PrimalPoint> m_master_primal_solution;
-    std::optional<DualPoint> m_last_master_solution;
+    std::optional<DualPoint> m_master_dual_solution;
     std::vector<DantzigWolfe::SubProblem::PhaseId> m_sub_problems_phases;
     double m_best_obj = -Inf;
     double m_best_bound = +Inf;
@@ -39,6 +41,7 @@ class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
     void solve_sub_problems_in_parallel();
     void analyze_sub_problems();
     void enrich_master();
+    bool check_numerical_stability();
     void pool_clean_up();
 
     void next_numerical_policy();
@@ -47,6 +50,8 @@ class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
     void log_master();
     void log_sub_problems();
     void log_end();
+
+    friend std::ostream& operator<<(std::ostream& t_os, idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::NumericalPolicy t_numerical_policy);
 public:
     ColumnGeneration(DantzigWolfeDecomposition& t_parent, bool t_use_farkas_for_infeasibility, double t_best_bound_stop);
 
@@ -69,4 +74,5 @@ class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
     void execute();
 };
 
+
 #endif //IDOL_COLUMNGENERATION_H
diff --git a/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h b/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h
index 3f1089d68..bab63a2be 100644
--- a/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h
+++ b/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h
@@ -19,6 +19,7 @@ class idol::Robust::MaxMinRelaxation : public OptimizerFactoryWithDefaultParamet
     std::unique_ptr<OptimizerFactory> m_master_optimizer_factory;
     std::unique_ptr<OptimizerFactory> m_deterministic_optimizer_factory;
     std::optional<bool> m_use_indicator;
+    std::list<PrimalPoint> m_initial_scenarios;
 protected:
     [[nodiscard]] Optimizer* create(const Model& t_model) const override;
     MaxMinRelaxation(const MaxMinRelaxation& t_src);
@@ -30,6 +31,7 @@ class idol::Robust::MaxMinRelaxation : public OptimizerFactoryWithDefaultParamet
     MaxMinRelaxation& with_master_optimizer(const OptimizerFactory& t_optimizer);
     MaxMinRelaxation& with_deterministic_optimizer(const OptimizerFactory& t_optimizer);
     MaxMinRelaxation& with_indicator(bool t_value);
+    MaxMinRelaxation& add_initial_scenario(PrimalPoint t_scenario);
 };
 
 #endif //IDOL_ROBUST_BBBB_H
\ No newline at end of file
diff --git a/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.h b/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.h
index 9349a6407..39e252ee8 100644
--- a/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.h
+++ b/lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.h
@@ -23,6 +23,7 @@ class idol::Optimizers::Robust::MaxMinRelaxation : public Algorithm {
     std::unique_ptr<OptimizerFactory> m_master_optimizer_factory;
     std::unique_ptr<OptimizerFactory> m_deterministic_optimizer_factory;
     const bool m_use_indicator;
+    const std::list<PrimalPoint> m_initial_scenarios;
 
     struct Formulation {
         Model model;
@@ -38,7 +39,8 @@ class idol::Optimizers::Robust::MaxMinRelaxation : public Algorithm {
                      const idol::Bilevel::Description& t_bilevel_description,
                      const OptimizerFactory& t_master_optimizer,
                      const OptimizerFactory& t_deterministic_optimizer,
-                     bool t_use_indicator
+                     bool t_use_indicator,
+                     const std::list<PrimalPoint>& t_initial_scenarios
     );
 
     [[nodiscard]] std::string name() const override { return "max-min relaxation"; }
@@ -54,6 +56,9 @@ class idol::Optimizers::Robust::MaxMinRelaxation : public Algorithm {
     [[nodiscard]] const OptimizerFactory& get_master_optimizer_factory() const { return *m_master_optimizer_factory; }
     [[nodiscard]] const OptimizerFactory& get_deterministic_optimizer_factory() const { return *m_deterministic_optimizer_factory; }
     [[nodiscard]] bool use_indicator() const { return m_use_indicator; }
+    [[nodiscard]] const Formulation& get_formulation() const { return *m_formulation; }
+
+    void build_model();
 protected:
     void add(const Var& t_var) override THROW_NOT_IMPLEMENTED
     void add(const Ctr& t_ctr) override THROW_NOT_IMPLEMENTED
@@ -80,7 +85,6 @@ class idol::Optimizers::Robust::MaxMinRelaxation : public Algorithm {
     void update_var_ub(const Var& t_var) override;
     void update_var_obj(const Var& t_var) override THROW_NOT_IMPLEMENTED
 
-    void build_model();
     void throw_if_no_formulation() const;
 };
 
diff --git a/lib/include/idol/robust/optimizers/column-and-constraint-generation/Formulation.h b/lib/include/idol/robust/optimizers/column-and-constraint-generation/Formulation.h
index 5f527c290..39c781adb 100644
--- a/lib/include/idol/robust/optimizers/column-and-constraint-generation/Formulation.h
+++ b/lib/include/idol/robust/optimizers/column-and-constraint-generation/Formulation.h
@@ -57,6 +57,7 @@ class idol::CCG::Formulation {
     const Model& original_model() const { return m_parent; }
     bool has_second_stage_epigraph() const { return m_second_stage_epigraph.has_value(); }
     bool has_second_stage_objective() const { return m_has_second_stage_objective; }
+    auto generated_scenarios() const { return ConstIteratorForward(m_generated_scenarios); }
 
 };
 
diff --git a/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.h b/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.h
index 8771b680d..222d4e99b 100644
--- a/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.h
+++ b/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.h
@@ -19,7 +19,7 @@ namespace idol::CVCCG {
 }
 
 class idol::CVCCG::Formulation {
-
+public:
     struct LinkingConstraint {
         Ctr ctr_in_uncertainty_set;
         Map<long int, std::pair<Var, Ctr>> critical_values;
@@ -54,7 +54,7 @@ class idol::CVCCG::Formulation {
         Uncertainty() = default;
         Uncertainty(const Ctr& t_ctr) : m_constraint_in_original_model(t_ctr) {}
     };
-
+private:
     const Optimizers::Robust::CriticalValueColumnAndConstraintGeneration& m_parent;
 
     // Analysis of the model
@@ -63,6 +63,7 @@ class idol::CVCCG::Formulation {
     bool m_all_linking_variables_are_binary = true;
     bool m_all_data_in_linking_constraints_is_integer = true;
     const bool m_use_cover_constraints = true;
+    double m_global_lower_bound = -Inf;
 
     Model m_master;
     Model m_sub_problem;
@@ -82,6 +83,8 @@ class idol::CVCCG::Formulation {
     void create_critical_value_variable_if_needed(const PrimalPoint& t_scenario);
     void create_critical_value_variable(const PrimalPoint& t_scenario, LinkingConstraint& t_linking);
     double compute_critical_value(const Ctr& t_ctr, const PrimalPoint& t_scenario) const;
+    double compute_trivial_lower_bound() const;
+    double compute_penalty(const LinExpr<Var, double>& t_row, CtrType t_type, double t_rhs) const;
 public:
     Formulation(const idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration& t_parent);
 
@@ -110,6 +113,7 @@ class idol::CVCCG::Formulation {
     void load_cut_from_pool();
 
     bool uses_indicator() const;
+
 };
 
 #endif //IDOL_CVCCG_FORMULATION_H
\ No newline at end of file
diff --git a/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h b/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h
index 6a58d553b..ff195ab5d 100644
--- a/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h
+++ b/lib/include/idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h
@@ -23,6 +23,7 @@ class idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration : pub
     std::unique_ptr<CVCCG::Formulation> m_formulation;
     std::list<Var> m_branching_candidates;
     unsigned int m_n_iterations = 0;
+    bool m_is_diving = false;
 public:
     CriticalValueColumnAndConstraintGeneration(const Model& t_model,
                                                const idol::Robust::Description& t_description,
@@ -45,6 +46,8 @@ class idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration : pub
     [[nodiscard]] const OptimizerFactory& get_deterministic_optimizer_factory() const { return *m_deterministic_optimizer_factory; }
     [[nodiscard]] bool use_indicator() const { return m_use_indicator; }
 
+    [[nodiscard]] const CVCCG::Formulation& get_formulation() const { return *m_formulation; }
+
     void declare_branching_on_unc_vars(std::list<Var> t_branching_candidates) { m_branching_candidates = std::move(t_branching_candidates); }
     [[nodiscard]] auto branching_candidates() const { return ConstIteratorForward(m_branching_candidates); }
     void set_unc_var_lb(const Var& t_var, double t_lb);
diff --git a/lib/src/mixed-integer/modeling/models/Model.cpp b/lib/src/mixed-integer/modeling/models/Model.cpp
index 8caf2748d..85ca2173d 100644
--- a/lib/src/mixed-integer/modeling/models/Model.cpp
+++ b/lib/src/mixed-integer/modeling/models/Model.cpp
@@ -1047,6 +1047,61 @@ idol::Model idol::Model::read_from_file(Env& t_env, const std::string& t_filenam
     throw Exception("Reading .mps or .lp files is only possible if Gurobi, GLPK or HiGHS is available.\nPlease install one of these solvers.");
 }
 
+void idol::Model::print_statistics(std::ostream& t_os) const {
+
+    unsigned int n_integer_vars = 0;
+    unsigned int n_binary_vars = 0;
+    unsigned int n_non_zeroes = 0;
+    double min_mat_val = Inf, min_bound_val = Inf;
+    double max_mat_val = -Inf, max_bound_val = -Inf;
+    const auto [min_obj_val, max_obj_val] = get_obj_expr().affine().linear().range();
+    const auto [min_rhs_val, max_rhs_val] = get_rhs_expr().range();
+
+    for (const auto& var : vars()) {
+        const auto type = get_var_type(var);
+        const auto& column = get_var_column(var);
+        const double lb = get_var_lb(var);
+        const double ub = get_var_ub(var);
+
+        if (type != Continuous) {
+            n_integer_vars++;
+            if (type == Binary) {
+                n_binary_vars++;
+            }
+        }
+
+        if (!is_neg_inf(lb) && !is_zero(lb, Tolerance::Sparsity)) {
+            const double abs_lb = std::abs(lb);
+            max_bound_val = std::max(max_bound_val, abs_lb);
+            min_bound_val = std::min(min_bound_val, abs_lb);
+        }
+
+        if (!is_pos_inf(ub) && !is_zero(ub, Tolerance::Sparsity)) {
+            const double abs_ub = std::abs(ub);
+            max_bound_val = std::max(max_bound_val, abs_ub);
+            min_bound_val = std::min(min_bound_val, abs_ub);
+        }
+
+        n_non_zeroes += column.size();
+
+        for (const auto& [ctr, coeff] : column) {
+            const double abs_coeff = std::abs(coeff);
+            min_mat_val = std::min(min_mat_val, abs_coeff);
+            max_mat_val = std::max(max_mat_val, abs_coeff);
+        }
+
+    }
+
+    t_os << "Model with " << m_constraints.size() << " row, " << m_variables.size() << " columns and " << n_non_zeroes << " nonzeroes" << std::endl;
+    t_os << "Variable types: " << m_variables.size() - n_integer_vars << " continuous, " << n_integer_vars << " integer (" << n_binary_vars << " binary)" << std::endl;
+    t_os << "Coefficient statistcs:\n"
+            "\tMatrix range\t[" << min_mat_val << ", " << max_mat_val << "]\n"
+            "\tObjective range\t[" << min_obj_val << ", " << max_obj_val << "]\n"
+            "\tBounds range\t[" << min_bound_val << ", " << max_bound_val << "]\n"
+            "\tRHS range\t\t[" << min_rhs_val << ", " << max_rhs_val << "]\n";
+
+}
+
 void idol::Model::build_rows() {
     for (const auto& ctr : m_constraints) {
         if (!m_env.version(*this, ctr).has_row()) {
diff --git a/lib/src/mixed-integer/optimizers/branch-and-bound/CutPool.cpp b/lib/src/mixed-integer/optimizers/branch-and-bound/CutPool.cpp
index 66a39b16d..89bec7aac 100644
--- a/lib/src/mixed-integer/optimizers/branch-and-bound/CutPool.cpp
+++ b/lib/src/mixed-integer/optimizers/branch-and-bound/CutPool.cpp
@@ -33,6 +33,7 @@ bool idol::CutPool::add_existing_cut_to_relaxation(const Ctr& t_cut, Model& t_re
 }
 
 bool idol::CutPool::add_cut(const TempCtr& t_cut, Model& t_relaxation) {
+
     Ctr cut(t_relaxation.env(), t_cut);
     m_all_cuts.emplace_back(cut);
 
diff --git a/lib/src/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.cpp b/lib/src/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.cpp
index 9c03fbecf..a5da0ce62 100644
--- a/lib/src/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.cpp
+++ b/lib/src/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.cpp
@@ -17,7 +17,7 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::execute() {
 
     m_status = Loaded;
     m_reason = NotSpecified;
-    m_last_master_solution.reset();
+    m_master_dual_solution.reset();
     m_iteration_count = 0;
     m_n_generated_columns = 0;
     m_solve_dual_master = true;
@@ -54,6 +54,8 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::execute() {
 
         log_sub_problems();
 
+        if (check_numerical_stability()) { continue; }
+
         if (check_stopping_criterion()) { break; }
 
         pool_clean_up();
@@ -93,7 +95,7 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::solve_dual_m
         m_best_obj = std::min(m_best_obj, iter_upper_bound);
 
         if (save_dual_solution) {
-            m_last_master_solution = save_dual(master);
+            m_master_dual_solution = save_dual(master);
         }
 
         return;
@@ -107,13 +109,13 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::solve_dual_m
         m_current_iteration_is_using_farkas = true;
 
         if (save_dual_solution) {
-            m_last_master_solution = save_farkas(master);
+            m_master_dual_solution = save_farkas(master);
         }
 
         return;
     }
 
-    m_last_master_solution.reset();
+    m_master_dual_solution.reset();
 
     m_is_terminated = true;
 
@@ -122,9 +124,10 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::solve_dual_m
 void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::update_sub_problems() {
 
     auto& formulation = m_parent.m_formulation;
+    assert(m_master_dual_solution.has_value());
     auto dual_values = m_current_iteration_is_using_farkas || m_numerical_policy >= NoDualSmoothing ?
-        m_last_master_solution.value() :
-        m_parent.m_stabilization->compute_smoothed_dual_solution(m_last_master_solution.value());
+        m_master_dual_solution.value() :
+        m_parent.m_stabilization->compute_smoothed_dual_solution(m_master_dual_solution.value());
 
     for (unsigned int i = 0, n = formulation.n_sub_problems() ; i < n ; ++i) {
         formulation.update_sub_problem_objective(i, dual_values, m_current_iteration_is_using_farkas);
@@ -185,11 +188,12 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::analyze_sub_
         return;
     }
 
-    const double iter_lower_bound = m_last_master_solution->objective_value() + sum_reduced_costs;
+    assert(m_master_dual_solution.has_value());
+    const double iter_lower_bound = m_master_dual_solution->objective_value() + sum_reduced_costs;
 
     if (m_best_bound <= iter_lower_bound) {
 
-        m_parent.m_stabilization->update_stability_center(m_last_master_solution.value());
+        m_parent.m_stabilization->update_stability_center(m_master_dual_solution.value());
         m_best_bound = std::min(m_best_obj, iter_lower_bound);
 
     }
@@ -282,7 +286,8 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::enrich_maste
 
             auto generator =  save_primal(model);
 
-            const double reduced_cost = formulation.compute_reduced_cost(i, m_last_master_solution.value(), generator);
+            assert(m_master_dual_solution.has_value());
+            const double reduced_cost = formulation.compute_reduced_cost(i, m_master_dual_solution.value(), generator);
 
             if (reduced_cost < -tol_red_cost) {
                 formulation.generate_column(i, std::move(generator));
@@ -310,46 +315,22 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::enrich_maste
     }
 }
 
-void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::initialize_sub_problem_phases() {
-
-    // TODO handle phase restart here
-
-    auto& formulation = m_parent.m_formulation;
-    unsigned int n_sub_problems = formulation.n_sub_problems();
-
-    m_sub_problems_phases.clear();
-    m_sub_problems_phases.reserve(n_sub_problems);
-    for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
-        const auto phase = m_parent.m_sub_problem_specifications[i].phases().begin();
-        if (formulation.sub_problem(i).has_optimizer()) { continue; }
-        formulation.sub_problem(i).use(**phase);
-        m_sub_problems_phases.emplace_back(phase);
-    }
-
-}
-
-void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::pool_clean_up() {
+bool idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::check_numerical_stability() {
 
     auto& formulation = m_parent.m_formulation;
     const auto n_sub_problems = formulation.n_sub_problems();
 
-    std::optional<PrimalPoint> primal_solution;
-    const auto& get_master_primal = [&]() {
-        if (primal_solution) {
-            return *primal_solution;
-        }
-        const auto& master = formulation.master();
-        if (const auto status = master.get_status() ; status == Optimal || status == Feasible) {
-            primal_solution = save_primal(master);
-        } else {
-            primal_solution = PrimalPoint();
-        }
-        return *primal_solution;
-    };
+    if (m_best_bound > m_best_obj + 1e-4) {
+        m_status = Fail;
+        m_reason = Numerical;
+        m_is_terminated = true;
+        return false;
+    }
 
     // Check max iteration without generating column
     if (m_n_iterations_without_generating_column >= m_max_n_iterations_without_generating_column) {
         next_numerical_policy();
+        return true;
     }
 
     // Check for cg being stuck
@@ -377,8 +358,49 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::pool_clean_u
 
     if (is_stuck) {
         next_numerical_policy();
+        return true;
     }
 
+    return false;
+}
+
+void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::initialize_sub_problem_phases() {
+
+    // TODO handle phase restart here
+
+    auto& formulation = m_parent.m_formulation;
+    unsigned int n_sub_problems = formulation.n_sub_problems();
+
+    m_sub_problems_phases.clear();
+    m_sub_problems_phases.reserve(n_sub_problems);
+    for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
+        const auto phase = m_parent.m_sub_problem_specifications[i].phases().begin();
+        if (formulation.sub_problem(i).has_optimizer()) { continue; }
+        formulation.sub_problem(i).use(**phase);
+        m_sub_problems_phases.emplace_back(phase);
+    }
+
+}
+
+void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::pool_clean_up() {
+
+    auto& formulation = m_parent.m_formulation;
+    const auto n_sub_problems = formulation.n_sub_problems();
+
+    std::optional<PrimalPoint> primal_solution;
+    const auto& get_master_primal = [&]() {
+        if (primal_solution) {
+            return *primal_solution;
+        }
+        const auto& master = formulation.master();
+        if (const auto status = master.get_status() ; status == Optimal || status == Feasible) {
+            primal_solution = save_primal(master);
+        } else {
+            primal_solution = PrimalPoint();
+        }
+        return *primal_solution;
+    };
+
     for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
         const auto& sub_problem_specifications = m_parent.m_sub_problem_specifications[i];
         const auto& column_pool = formulation.column_pool(i);
@@ -410,6 +432,11 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::next_numeric
         for (unsigned int i = 0 ; i < n_sub_problems ; ++i) {
             m_parent.m_formulation.clean_up(i, .0, primals, false);
         }
+
+        m_status = Infeasible;
+        m_best_bound = -Inf;
+        m_best_obj = Inf;
+
         return;
     }
 
@@ -435,7 +462,7 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::log_init() {
 
 void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::log_master() {
 
-    if (!m_parent.get_param_logs() || !m_parent.m_logger || !m_last_master_solution.has_value()) {
+    if (!m_parent.get_param_logs() || !m_parent.m_logger || !m_master_dual_solution.has_value()) {
         return;
     }
 
@@ -445,9 +472,9 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::log_master()
                 m_iteration_count,
                 m_parent.time().count(),
                 m_status,
-                m_last_master_solution->status(),
-                m_last_master_solution->reason(),
-                m_last_master_solution->objective_value(),
+                m_master_dual_solution->status(),
+                m_master_dual_solution->reason(),
+                m_master_dual_solution->objective_value(),
                 formulation.master().optimizer().time().count(),
                 m_best_bound,
                 m_best_obj,
@@ -493,3 +520,13 @@ void idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::log_end() {
 
     m_parent.m_logger->log_end();
 }
+
+std::ostream& idol::Optimizers::operator<<(std::ostream& t_os, DantzigWolfeDecomposition::ColumnGeneration::NumericalPolicy t_numerical_policy) {
+    switch (t_numerical_policy) {
+    case idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::Default: return t_os << "default";
+    case idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::ColumnPoolCleanUp: return t_os << "clean up";
+    case idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::NoDualSmoothing: return t_os << "no dual smoothing";
+    case idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::Failure: return t_os << "failure";
+    }
+    throw idol::Exception("Enum out of bounds.");
+}
diff --git a/lib/src/mixed-integer/optimizers/presolve/StandardScaling.cpp b/lib/src/mixed-integer/optimizers/presolve/StandardScaling.cpp
index 5bd26bf10..70476bc40 100644
--- a/lib/src/mixed-integer/optimizers/presolve/StandardScaling.cpp
+++ b/lib/src/mixed-integer/optimizers/presolve/StandardScaling.cpp
@@ -20,7 +20,9 @@ bool idol::Presolvers::StandardScaling::execute(Model& t_model) {
             continue;
         }
 
-        const double closest_power_of_2 = std::exp2(std::round(std::log2(infinity_norm)));
+        int e = 0;
+        std::frexp(infinity_norm, &e);
+        double closest_power_of_2 = std::ldexp(1.0, e - 1); // scale = 2^(e-1) or 2^e depending on your normalization choice
         if (closest_power_of_2 != 1.) {
             t_model.set_ctr_row(ctr, row / closest_power_of_2);
             t_model.set_ctr_rhs(ctr, rhs / closest_power_of_2);
diff --git a/lib/src/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.cpp b/lib/src/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.cpp
index be8040ece..73febd881 100644
--- a/lib/src/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.cpp
+++ b/lib/src/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.cpp
@@ -27,7 +27,8 @@ idol::Optimizer* idol::Robust::MaxMinRelaxation::create(const Model& t_model) co
         m_bilevel_description,
         *m_master_optimizer_factory,
         *m_deterministic_optimizer_factory,
-        m_use_indicator.value_or(false)
+        m_use_indicator.value_or(false),
+        m_initial_scenarios
     );
 
     return result;
@@ -39,7 +40,8 @@ idol::Robust::MaxMinRelaxation::MaxMinRelaxation(const MaxMinRelaxation& t_src)
     m_bilevel_description(t_src.m_bilevel_description),
     m_use_indicator(t_src.m_use_indicator),
     m_master_optimizer_factory(t_src.m_master_optimizer_factory ? t_src.m_master_optimizer_factory->clone() : nullptr),
-    m_deterministic_optimizer_factory(t_src.m_deterministic_optimizer_factory ? t_src.m_deterministic_optimizer_factory->clone() : nullptr) {
+    m_deterministic_optimizer_factory(t_src.m_deterministic_optimizer_factory ? t_src.m_deterministic_optimizer_factory->clone() : nullptr),
+    m_initial_scenarios(t_src.m_initial_scenarios) {
 
 }
 
@@ -79,3 +81,10 @@ idol::Robust::MaxMinRelaxation& idol::Robust::MaxMinRelaxation::with_indicator(b
 
     return *this;
 }
+
+idol::Robust::MaxMinRelaxation& idol::Robust::MaxMinRelaxation::add_initial_scenario(PrimalPoint t_scenario) {
+
+    m_initial_scenarios.push_back(std::move(t_scenario));
+
+    return *this;
+}
diff --git a/lib/src/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.cpp b/lib/src/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.cpp
index bccdfe1a9..11b73e2d3 100644
--- a/lib/src/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.cpp
+++ b/lib/src/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.cpp
@@ -11,12 +11,14 @@ idol::Optimizers::Robust::MaxMinRelaxation::MaxMinRelaxation(const Model& t_mode
                                                              const idol::Bilevel::Description& t_bilevel_description,
                                                              const OptimizerFactory& t_master_optimizer,
                                                              const OptimizerFactory& t_deterministic_optimizer,
-                                                             bool t_use_indicator) : Algorithm(t_model),
+                                                             bool t_use_indicator,
+                                                             const std::list<PrimalPoint>& t_initial_scenarios) : Algorithm(t_model),
                                                                                      m_description(t_description),
                                                                                      m_bilevel_description(t_bilevel_description),
                                                                                      m_master_optimizer_factory(t_master_optimizer.clone()),
                                                                                      m_deterministic_optimizer_factory(t_deterministic_optimizer.clone()),
-                                                                                     m_use_indicator(t_use_indicator)
+                                                                                     m_use_indicator(t_use_indicator),
+                                                                                     m_initial_scenarios(t_initial_scenarios)
 {}
 
 double idol::Optimizers::Robust::MaxMinRelaxation::get_var_primal(const Var& t_var) const {
@@ -33,7 +35,7 @@ void idol::Optimizers::Robust::MaxMinRelaxation::hook_optimize() {
     if (!m_formulation) {
         build_model();
     }
-    
+
     auto& reformulation = m_formulation->model;
     
     // Parameters
@@ -59,15 +61,6 @@ void idol::Optimizers::Robust::MaxMinRelaxation::hook_optimize() {
     set_best_bound(-m_formulation->model.get_best_bound());
     set_best_obj(-m_formulation->model.get_best_obj());
 
-    /*
-    for (const auto& var : parent().vars()) {
-        if (var.name()[0] == 'x') {
-            std::cout << var.name() << " = " << get_var_primal(var) << std::endl;
-        }
-    }
-    std::cout << "****************" << std::endl;
-    */
-
 }
 
 idol::Optimizers::Robust::MaxMinRelaxation::Formulation::Formulation(Model&& t_hpr, idol::Robust::Description&& t_description)
@@ -125,6 +118,55 @@ void idol::Optimizers::Robust::MaxMinRelaxation::build_model() {
         }
     }
 
+    if (!m_initial_scenarios.empty()) {
+
+        if (get_param_logs()) {
+            std::cout << "Adding " << m_initial_scenarios.size() << " initial scenarios..." << std::endl;
+        }
+
+        unsigned int k = 0;
+        for (const auto& scenario : m_initial_scenarios) {
+
+            std::vector<std::optional<Var>> vars;
+            vars.resize(original_model.vars().size());
+
+            for (const auto& var : original_model.vars()) {
+                if (m_bilevel_description.is_upper(var)) {
+                    continue;
+                }
+                const double lb = original_model.get_var_lb(var);
+                const double ub = original_model.get_var_ub(var);
+                const auto type = original_model.get_var_type(var);
+                const auto copy = uncertainty_set.add_var(lb, ub, type, 0, var.name() + "_" + std::to_string(k));
+                vars[original_model.get_var_index(var)] = copy;
+            }
+
+            for (const auto& ctr : original_model.ctrs()) {
+                const auto& row = original_model.get_ctr_row(ctr);
+                const auto rhs = original_model.get_ctr_rhs(ctr);
+                const auto type = original_model.get_ctr_type(ctr);
+
+                LinExpr copied_row;
+                for (const auto& [var, coeff] : row) {
+                    const auto index = original_model.get_var_index(var);
+                    copied_row += coeff * vars.at(index).value_or(var);
+                }
+
+                uncertainty_set.add_ctr(TempCtr(std::move(copied_row), type, rhs - evaluate(m_description.uncertain_rhs(ctr), scenario)), ctr.name() + "_" + std::to_string(k));
+            }
+
+            LinExpr copied_obj;
+            for (const auto& [var, coeff] : original_model.get_obj_expr().affine().linear()) {
+                const auto index = original_model.get_var_index(var);
+                copied_obj += coeff * vars.at(index).value_or(var);
+            }
+            uncertainty_set.add_ctr(original_model.get_obj_expr().affine().linear() >= copied_obj, "obj_" + std::to_string(k));
+
+            ++k;
+        }
+
+    }
+
     idol::Robust::Description description(uncertainty_set);
     const auto one = hpr.add_var(1, 1, Continuous, -original_model.get_obj_expr().affine().constant(), "__one");
     description.set_uncertain_obj(one, -original_model.get_obj_expr().affine().linear());
diff --git a/lib/src/robust/optimizers/column-and-constraint-generation/Formulation.cpp b/lib/src/robust/optimizers/column-and-constraint-generation/Formulation.cpp
index 918cfdd3b..04136f00e 100644
--- a/lib/src/robust/optimizers/column-and-constraint-generation/Formulation.cpp
+++ b/lib/src/robust/optimizers/column-and-constraint-generation/Formulation.cpp
@@ -226,8 +226,9 @@ void idol::CCG::Formulation::add_scenario_to_master(const idol::Point<idol::Var>
             }
         }
 
-        m_generated_scenarios.emplace_back(t_scenario);
     }
+
+    m_generated_scenarios.emplace_back(t_scenario);
 }
 
 void idol::CCG::Formulation::copy_bilevel_description(const ::idol::Bilevel::Description& t_src, const ::idol::Bilevel::Description& t_dest) const {
diff --git a/lib/src/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.cpp b/lib/src/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.cpp
index a0cc958c4..9e7ef213c 100644
--- a/lib/src/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.cpp
+++ b/lib/src/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.cpp
@@ -145,6 +145,8 @@ void idol::Optimizers::Robust::ColumnAndConstraintGeneration::hook_optimize() {
 
         solve_adversarial_problem();
 
+        if (m_n_iterations >= get_param_iteration_limit()) { terminate(); }
+
         if (is_terminated()) { break; }
 
         ++m_n_iterations;
@@ -295,7 +297,9 @@ void idol::Optimizers::Robust::ColumnAndConstraintGeneration::check_termination_
     }
 
     if (get_best_bound() > get_best_obj() + get_tol_mip_absolute_gap()) {
-        std::cerr << "The current best bound is larger than current best obj. This should should not happen. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cerr << "The current best bound is larger than current best obj. This should should not happen. Terminating..." << std::endl;
+        }
         set_status(Fail);
         set_reason(Numerical);
         terminate();
@@ -305,21 +309,27 @@ void idol::Optimizers::Robust::ColumnAndConstraintGeneration::check_termination_
     set_status(is_pos_inf(get_best_obj()) ? Infeasible : Feasible);
 
     if (get_remaining_time() == 0) {
-        std::cout << "The time limit has been reached. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cout << "The time limit has been reached. Terminating..." << std::endl;
+        }
         set_reason(TimeLimit);
         terminate();
         return;
     }
 
     if (m_n_iterations > get_param_iteration_limit()) {
-        std::cout << "The iteration limit has been reached. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cout << "The iteration limit has been reached. Terminating..." << std::endl;
+        }
         set_reason(IterLimit);
         terminate();
         return;
     }
 
     if (Algorithm::get_relative_gap() <= get_tol_mip_relative_gap() || Algorithm::get_absolute_gap() <= get_tol_mip_absolute_gap()) {
-        std::cout << "The optimality gap has been closed. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cout << "The optimality gap has been closed. Terminating..." << std::endl;
+        }
         set_status(Optimal);
         set_reason(Proved);
         terminate();
diff --git a/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.cpp b/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.cpp
index cebbd863e..c1499c10b 100644
--- a/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.cpp
+++ b/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Formulation.cpp
@@ -16,6 +16,7 @@ idol::CVCCG::Formulation::Formulation(const idol::Optimizers::Robust::CriticalVa
     check_assumptions();
     initialize_master();
     initialize_sub_problem();
+    m_global_lower_bound = compute_trivial_lower_bound();
 
 }
 
@@ -111,7 +112,7 @@ void idol::CVCCG::Formulation::initialize_master() {
     }
 
     // Create cover constraints
-    if (m_use_cover_constraints) {
+    if (!m_parent.use_indicator() && m_use_cover_constraints) {
         for (auto& linking : m_linking_constraints) {
             linking.cover_constraint = m_master.add_ctr(LinExpr(), LessOrEqual, 1, "__cover_" + linking.ctr_in_uncertainty_set.name());
         }
@@ -392,8 +393,107 @@ double idol::CVCCG::Formulation::compute_critical_value(const Ctr& t_ctr, const
     return result + 1;
 }
 
+double idol::CVCCG::Formulation::compute_trivial_lower_bound() const {
+
+    const auto& model = m_parent.parent();
+    const auto& uncertainty_set = m_parent.description().uncertainty_set();
+    const auto& objective = m_parent.parent().get_obj_expr();
+
+    assert(!objective.has_quadratic());
+
+    double result = objective.affine().constant();
+
+    for (const auto& [var, coeff] : objective.affine().linear()) {
+        if (coeff > 0) {
+            const double lb = model.get_var_lb(var);
+            assert(!is_neg_inf(lb));
+            result += coeff * lb;
+        } else {
+            const double ub = model.get_var_ub(var);
+            assert(!is_pos_inf(ub));
+            result += coeff * ub;
+        }
+    }
+
+    for (const auto& [var, unc_coefff] : m_parent.description().uncertain_objs()) {
+        for (const auto& [unc_var, coeff] : unc_coefff) {
+
+            // Underestimate the following
+            const double lb_var = model.get_var_lb(var);
+            const double ub_var = model.get_var_ub(var);
+            const double lb_unc_var = uncertainty_set.get_var_lb(unc_var);
+            const double ub_unc_var = uncertainty_set.get_var_ub(unc_var);
+
+            const double v1 = lb_var * lb_unc_var;
+            const double v2 = lb_var * ub_unc_var;
+            const double v3 = ub_var * lb_unc_var;
+            const double v4 = ub_var * ub_unc_var;
+
+            const double min_prod = std::min({v1, v2, v3, v4});
+            const double max_prod = std::max({v1, v2, v3, v4});
+
+            if (coeff >= 0) {
+                result += coeff * min_prod;
+            } else {
+                result += coeff * max_prod;
+            }
+
+        }
+    }
+
+    return result;
+}
+
+double idol::CVCCG::Formulation::compute_penalty(const LinExpr<Var, double>& t_row, CtrType t_type, double t_rhs) const {
+
+    // a^\top x \le b => M = max(0, (max_{x\in [l,u]} a^\top x ) - b) ==> a^\top x \le b + Mz
+    // a^T x >= b => M = max(0, b - (min_{x in [l,u]} a^T x)) ==> a^\top x >= b -Mz
+
+    const auto& model = m_parent.parent();
+
+    if (t_type == LessOrEqual) {
+
+        double max_activity = 0;
+
+        for (const auto& [var, coeff] : t_row) {
+            if (coeff > 0) {
+                const double ub = model.get_var_ub(var);
+                assert(!is_pos_inf(ub));
+                max_activity += coeff * ub;
+            } else {
+                const double lb = m_epigraph_variable && m_epigraph_variable->id() == var.id() ? m_master.get_var_lb(var) : model.get_var_lb(var);
+                assert(!is_neg_inf(lb));
+                max_activity += coeff * lb;
+            }
+        }
+
+        return std::max(0., max_activity - t_rhs);
+    }
+
+    double min_activity = 0;
+
+    for (const auto& [var, coeff] : t_row) {
+        if (coeff > 0) {
+            const double lb = model.get_var_lb(var);
+            assert(!is_neg_inf(lb));
+            min_activity += coeff * lb;
+        } else {
+            const double ub = model.get_var_ub(var);
+            assert(!is_pos_inf(ub));
+            min_activity += coeff * ub;
+        }
+    }
+
+    return std::max(0., t_rhs - min_activity);
+
+}
+
 void idol::CVCCG::Formulation::create_critical_value_variable_if_needed(const PrimalPoint& t_scenario) {
 
+    if (m_parent.use_indicator()) {
+        return;
+    }
+
     for (auto& linking : m_linking_constraints) {
 
         const double critical_value = compute_critical_value(linking.ctr_in_uncertainty_set, t_scenario);
@@ -480,66 +580,55 @@ void idol::CVCCG::Formulation::add_scenario_to_master(const std::list<GeneratedS
     const auto& description = m_parent.description();
     const auto& [scenario, master_solution] = *t_iterator_in_pool;
 
-    const LinExpr<Var, LinExpr<Var>>* unc_coeffs = nullptr;
+    const LinExpr<Var, LinExpr<Var>>* unc_row_coeffs = nullptr;
     const LinExpr<Var>* unc_constant = nullptr;
-    CtrType type = LessOrEqual;
+    CtrType type;
 
     LinExpr<Var> row;
     double constant = 0.;
 
     if (t_uncertainty.is_constraint()) {
-        const auto& ctr = t_uncertainty.ctr();
 
+        // We have an uncertain constraint : a^\top x \le b or a^T x >= b
+
+        const auto& ctr = t_uncertainty.ctr();
         row = model.get_ctr_row(ctr);
         type = model.get_ctr_type(ctr);
         constant = model.get_ctr_rhs(ctr);
-
-        unc_coeffs = &description.uncertain_mat_coeffs(ctr);
+        unc_row_coeffs = &description.uncertain_mat_coeffs(ctr);
         unc_constant = &description.uncertain_rhs(ctr);
+
     } else {
-        row = model.get_obj_expr().affine().linear();
-        unc_coeffs = &description.uncertain_objs();
-        constant = model.get_obj_expr().affine().constant();
+
+        // We rewrite the "objective" in epigraph as c^\top x - theta \le -c_0
 
         if (!m_epigraph_variable) {
-            m_epigraph_variable = m_master.add_var(-Inf, Inf, Continuous, 1, "__epigraph");
+            m_epigraph_variable = m_master.add_var(m_global_lower_bound, Inf, Continuous, 1, "__epigraph");
         }
 
+        row = model.get_obj_expr().affine().linear() - *m_epigraph_variable;
+        type = LessOrEqual;
+        constant = -model.get_obj_expr().affine().constant();
+        unc_row_coeffs = &description.uncertain_objs();
+        unc_constant = nullptr;
+
     }
 
-    for (const auto& [var, unc_coeff] : *unc_coeffs) {
-        row += evaluate(unc_coeff, scenario) * var;
+    // Add uncertain coefficients terms
+    if (unc_row_coeffs) {
+        for (const auto& [var, unc_coeff] : *unc_row_coeffs) {
+            row += evaluate(unc_coeff, scenario) * var;
+        }
     }
+
+    // Add uncertain constant terms
     if (unc_constant) {
         for (const auto& [unc_param, coeff] : *unc_constant) {
             constant += coeff * scenario.get(unc_param);
         }
     }
 
-    double penalty = -constant;
-    for (const auto& [var, coeff] : row) {
-        if (coeff > 0) {
-            const double ub = model.get_var_ub(var);
-            if (is_pos_inf(ub)) {
-                throw Exception("Found variable with infinite bound during big-M computation.");
-            }
-            penalty += coeff * ub;
-        } else {
-            const double lb = model.get_var_lb(var);
-            if (is_neg_inf(lb)) {
-                throw Exception("Found variable with infinite bound during big-M computation.");
-            }
-            penalty += coeff * lb;
-        }
-    }
-
-    if (type == LessOrEqual) {
-        penalty *= -1.;
-    }
-
-    if (!t_uncertainty.is_constraint()) {
-        row -= *m_epigraph_variable;
-    }
+    const double penalty = (type == LessOrEqual ? -1. : 1.) * compute_penalty(row, type, constant);
 
     if (uses_indicator()) { // Add indicator
 
diff --git a/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.cpp b/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.cpp
index 3bce001e9..0fc16a7ea 100644
--- a/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.cpp
+++ b/lib/src/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.cpp
@@ -4,6 +4,7 @@
 #include "idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h"
 
 #include "idol/general/optimizers/logs.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/ColumnAndConstraintGeneration.h"
 
 idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::CriticalValueColumnAndConstraintGeneration(const Model& t_model, const idol::Robust::Description& t_description, const OptimizerFactory& t_master_optimizer, const OptimizerFactory& t_deterministic_optimizer, bool t_use_indicator) :
     Algorithm(t_model),
@@ -47,6 +48,17 @@ void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::hook_
         if (is_terminated()) { break; }
 
     } while (true);
+
+    if (get_param_logs()) {
+        if (!m_branching_candidates.empty()) {
+            unsigned int n_fractional = 0;
+            for (const auto& var : m_branching_candidates) {
+                const double val = get_var_primal(var);
+                n_fractional += !is_integer(val, get_tol_integer());
+            }
+            std::cout << "Solution diversity: " << (100.0 * n_fractional / static_cast<double>(m_branching_candidates.size())) << " %." << std::endl;
+        }
+    }
 }
 
 void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::solve_master_problem() {
@@ -81,7 +93,9 @@ void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::check
     }
 
     if (get_best_bound() > get_best_obj() + get_tol_mip_absolute_gap()) {
-        std::cerr << "The current best bound is larger than current best obj. This should should not happen. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cerr << "The current best bound is larger than current best obj. This should should not happen. Terminating..." << std::endl;
+        }
         set_status(Fail);
         set_reason(Numerical);
         terminate();
@@ -91,21 +105,27 @@ void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::check
     set_status(is_pos_inf(get_best_obj()) ? Infeasible : Feasible);
 
     if (get_remaining_time() == 0) {
-        std::cout << "The time limit has been reached. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cout << "The time limit has been reached. Terminating..." << std::endl;
+        }
         set_reason(TimeLimit);
         terminate();
         return;
     }
 
     if (m_n_iterations > get_param_iteration_limit()) {
-        std::cout << "The iteration limit has been reached. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cout << "The iteration limit has been reached. Terminating..." << std::endl;
+        }
         set_reason(IterLimit);
         terminate();
         return;
     }
 
     if (Algorithm::get_relative_gap() <= get_tol_mip_relative_gap() || Algorithm::get_absolute_gap() <= get_tol_mip_absolute_gap()) {
-        std::cout << "The optimality gap has been closed. Terminating..." << std::endl;
+        if (get_param_logs()) {
+            std::cout << "The optimality gap has been closed. Terminating..." << std::endl;
+        }
         set_status(Optimal);
         set_reason(Proved);
         terminate();
@@ -130,25 +150,32 @@ void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::solve
         auto& sub_problem = m_formulation->sub_problem();
         sub_problem.optimize();
 
-        auto scenario = save_primal(sub_problem);
+        const auto status = sub_problem.get_status();
 
-        if (scenario.status() != Optimal) {
-            set_status(scenario.status());
-            set_reason(scenario.reason());
+        if (status != Optimal) {
+            set_status(status);
+            set_reason(sub_problem.get_reason());
             terminate();
             return;
         }
 
+        auto scenario = save_primal(sub_problem);
+
         if (!uncertainty.is_constraint()) {
+
             const double LB = get_best_bound();
             const double UB = LB - scenario.objective_value();
             if (relative_gap(LB, UB) > get_tol_mip_absolute_gap() && absolute_gap(LB, UB) > get_tol_mip_absolute_gap()) {
                 scenarios.emplace_back(std::move(scenario));
-                is_feasible = false;
             }
-        } else if (!m_formulation->master_provides_a_valid_bound() || -scenario.objective_value() > get_tol_feasibility()) {
+
+            continue;
+        }
+
+        if (!m_formulation->master_provides_a_valid_bound() || -scenario.objective_value() > get_tol_feasibility()) {
             scenarios.emplace_back(std::move(scenario));
             is_feasible = false;
+            continue;
         }
 
     }
@@ -163,7 +190,7 @@ void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::solve
 
     if (is_feasible && m_formulation->master_provides_a_valid_bound()) {
         assert(scenarios.size() == 1);
-        set_best_obj(master_solution.get(m_formulation->epigraph_variable()) - scenarios.front().objective_value());
+        set_best_obj(std::min(get_best_obj(), master_solution.get(m_formulation->epigraph_variable()) - scenarios.front().objective_value()));
     }
 
     if (scenarios.size() != 1) {
@@ -238,6 +265,7 @@ void idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration::hook_
     Algorithm::hook_before_optimize();
 
     m_n_iterations = 0;
+    m_is_diving = false;
     set_best_obj(Inf);
     set_best_bound(-Inf);