column_generation.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import sys
import types
from amplpy import AMPL, Parameter

# Example description
# This example uses Gilmore-Gomory column generation procedure for 
# the cutting-stock (roll trim) problem.
# Keeps the knapsack problem as an AMPL model, and adds columns
# to the cutting-stock via ampls; finally it imports the results in AMPL.

import amplpy_copt as ampls
SOLVER = "copt"

WIDTHS = [20, 45, 50, 55, 75]
ORDERS = [48, 35, 24, 10, 8]
RAW_WIDTH = 110

# We represent a pattern as a dictionary width -> numberofcuts
# where we only store the nonzeroes
def generate_default_patterns(widths : list):
    return [{w : RAW_WIDTH // w} for w in widths]
 
def add_patterns(ampl: AMPL, patterns : list):
    
    npat = int(ampl.param["nPatterns"].value())
    ampl.param["nPatterns"].set(len(patterns))
    ampl.eval("update data rolls;")
    ampl.param['rolls'] = {
        (w, 1+i): n
        for i in range(npat,len(patterns)) 
        for w, n in patterns[i].items()
    }


def cutting_stock_model():
    """Create an instance of AMPL and a model"""
    a = AMPL()
    a.option["presolve"]=0
    a.option["solver"]=SOLVER

    a.eval("""param nPatterns integer >= 0;
    set PATTERNS = 1..nPatterns;  
    set WIDTHS;                 
    param order{ WIDTHS } >= 0;
    param rawWidth;
    param rolls{ WIDTHS,PATTERNS } >= 0, default 0;

    var Cut{ PATTERNS } integer >= 0;

    minimize TotalRawRolls : sum{ p in PATTERNS } Cut[p] + to_come;
    subject to OrderLimits{ w in WIDTHS }:
    sum{ p in PATTERNS } rolls[w, p] * Cut[p]+ to_come >= order[w];""") 

    a.param["nPatterns"]=0
    a.set["WIDTHS"]=WIDTHS
    a.param["order"]=ORDERS
    a.param["rawWidth"]=RAW_WIDTH

    return a

def generate_pattern(ampl: AMPL, duals : list):
    ampl.param["price"]=duals
    ampl.solve()
    ampl.eval("display price;")
    reduced_cost=ampl.obj['Reduced_Cost'].value()
    if reduced_cost < -0.00001:
        return ampl.var["Use"].getValues().to_dict()
    return None

def knapsack_model():
    a = AMPL()
    a.option["presolve"]=0
    a.option["solver"]=SOLVER
    a.eval("""set WIDTHS; 
    param rawWidth;
    param price {WIDTHS} default 0.0;
    var Use {WIDTHS} integer >= 0;
    minimize Reduced_Cost: 1 - sum{ i in WIDTHS } price[i] * Use[i];
    subject to Width_Limit: sum{ i in WIDTHS } i * Use[i] <=rawWidth;""")
    a.set["WIDTHS"]=WIDTHS
    a.param["rawWidth"]=RAW_WIDTH
    return a

USE_AMPLS_ENTITY_RECORDING = True
def run_example():
    # Decalare the two models in AMPL
    cs = cutting_stock_model()
    knap = knapsack_model()

    # Generate starting patterns (each patterns simply cuts the roll
    # all at the same width)
    patterns = generate_default_patterns(WIDTHS)
    # Add them to the AMPL instance
    add_patterns(cs, patterns)
    
    # Export the (relaxed) cutting stock model to ampls
    cs.option["relax_integrality"]=1
    options = ["outlev=1", "pre:maxrounds=0"] if SOLVER=="scip" else None
    ampls_cs = cs.to_ampls(SOLVER, options)

    

    while True: # Column generation happens in the solver
        # Optimize the cutting stock model, get the dual vector,
        # use it in the knapsack model (in AMPL) to generate 
        # a new pattern 
        ampls_cs.optimize()
        duals = ampls_cs.get_dual_vector()
        print(duals)
        p=generate_pattern(knap, duals)
        
        if p is None: break # No new pattern has been found, finish

        # Keep track of the patterns we add
        patterns.append(p)
        
        # Create the variable indices and coefficients
        index=0
        indices=[]
        coeffs=[]
        for _, c in p.items():
            if c != 0:
                indices.append(index)
                coeffs.append(c)
            index+=1
        # Add variable in the cutting stock model
        # Note the last parameter: we add the variable as "relaxed", so that it is 
        v=ampls_cs.addVariable(indices, coeffs, 0, 10000000,
                                1, ampls.VarType.Integer, True)
        
        if USE_AMPLS_ENTITY_RECORDING:
            ampls_cs.record(v)
    
    # Add all the patterns that has been generated in the loop above to
    # the AMPL version of the cutting stock model, then solve the 
    # integer problem to get the final result

    if USE_AMPLS_ENTITY_RECORDING:
        cs.import_ampls_solution(ampls_cs, import_entities=True)
    else:
        add_patterns(cs, patterns)
        cs.import_ampls_solution(ampls_cs)
    cs.option["relax_integrality"]=0
    cs.solve()
    cs.display("TotalRawRolls, rolls;")
    cs.display("{i in 1.._nvars} _varname[i]");
    assert cs.get_current_objective().value() == 47
    
run_example()