Branch-n-price implementation complications

import gurobipy as gu
import math
import time

class Problem:
    def __init__(self, dfData, Cover):
        self.I = dfData['I'].dropna().astype(int).unique().tolist()
        self.T = dfData['T'].dropna().astype(int).unique().tolist()
        self.K = dfData['K'].dropna().astype(int).unique().tolist()
        self.cover = Cover
        self.model = gu.Model("MasterProblem")
        self.Min_WD = 3
        self.Max_WD = 5
        self.alpha = {1: 0.87, 2: 0.83, 3: 1, 4: 0.91, 5: 0.88, 6: 1, 7: 1, 8: 0.91, 9: 1, 10: 0.84, 11: 0.82, 12: 0.92, 13: 1, 14: 0.95}
        self.M = 100
        self.F_S = [(3, 1), (3, 2), (2, 1)]

    def buildLinModel(self):
        self.t0 = time.time()
        self.generateVariables()
        self.genGenCons()
        self.generateObjective()
        self.model.update()


    def generateVariables(self):
        self.x = self.model.addVars(self.I, self.T, self.K, vtype=gu.GRB.BINARY, name="x")
        self.y = self.model.addVars(self.I, self.T, vtype=gu.GRB.BINARY, name="y")
        self.u = self.model.addVars(self.T, self.K, lb = 0, vtype=gu.GRB.CONTINUOUS, name="u")
        self.moti = self.model.addVars(self.I, self.T, self.K, vtype=gu.GRB.CONTINUOUS, lb=0, ub=1, name="perf")
        self.mood = self.model.addVars(self.I, self.T, vtype=gu.GRB.BINARY, name="mood")


    def genGenCons(self):
        for i in self.I:
            for t in range(1, len(self.T) - self.Max_WD + 1):
                self.model.addLConstr(
                    gu.quicksum(self.y[i, u] for u in range(t, t + 1 + self.Max_WD)) <= self.Max_WD)
            for t in range(2, len(self.T) - self.Min_WD + 1):
                self.model.addLConstr(
                    gu.quicksum(self.y[i, u] for u in range(t + 1, t + self.Min_WD + 1)) >= self.Min_WD * (
                                self.y[i, t + 1] - self.y[i, t]))
            for k1, k2 in self.F_S:
                for t in range(1, len(self.T)):
                    self.model.addLConstr(self.x[i, t, k1] + self.x[i, t + 1, k2] <= 1)
            for t in self.T:
                self.model.addLConstr(gu.quicksum(self.x[i, t, k] for k in self.K) == self.y[i, t])
                for k in self.K:
                    self.model.addLConstr(self.moti[i, t, k] == self.alpha[t] * self.x[i, t, k], name=f"moti_{t}_{k}")
        for t in self.T:
            for k in self.K:
                self.model.addLConstr(
                    gu.quicksum(self.x[i, t, k] for i in self.I) + self.u[t, k] >= self.cover[t, k])
        #for i in self.I:
            #for t in self.T:
                #self.model.addConstr(self.alpha[t] * gu.quicksum(self.x[i, t, k] for k in self.K) + self.mood[i, t] == 1, name="Mood_constraint")
                #for k in self.K:
                    #self.model.addConstr(-self.M * (1 - self.x[i, t, k]) <= self.moti[i, t, k] - self.mood[i, t], name="Motivation_lower_bound")
                    #self.model.addConstr(self.moti[i, t, k] - self.mood[i, t] <= self.M * (1 - self.x[i, t, k]), name="Motivation_upper_bound")
                    #self.model.addConstr(self.moti[i, t, k] <= self.x[i, t, k], name="Motivation_upper_limit")
        self.model.update()

    def generateObjective(self):
        self.model.setObjective(gu.quicksum(self.u[t, k] for k in self.K for t in self.T), sense=gu.GRB.MINIMIZE)

    def solveModel(self):
        self.model.optimize()

import gurobipy as gu
import statistics
import numpy as np
import math


