cutting stock 列生成 python实现

import numpy as np
from docplex.mp.model import Model as CplexModel


class ExactSolution:
    def __init__(self, mode_name, mode_corporate="cplex"):
        if mode_corporate == "cplex":
            self.model = CplexModel(mode_name)
        # elif mode_corporate == "scip":
        #     pass
        # todo 待补充
        else:
            raise AttributeError("the mode_corporate is not correct!")

    def __call__(self):
        self.add_decision_var()
        self.add_constraint()
        self.set_obj()
        # start solve
        self.optimize()
        return None

    def add_decision_var(self):
        raise NotImplementedError

    def add_constraint(self):
        raise NotImplementedError

    def set_obj(self, *args):
        raise NotImplementedError

    def optimize(self):
        raise NotImplementedError


class MasterProblem(ExactSolution):
    def __init__(self, model_name, lengths, quantities, max_length, mode_corporate="cplex"):
        self.lengths = lengths
        self.quantities = quantities
        self.max_length = max_length
        self.constr_coeff = np.eye(len(self.lengths)).astype(int) * [
            self.max_length // len_ for len_ in self.lengths]
        super().__init__(model_name, mode_corporate)

    def add_decision_var(self):
        self.x_var = {}
        for ix in range(len(self.lengths)):
            self.x_var[ix] = self.model.continuous_var(lb=0, name="x_%s" % ix)
        self.cut_method_count = len(self.lengths)

    def add_constraint(self):
        self.model.add_constraints(self.model.sum(self.x_var[j] * self.constr_coeff[i, j]
                                                  for j in range(self.cut_method_count)) >= self.quantities[i]
                                   for i in range(len(self.lengths)))

    def set_obj(self):
        self.model.minimize(self.model.sum(self.x_var[i] for i in range(self.cut_method_count)))

    def optimize(self):
        self.model.solve(log_output=False)
        if not self.model.solution:
            raise RuntimeError("can't get solution!")

    def get_dual_vars(self):
        self.model.solution.ensure_dual_values(self.model, self.model.get_engine())
        return self.model.solution.get_dual_values(self.model.iter_constraints())

    def update_model_and_solve(self, new_add_y):
        self.update_model_add_var()
        self.update_model_mdf_constraint(new_add_y)
        # self.model.clean_before_solve = True
        self.set_obj()
        self.optimize()

    def update_model_add_var(self):
        self.x_var[self.cut_method_count] = self.model.continuous_var(lb=0, name="addx_%s" % self.cut_method_count)
        self.cut_method_count += 1

    def update_model_mdf_constraint(self, new_col):
        self.constr_coeff = np.hstack((self.constr_coeff, np.array([new_col]).astype(int).transpose()))
        self.model.clear_constraints()
        self.add_constraint()


class SubProblem(ExactSolution):
    def __init__(self, model_name, lengths, quantities, max_length, mode_corporate="cplex"):
        self.lengths = lengths
        self.quantities = quantities
        self.max_length = max_length
        super().__init__(model_name, mode_corporate)
        self.add_decision_var()
        self.add_constraint()

    def add_decision_var(self):
        self.x_var = {}
        for ix in range(len(self.lengths)):
            self.x_var[ix] = self.model.integer_var(lb=0, name="y_%s" % ix)

    def add_constraint(self):
        self.model.add_constraint(
            self.model.sum(self.lengths[i] * self.x_var[i] for i in range(len(self.lengths))) <= self.max_length)

    def set_obj(self, multiplies):
        # self.model.maximize(self.model.sum(multiplies[i] * self.x_var[i] for i in range(len(self.lengths))))
        self.model.minimize(1 - self.model.sum(multiplies[i] * self.x_var[i] for i in range(len(self.lengths))))

    def optimize(self):
        self.model.solve(log_output=False)
        if not self.model.solution:
            raise RuntimeError("can't get solution!")

    def __call__(self, *args, **kwargs):
        # cutting stock的subproblem有点特殊 需要单独设置obj后才能跑
        raise RuntimeError("please call by other method.")

    def get_solution(self):
        return [self.model.get_var_by_index(i).solution_value for i in range(len(self.lengths))]

    def get_reduced_cost(self):
        return self.model.objective_value


class CuttingStock:
    def __init__(self):
        pass

    def reset(self):
        self.reduced_cost = float('-inf')

    def run(self, lengths, quantities, max_length, master_model_name='master', sub_model_name='sub'):
        self.reset()
        self.get_linear_result(lengths, quantities, max_length, master_model_name, sub_model_name)
        self.get_round_result()

    def get_linear_result(self, lengths, quantities, max_length, master_model_name, sub_model_name):
        master_prob = MasterProblem(master_model_name, lengths, quantities, max_length)
        master_prob()
        sub_prob = SubProblem(sub_model_name, lengths, quantities, max_length)

        while self.reduced_cost < 0:
            pi = master_prob.get_dual_vars()
            sub_prob.set_obj(pi)
            sub_prob.optimize()
            y = sub_prob.get_solution()
            self.reduced_cost = sub_prob.get_reduced_cost()
            print("the reduce_cost is %s, the obj_value is %s "
                  % (self.reduced_cost, master_prob.model.solution.objective_value))
            master_prob.update_model_and_solve(y)
        # print("final result:   ", list(master_prob.model.iter_variables()))

    def get_round_result(self):
        # todo
        pass


if __name__ == "__main__":
    # max_length = 20  # width of large roll
    # lengths = [3, 7, 9, 16]
    # quantities = [25, 30, 14, 8]
    max_length = 16  # width of large roll
    lengths = [3, 6, 7]
    quantities = [25, 20, 18]
    CuttingStock().run(lengths, quantities, max_length)

这个结果还是线性的 还需要转成整数的处理

另外， 如果需求是3种长度，那结果一定是3种砍法方案的组合； 需求是4种，结果也是4种。好像觉得这个建模有点不准确？极端情况， 如果需要 3，6，7  三种长度的各20根，总长为16，解出来总是不对的..