# Cutting stock model with column generation. # # (C) 2020-2025 Fair Isaac Corporation import xpress as xp import numpy as np # Problem data: # widths = [17.0, 21.0, 22.5, 24.0, 29,5] # demand = [150, 96, 48, 108, 227] # max_width = 94 widths = [6, 11, 17, 21, 24, 28, 30, 33, 42, 49, 56, 69, 74, 87, 91] demand = [9, 6, 20, 30, 17, 19, 25, 12, 8, 20, 5, 14, 15, 18, 10] max_width = 100 eps = 1e-6 max_cols = 100 def solve_knapsack_problem(n_item_types, profit, resource_consumed, total_resource, demand): """ Solve knapsack subproblem. :param n_item_types: Total number of items :param profit: Profit of each item. 1-d numpy array :param resource_consumed: Resource consumed by each item. 1-d numpy array :param total_resource: total available resource :param demand: Demand for each item. 1-d numpy array :return: xbest: Number of items of type i used in an optimal solution zbest: Value of the optimal solution """ N = range(n_item_types) p = xp.problem(name="knapsack") p.controls.outputlog = 0 # Create a NumPy array of variables x = xp.array([p.addVariable(name='use_{}'.format(i), lb=0, ub=demand[i], vartype=xp.integer) for i in N]) p.addConstraint(xp.Dot(resource_consumed, x) <= total_resource) p.setObjective(xp.Dot(profit, x), sense=xp.maximize) p.optimize() return p.attributes.objval, p.getSolution() def cutting_stock_model(): """ Solve a cutting stock problem by defining main problem and subproblem. """ n_widths = len(widths) N = range(n_widths) p = xp.problem(name='CutStock') patterns = np.zeros((n_widths, n_widths)) for j in N: patterns[j][j] = int(max_width / widths[j]) # Create array of variables v_patterns = xp.array([ p.addVariable(name='pat_{}'.format(i), lb=0, ub=int(float(demand[i] / patterns[i][i]) + 1)) for i in N ]) p.setObjective(xp.Dot(np.ones(n_widths), v_patterns), sense=xp.minimize) p.addConstraint(xp.Dot(patterns, v_patterns) >= demand) p.controls.outputlog = 0 for npass in range(max_cols): print('=============Iteration {}==========='.format(npass)) p.optimize() cs_lp_sol = p.getSolution() duals = p.getDuals() obj_val = p.attributes.objval z, x_best = solve_knapsack_problem(n_widths, np.array(duals), np.array(widths), max_width, demand) if z < 1 + eps: print('No profitable column found') break print('New pattern found with marginal cost {}'.format(z - 1)); print(x_best) cur_pat = p.addVariable(name='pat_{}'.format(n_widths + npass), vartype=xp.continuous) for i in N: p.chgCoef(i, cur_pat, x_best[i]) p.chgObj([cur_pat], [1]) print("LP Solution after column generation. Obj: {}, Sol = {}".format(obj_val, cs_lp_sol)) n_final_patterns = p.getAttrib('cols') print('Number of patterns in the final LP = {}'.format(n_final_patterns)) # Change the variables to integers. p.chgColType([i for i in range(n_final_patterns)], ['I' for _ in range(n_final_patterns)]) p.optimize() integer_obj_val = p.attributes.objval int_sol = p.getSolution() print("MIP Solution with final columns. Obj: {}, Sol = {}".format(integer_obj_val, int_sol)) if __name__ == '__main__': cutting_stock_model()