class MasterProblem:
    def __init__(self, df, Coverage, max_iteration, current_iteration, start):
        self.max_iteration = max_iteration
        self.nurses = df['I'].dropna().astype(int).unique().tolist()
        self.days = df['T'].dropna().astype(int).unique().tolist()
        self.shifts = df['K'].dropna().astype(int).unique().tolist()
        self._current_iteration = current_iteration
        self.roster = [i for i in range(1, self.max_iteration + 2)]
        self.cover = Coverage
        self.model = gu.Model("MasterProblem")
        self.cons_cover = {}
        self.cons_lmbda = {}
        self.cover_values = [self.cover[key] for key in self.cover.keys()]
        self.start = start
        self.all_schedules = {}
        self.branch_constraints = []
        self.branch_constraints_mp = []
        self.branch_history = {}

    def buildModel(self):
        self.generateVariables()
        self.generateConstraints()
        self.generateObjective()
        self.model.update()

    def generateVariables(self):
        self.u = self.model.addVars(self.days, self.shifts, vtype=gu.GRB.CONTINUOUS, lb=0, name='u')
        self.lmbda = self.model.addVars(self.roster, vtype=gu.GRB.INTEGER, lb=0, name='lmbda')

    def generateConstraints(self):
        self.cons_lmbda = self.model.addLConstr(len(self.nurses) == gu.quicksum(self.lmbda[r] for r in self.roster), name="lmb")
        for t in self.days:
            for s in self.shifts:
                self.cons_cover[t, s] = self.model.addConstr(gu.quicksum(self.lmbda[r] for r in self.roster) + self.u[t, s] >= self.cover[t, s], "cover(" + str(t) + "," + str(s) + ")")
        return self.cons_lmbda, self.cons_cover

    def generateObjective(self):
        self.model.setObjective(gu.quicksum(self.u[t, s] for t in self.days for s in self.shifts), sense=gu.GRB.MINIMIZE)

    def getDuals(self):
        dual_sum_left = 0
        dual_sum_right = 0
        formula_parts_left = []
        formula_parts_right = []

        for constr in self.model.getConstrs():
            if 'left' in constr.ConstrName:
                value = constr.Pi
                formula_parts_left.append(f"{value} + ")
                print(f'Left {constr, value}')
                dual_sum_left += value*(-1)
            elif 'right' in constr.ConstrName:
                value = constr.Pi
                formula_parts_right.append(f"{value} - ")
                print(f'Right {constr, value}')
                dual_sum_right += value
            else:
                dual_sum_left, dual_sum_right = 0, 0
        print(f'Left, Right: {dual_sum_left, dual_sum_right}')
        formula = " ".join(formula_parts_left) + " - " + " ".join(formula_parts_right)
        print(f'Dual formula {formula}')
        dual_sum_total = dual_sum_left + dual_sum_right

        print(dual_sum_total)
        return {(d, s): self.cons_cover[d, s].Pi for d in self.days for s in self.shifts}, self.cons_lmbda.Pi, dual_sum_left, dual_sum_right

    def startSol(self, start_vals):
        for t in self.days:
            for s in self.shifts:
                if (t, s) in start_vals:
                    value_cons = start_vals[t, s]
                else:
                    value_cons = 0

                if (t, s) in start_vals:
                    self.model.chgCoeff(self.cons_cover[t, s], self.lmbda[1], value_cons)
        self.model.update()

    def initCoeffs(self):
        for t in self.days:
            for s in self.shifts:
                for r in self.roster:
                    self.model.chgCoeff(self.cons_cover[t, s], self.lmbda[r], 0)
        self.model.update()

    def addCol(self, itr, schedules_perf, index=None):
        for t in self.days:
            for s in self.shifts:
                if (t, s, itr + 1) in schedules_perf:
                    value_cons = schedules_perf[t, s, itr + 1]
                else:
                    value_cons = 0

                if (t, s, itr + 1) in schedules_perf:
                    self.model.chgCoeff(self.cons_cover[t, s], self.lmbda[itr + 1], value_cons)
        self.model.update()

    def finalSolve(self, timeLimit):
        self.model.setParam('TimeLimit', timeLimit)
        self.model.Params.OutputFlag = 1
        self.model.setAttr("vType", self.lmbda, gu.GRB.INTEGER)
        self.model.update()
        self.model.optimize()
        if self.model.status == gu.GRB.OPTIMAL:
            print("***** Optimal solution found *****")
        else:
            print("***** No optimal solution found *****")


    def solveModel(self, timeLimit):
        try:
            self.model.setParam('TimeLimit', timeLimit)
            self.model.Params.OutputFlag = 0
            self.model.Params.IntegralityFocus = 1
            self.model.Params.MIPGap = 1e-5
            self.model.setParam('ConcurrentMIP', 2)
            self.model.optimize()
        except gu.GurobiError as e:
            print('Error code ' + str(e.errno) + ': ' + str(e))

    def solveRelaxModel(self):
        try:
            self.model.Params.OutputFlag = 0
            self.model.Params.MIPGap = 1e-6
            self.model.Params.Method = 2
            self.model.Params.Crossover = 0
            for v in self.model.getVars():
                v.setAttr('vtype', 'C')
                v.setAttr('lb', 0.0)
            self.model.optimize()
        except gu.GurobiError as e:
            print('Error code ' + str(e.errno) + ': ' + str(e))

    def addSchedule(self, schedule):
        self.all_schedules.update(schedule)

    def exportSchedules(self):
        return self.all_schedules


    def addColsFromDict(self, all_schedules):
        for (t, s, it), value_cons in all_schedules.items():
            if (t, s) in self.cons_cover:
                self.model.chgCoeff(self.cons_cover[t, s], self.lmbda[it], value_cons)
        self.model.update()


    def add_branch_constraint(self, day, shift, direction, bound_value, schedules):
        """
        Adds a branching constraint to the model and tracks it in branch_constraints.

        Args:
            day (int): Der Tag der Constraint.
            shift (int): Die Schicht der Constraint.
            direction (str): 'left' für <= oder 'right' für >=.
            bound_value (float): Der Wert, der als Grenze genutzt wird.
            relevant (list): Liste der relevanten Variablen.
        """



        expr = gu.quicksum(self.lmbda[r] * schedules.get((day, shift, r), 0) for r in self.roster)


        if direction == 'left':
            constr = self.model.addConstr(expr <= math.floor(bound_value),
                                          name=f"branch_left_{day}_{shift}")
        else:
            constr = self.model.addConstr(expr >= math.ceil(bound_value),
                                          name=f"branch_right_{day}_{shift}")

        self.branch_constraints.append({
            'day': day,
            'shift': shift,
            'direction': direction,
            'bound': bound_value,
            'relevant': schedules,
            'constraint': constr
        })

        self.model.update()

    def propagate_constraints(self, parent_problem):
        """Propagates branching constraints from parent node"""
        for constraint_info in parent_problem.branch_constraints:
            self.add_branch_constraint(
                constraint_info['day'],
                constraint_info['shift'],
                constraint_info['direction'],
                constraint_info['bound'],
                constraint_info['relevant']
            )

        self.model.update()

    def add_branch_constraint_mp(self, ind, bound_value, direction):
        """
        Adds a branching constraint to the model and tracks it in branch_constraints.

        Args:
            day (int): Der Tag der Constraint.
            shift (int): Die Schicht der Constraint.
            direction (str): 'left' für <= oder 'right' für >=.
            bound_value (float): Der Wert, der als Grenze genutzt wird.
            relevant (list): Liste der relevanten Variablen.
        """


        # Neues Constraint hinzufügen
        if direction == 'left':
            constr = self.model.addConstr(-1 * self.lmbda[ind] >= -1 * math.floor(bound_value),
                                          name=f"branch_left_{ind}")
        else:
            constr = self.model.addConstr(self.lmbda[ind] >= math.ceil(bound_value),
                                          name=f"branch_right_{ind}")

        # Neues Constraint zu `branch_constraints` hinzufügen
        self.branch_constraints.append({
            'day': None,
            'shift': None,
            'direction': direction,
            'bound': bound_value,
            'relevant': None,
            'constraint': constr
        })

        self.model.update()

    def find_and_print_fractional_betas(self, schedule, lmb_vals):

        beta_values = {}
        for d in self.days:
            for s in self.shifts:
                beta_d_s = sum(schedule.get((d, s, k), 0) * lmb_vals.get((k), 0) for k in self.roster)
                beta_values[d, s] = beta_d_s
                print(f'Beta{d}{s} = {beta_d_s}')

        fractional_betas = {}
        for (d, s), beta in beta_values.items():
            if abs(beta - round(beta)) > 1e-6:
                fractional_betas[(d, s)] = beta

        if not fractional_betas:
            None
        else:
            for (d, s), beta in fractional_betas.items():
                fractional_part = abs(beta - round(beta))
                print(f"beta_{d}_{s} = {beta}, Fraktionsanteil = {fractional_part}")

        if not fractional_betas:
            most_fractional_beta = None
        else:
            most_fractional_beta = max(fractional_betas.items(), key=lambda x: abs(x[1] - round(x[1])))
            print(
                f"beta_{most_fractional_beta[0][0]}_{most_fractional_beta[0][1]} = {most_fractional_beta[1]}")

        k_values = []
        for (key_t, key_s, key_k), value in schedule.items():
            if key_t == most_fractional_beta[0][0] and key_s == most_fractional_beta[0][1] and value > 0:
                k_values.append(key_k)

        return most_fractional_beta[0][0], most_fractional_beta[0][1], most_fractional_beta[1], [k for k in k_values if lmb_vals.get(k, 0) > 0]

    def getLambdas(self):
        return self.model.getAttr("X", self.lmbda)

import gurobipy as gu
import math

class Subproblem:
    def __init__(self, duals_i, duals_ts, df, iteration, branch_duals_left, branch_duals_right):
        itr = iteration + 1
        self.days = df['T'].dropna().astype(int).unique().tolist()
        self.shifts = df['K'].dropna().astype(int).unique().tolist()
        self.duals_i = duals_i
        self.duals_ts = duals_ts
        self.duals_branch_left = branch_duals_left
        self.duals_branch_right = branch_duals_right
        self.model = gu.Model("Subproblem")
        self.itr = itr
        self.F_S = [(3, 1), (3, 2), (2, 1)]
        self.Min_WD = 3
        self.Max_WD = 5
        self.M = 100
        self.alpha = {1: 0.87, 2: 0.83, 3: 0.92, 4: 0.91, 5: 0.88, 6: 1, 7: 0.93, 8: 0.91, 9: 0.93, 10: 0.84, 11: 0.82, 12: 0.92, 13: 1, 14: 0.95}
        self.branch_constraints = []
        self.branch_constraints_mp = []


    def buildModel(self):
        self.generateVariables()
        self.generateConstraints()
        self.generateObjective()
        self.model.update()

    def generateVariables(self):
        self.x = self.model.addVars(self.days, self.shifts, vtype=gu.GRB.BINARY, name="x")
        self.y = self.model.addVars(self.days, vtype=gu.GRB.BINARY, name="y")
        self.mood = self.model.addVars(self.days, vtype=gu.GRB.BINARY, name="mood")
        self.moti = self.model.addVars(self.days, self.shifts, [self.itr], vtype=gu.GRB.CONTINUOUS, lb=0, ub=1, name="moti")

    def generateConstraints(self):
        for t in self.days:
            self.model.addLConstr(gu.quicksum(self.x[t, k] for k in self.shifts) == self.y[t], name = f"day_assignment_{t}")
            for s in self.shifts:
                self.model.addConstr(self.alpha[t] * gu.quicksum(self.x[t, s] for s in self.shifts)  == self.moti[t, s, self.itr], name="Mood_constraint")
        for k1, k2 in self.F_S:
            for t in range(1, len(self.days)):
                self.model.addLConstr(self.x[t, k1] + self.x[t + 1, k2] <= 1, name = f"rotation_{t}")
        for t in range(1, len(self.days) - self.Max_WD + 1):
            self.model.addLConstr(gu.quicksum(self.y[u] for u in range(t, t + 1 + self.Max_WD)) <= self.Max_WD, name = f"max_days_{t}")
        for t in range(2, len(self.days) - self.Min_WD + 1):
            self.model.addLConstr(gu.quicksum(self.y[u] for u in range(t + 1, t + self.Min_WD + 1)) >= self.Min_WD * (self.y[t + 1] - self.y[t]), name = f"min_days_{t}")
        self.model.update()

    def generateObjective(self):
        self.model.setObjective(0 - gu.quicksum(self.moti[t, s, self.itr] * self.duals_ts[t, s] for t in self.days for s in self.shifts) - self.duals_i + self.duals_branch_left - self.duals_branch_right, sense=gu.GRB.MINIMIZE)

    def solveSP(self):
        self.model.Params.OutputFlag = 1
        self.model.optimize()
        if self.model.status != gu.GRB.OPTIMAL:
            self.model.computeIIS()
            print('\nThe following constraints and variables are in the IIS:')
            for c in self.model.getConstrs():
                if c.IISConstr: print(f'\t{c.constrname}: {self.model.getRow(c)} {c.Sense} {c.RHS}')

            for v in self.model.getVars():
                if v.IISLB: print(f'\t{v.varname} ≥ {v.LB}')
                if v.IISUB: print(f'\t{v.varname} ≤ {v.UB}')
            else:
                None


    def add_branch_constraint(self, day, shift, direction):
        """Adds branbetang constraint to subproblem."""

        if direction == 'left':
            self.model.addLConstr(0 == self.x[day, shift], name=f"sp_branch_left_{day}_{shift}")
        else:
            self.model.addConstr(1 == self.x[day, shift], name=f"sp_branch_right_{day}_{shift}")

        self.branch_constraints.append({
            'day': day,
            'shift': shift,
            'direction': direction
        })

        self.model.update()

    def propagate_constraints(self, parent_subproblem):
        """Propagates all branbetang constraints from parent"""
        for constraint_info in parent_subproblem.branch_constraints:
            self.add_branch_constraint(
                constraint_info['day'],
                constraint_info['shift'],
                constraint_info['direction']
            )


import pandas as pd

# **** Prerequisites ****
I, T, K = list(range(1,26)), list(range(1,15)), list(range(1,4))
data = pd.DataFrame({
    'I': I + [np.nan] * (max(len(I), len(T), len(K)) - len(I)),
    'T': T + [np.nan] * (max(len(I), len(T), len(K)) - len(T)),
    'K': K + [np.nan] * (max(len(I), len(T), len(K)) - len(K))
})

# Parameter
time_Limit, max_itr, output_len, mue, threshold = (4, 100, 100, 1e-4, 5e-4)
demand_dict = {(1, 1): 9, (1, 2): 11, (1, 3): 4, (2, 1): 10, (2, 2): 9, (2, 3): 9, (3, 1): 8, (3, 2): 16, (3, 3): 2, (4, 1): 4, (4, 2): 21, (4, 3): 3, (5, 1): 2, (5, 2): 7, (5, 3): 13, (6, 1): 7, (6, 2): 5, (6, 3): 13, (7, 1): 14, (7, 2): 1, (7, 3): 5, (8, 1): 7, (8, 2): 12, (8, 3): 3, (9, 1): 3, (9, 2): 11, (9, 3): 8, (10, 1): 1, (10, 2): 13, (10, 3): 6, (11, 1): 19, (11, 2): 8, (11, 3): 4, (12, 1): 11, (12, 2): 4, (12, 3): 15, (13, 1): 6, (13, 2): 6, (13, 3): 12, (14, 1): 4, (14, 2): 10, (14, 3): 9}
dir, dir1 = 'right', 'right'

# **** Compact Solver ****
problem = Problem(data, demand_dict)
problem.buildLinModel()
problem.solveModel()

# **** Column Generation ****
modelImprovable, reached_max_itr = True, False

# Get Starting Solutions
problem_start = Problem(data, demand_dict)
problem_start.buildLinModel()
problem_start.model.Params.MIPGap = 0.8
problem_start.model.update()
problem_start.model.optimize()
start_values_moti = {(t, s): problem_start.x[1, t, s].x for t in T for s in K}
print(start_values_moti)
while True:
    # Initialize iterations
    itr = 0
    last_itr = 0

    # Create empty results lists
    master = MasterProblem(data, demand_dict, max_itr, itr, start_values_moti)
    master.buildModel()
    master.initCoeffs()
    print("*{:^{output_len}}*".format("Restricted Master Problem successfully built!", output_len=output_len))

    # Initialize and solve relaxed model
    master.startSol(start_values_moti)
    master.addSchedule({(t, s, 1): problem_start.moti[1, t, s].x for t in T for s in K})
    master.solveRelaxModel()

    while (modelImprovable) and itr < max_itr:
        # Start
        itr += 1

        # Solve RMP and get duals
        master.current_iteration = itr + 1
        master.solveRelaxModel()
        duals_ts, duals_i, duals_bl, duals_br = master.getDuals()

        # Solve SPs
        modelImprovable = False
        subproblem = Subproblem(duals_i, duals_ts, data, itr, duals_bl, duals_br)
        subproblem.buildModel()
        subproblem.solveSP()
        last_itr = itr + 1

        # Generate and add columns with reduced cost
        if subproblem.model.objval < -threshold:
            Schedules = subproblem.model.getAttr("X", subproblem.moti)
            master.addCol(itr, Schedules)
            master.addSchedule(Schedules)
            master.model.update()
            modelImprovable = True

        if not modelImprovable:
            day, shift, value, filtered_list = master.find_and_print_fractional_betas(master.exportSchedules(), master.getLambdas())
            print(f'Branching Info: {day, shift, value, filtered_list}')
            print("*{:^{output_len}}*".format("No more improvable columns found.", output_len=output_len))
            break

    if modelImprovable and itr == max_itr:
        print("*{:^{output_len}}*".format("More iterations needed. Increase max_itr and restart the process.", output_len=output_len))
        max_itr *= 2
    else:
        break

# Solve Master Problem with integrality restored
master.finalSolve(time_Limit)

print('First Branch - Left')

modelImprovable1, reached_max_itr1 = True, False
while True:
    # Initialize iterations
    itr1, last_itr1 = last_itr - 1, 0


    # Create empty results lists
    master1 = MasterProblem(data, demand_dict, max_itr, itr1, None)
    master1.buildModel()
    master1.initCoeffs()
    master1.addColsFromDict(master.exportSchedules())
    master1.addSchedule(master.exportSchedules())
    master1.propagate_constraints(master)
    master1.add_branch_constraint(day, shift, dir, value, master.exportSchedules())
    master1.model.write('master.lp')

    print("*{:^{output_len}}*".format("Restricted Master Problem successfully built!", output_len=output_len))

    # Initialize and solve relaxed model
    master1.solveRelaxModel()

    while (modelImprovable1) and itr1 < max_itr:
        # Start
        itr1 += 1

        # Solve RMP and get duals
        master1.current_iteration = itr1 + 1
        master1.solveRelaxModel()
        duals_ts1, duals_i1, duals_bl1, duals_br1 = master1.getDuals()

        # Solve SPs
        modelImprovable1 = False
        print("*{:^{output_len}}*".format(f"Begin Column Generation Iteration {itr1}", output_len=output_len))

        # Build SP
        subproblem1 = Subproblem(duals_i1, duals_ts1, data, itr1, duals_bl1, duals_br1)
        subproblem1.propagate_constraints(subproblem)
        subproblem1.buildModel()
        subproblem1.add_branch_constraint(day, shift, dir)
        subproblem1.solveSP()
        redCost1 = subproblem1.model.objval
        last_itr1 = itr1 + 1

        # Generate and add columns with reduced cost
        if redCost1 < -threshold:
            Schedules1 = subproblem1.model.getAttr("X", subproblem1.moti)
            master1.addCol(itr1, Schedules1)
            master1.addSchedule(Schedules1)
            master1.model.update()
            modelImprovable1 = True

        if not modelImprovable1:
            day1, shift1, value1, filtered_list1 = master1.find_and_print_fractional_betas(master1.exportSchedules(), master1.getLambdas())
            print(f'Branching Info1: {day1, shift1, value1, filtered_list1}')
            print("*{:^{output_len}}*".format("No more improvable columns found.", output_len=output_len))
            break

    if modelImprovable1 and itr1 == max_itr:
        break
    else:
        break

# Solve Master Problem with integrality restored
master1.finalSolve(time_Limit)

print('Second Branch - Left')

modelImprovable2, reached_max_itr2 = True, False
while True:
    # Initialize iterations
    itr2, last_itr2 = last_itr1 - 1, 0

    # Create empty results lists
    master2 = MasterProblem(data, demand_dict, max_itr, itr2, None)
    master2.buildModel()
    master2.initCoeffs()
    master2.addColsFromDict(master1.exportSchedules())
    master2.addSchedule(master1.exportSchedules())
    master2.propagate_constraints(master1)
    master2.add_branch_constraint(day1, shift1, dir1, value1, master1.exportSchedules())

    print("*{:^{output_len}}*".format("Restricted Master Problem successfully built!", output_len=output_len))

    # Initialize and solve relaxed model
    master2.solveRelaxModel()

    # Retrieve dual values

    while (modelImprovable2) and itr2 < max_itr:
        # Start
        itr2 += 1

        # Solve RMP and get duals
        master2.current_iteration = itr2 + 1
        master2.solveRelaxModel()
        duals_ts2, duals_i2, duals_bl2, duals_br2 = master2.getDuals()


        # Solve SPs
        modelImprovable2 = False
        print("*{:^{output_len}}*".format(f"Begin Column Generation Iteration {itr2}", output_len=output_len))
        subproblem2 = Subproblem(duals_i2, duals_ts2, data, itr2, duals_bl2, duals_br2)
        subproblem2.buildModel()
        subproblem2.propagate_constraints(subproblem1)
        subproblem2.add_branch_constraint(day1, shift1, dir1)
        subproblem2.model.update()

        # Save time to solve SP
        subproblem2.solveSP()
        last_itr2 = itr2 + 1

        # Generate and add columns with reduced cost
        if subproblem2.model.objval < -threshold:
            Schedules2 = subproblem2.model.getAttr("X", subproblem2.moti)
            master2.addCol(itr2, Schedules2)
            print(f'New columdn in {itr2}: {Schedules2}')
            master2.addSchedule(Schedules2)
            master2.model.update()
            modelImprovable2 = True

        if not modelImprovable2:
            day2, shift2, value2, filtered_list2 = master2.find_and_print_fractional_betas(master2.exportSchedules(), master2.getLambdas())
            print(f'Branching Info2: {day2, shift2, value2, filtered_list2}')
            print("*{:^{output_len}}*".format("No more improvable columns found.", output_len=output_len))
            break

    if modelImprovable2 and itr2 == max_itr:
        break
    else:
        break

# Solve Master Problem with integrality restored
master2.finalSolve(time_Limit)