Performance issue in CP-SAT

Question

Dears, I'm having performance issues on the code below. I know it is challenging to forecast 435 variables, but if i find the way to run the model below for a big amount of variables it will be great for me, otherwise if you can suggest alternative way i would be grateful. Here the code:

from __future__ import division
from __future__ import print_function

from ortools.sat.python import cp_model
from ortools.sat.python.cp_model import IntVar
import time
import datetime


def cp_sat_program():

    # We have 2 products Gold, Silver each of them has a profit and a bad debt value.
    # We have a dataset with a list of customers each one with a PD value.
    # We want to maximize the profit of the portfolio of customers keeping the expected loss under a threshold
    # and the number of reject less than 55%

    time_start = time.time()
    date_start = datetime.datetime.fromtimestamp(time_start).strftime('%Y-%m-%d %H:%M:%S')

    # Creates the constant variables.
    scaling = 1000  # i'm scaling since the algorithm works with integer only
    profit_gold = int(0.09 * scaling)
    profit_silver = int(0.023 * scaling)
    profit_reject = int(0 * scaling)
    customers_pd = [0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24,
                    25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94,
                    0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01,
                    25.26, 0.01, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94,
                    0.01, 0.24, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01,
                    25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94,
                    0.24, 0.01,0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24,
                    25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94,
                    0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24
                   ]

    num_customers = len(customers_pd)
    all_customers_ids = range(num_customers)

    customers_pd_scaled = {}
    # I'm scaling the PD of the customers because the system works with integers
    for i in all_customers_ids:
        customers_pd_scaled[i] = int(customers_pd[i] * scaling)
    # Bad Debt for each product:
    bad_debt_gold = int(profit_gold*0.63*scaling)
    bad_debt_silver = int(profit_silver*0.62*scaling)
    bad_debt_reject = int(0 * scaling)

    # Creates the model.
    model = cp_model.CpModel()
    # Creates the model variables
    suggested_products = {}
    for p in all_customers_ids:
        suggested_products[p] = [model.NewBoolVar('c' + str(p) + 'Gold'), model.NewBoolVar('c' + str(p) + 'Silver'),
                                 model.NewBoolVar('c' + str(p) + 'Reject')]

    scaled_r = model.NewIntVar(0, scaling * 1000, 'rejects')
    denom = model.NewIntVar(1, 3 * 1000, 'total Requests')
    division = model.NewIntVar(0, scaling * 1000, 'reject/total Requests')
    num_products = range(len(suggested_products[0]))

    # Creates the constraints.
    for i in all_customers_ids:
        model.Add(sum(suggested_products[i][b] for b in num_products) == 1)
    num_rejects = sum(suggested_products[i][2] for i in all_customers_ids)
    model.Add(scaled_r == (num_rejects * scaling))  # i'm scaling since the algorithm works with integer only
    model.Add(denom == num_customers)
    model.AddDivisionEquality(division, scaled_r, denom)  # calculate the division
    model.Add(division <= int(0.55 * scaling))  # percentage of rejects less than 55%

    # expected loss less or equal than 250.000:
    model.Add(sum(suggested_products[i][0] * customers_pd_scaled[i] * bad_debt_gold for i in all_customers_ids) +
              sum(suggested_products[i][1] * customers_pd_scaled[i] * bad_debt_silver for i in all_customers_ids) +
              sum(suggested_products[i][2] * customers_pd_scaled[i] * bad_debt_reject for i in all_customers_ids)
              <= int(250000 * scaling))

    # goal
    model.Maximize(sum(suggested_products[i][0] * profit_gold for i in all_customers_ids) +
                   sum(suggested_products[i][1] * profit_silver for i in all_customers_ids) +
                   sum(suggested_products[i][2] * profit_reject for i in all_customers_ids)
                   )
    # Creates a solver and solves the model.
    solver = cp_model.CpSolver()
    solver.parameters.num_search_workers = 8
    status = solver.Solve(model)
    goal = solver.ObjectiveValue()
    if status != cp_model.INFEASIBLE:
        for i in all_customers_ids:
            print('CUSTOMER NUMBER ' + str(i) + ': Gold =' + str(solver.Value(suggested_products[i][0]))
                  + '  Silver =' + str(solver.Value(suggested_products[i][1]))
                  + '  Reject =' + str(solver.Value(suggested_products[i][2])))

        print('Goal = %i' % goal)
        print("status = %s" % solver.StatusName(status))
    else:
        print("status = %s" % solver.StatusName(status))


    time_end = time.time()
    date_end = datetime.datetime.fromtimestamp(time_end).strftime('%Y-%m-%d %H:%M:%S')
    print('Date Start: ' + date_start)
    print('Date End: ' + date_end)

cp_sat_program()

If the value of the limit for that constrain below is 330000 it is very fast if the limit is 250000 it will never respond after one hour:

model.Add(sum(suggested_products[i][0] * customers_pd_scaled[i] * bad_debt_gold for i in all_customers_ids) +
          sum(suggested_products[i][1] * customers_pd_scaled[i] * bad_debt_silver for i in all_customers_ids) +
          sum(suggested_products[i][2] * customers_pd_scaled[i] * bad_debt_reject for i in all_customers_ids)
          <= int(330000 * scaling))

Answer 1

If you really want to solve with a solver, here is a cleaned up version

from __future__ import division
from __future__ import print_function

from ortools.sat.python import cp_model
from ortools.sat.python.cp_model import IntVar
import time
import datetime


def cp_sat_program():

    # We have 2 products Gold, Silver each of them has a profit and a bad debt value.
    # We have a dataset with a list of customers each one with a PD value.
    # We want to maximize the profit of the portfolio of customers keeping the expected loss under a threshold
    # and the number of reject less than 55%

    time_start = time.time()
    date_start = datetime.datetime.fromtimestamp(time_start).strftime('%Y-%m-%d %H:%M:%S')

    # Creates the constant variables.
    scaling = 1000  # i'm scaling since the algorithm works with integer only

    profit_gold = int(0.09 * scaling)
    profit_silver = int(0.023 * scaling)
    profit_reject = int(0 * scaling)

    customers_pd = [0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24,
                    25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94,
                    0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01,
                    25.26, 0.01, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94,
                    0.01, 0.24, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01,
                    25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94,
                    0.24, 0.01,0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24,
                    25.26, 5.94, 0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94,
                    0.24, 0.01, 25.26, 0.01, 25.26, 5.94, 0.01, 0.24, 0.01, 25.26, 5.94, 0.01, 25.26, 5.94, 0.24
                   ]

    num_customers = len(customers_pd)
    all_customers_ids = range(num_customers)

    customers_pd_scaled = [int(pd * scaling) for pd in customers_pd]

    # Bad Debt for each product:
    bad_debt_gold = int(profit_gold*0.63*scaling)
    bad_debt_silver = int(profit_silver*0.62*scaling)
    bad_debt_reject = int(0 * scaling)

    # Creates the model.
    model = cp_model.CpModel()

    # Creates the model variables
    suggested_products = {}
    for p in all_customers_ids:
        suggested_products[p] = (model.NewBoolVar('c' + str(p) + 'Gold'), model.NewBoolVar('c' + str(p) + 'Silver'),
                                 model.NewBoolVar('c' + str(p) + 'Reject'))

    num_products = range(len(suggested_products[0]))

    # Creates the constraints.
    for i in all_customers_ids:
        model.Add(sum(suggested_products[i][b] for b in num_products) == 1)
        model.Add(suggested_products[i][1] == 0)

    # At most 55% of rejects
    model.Add(sum(suggested_products[i][2] for i in all_customers_ids) <= int(num_customers * 0.55))

    # expected loss less or equal than 250.000:
    model.Add(sum(suggested_products[i][0] * customers_pd_scaled[i] * bad_debt_gold for i in all_customers_ids) +
              sum(suggested_products[i][1] * customers_pd_scaled[i] * bad_debt_silver for i in all_customers_ids)
              <= int(250000 * scaling))

    # goal
    model.Maximize(sum(suggested_products[i][0] * profit_gold for i in all_customers_ids) +
                   sum(suggested_products[i][1] * profit_silver for i in all_customers_ids))
    # Creates a solver and solves the model.
    solver = cp_model.CpSolver()
    solver.parameters.log_search_progress = True

    status = solver.Solve(model)
    goal = solver.ObjectiveValue()
    if status != cp_model.INFEASIBLE:
        for i in all_customers_ids:
            print('CUSTOMER NUMBER ' + str(i) + ': Gold =' + str(solver.Value(suggested_products[i][0]))
                  + '  Silver =' + str(solver.Value(suggested_products[i][1]))
                  + '  Reject =' + str(solver.Value(suggested_products[i][2])))

        print('Goal = %i' % goal)
        print("status = %s" % solver.StatusName(status))
    else:
        print("status = %s" % solver.StatusName(status))


    time_end = time.time()
    date_end = datetime.datetime.fromtimestamp(time_end).strftime('%Y-%m-%d %H:%M:%S')
    print('Date Start: ' + date_start)
    print('Date End: ' + date_end)

cp_sat_program()

It returns an assignment with value 6030 instantaneously.

Answer 2

I believe you can craft a near optimal solution.

First profit_gold >> profit_silver, while bad_debt_silver ~= bad_debt_gold. So you can assume that you will only assign gold.

Thus the problem now becomes assigning the maximum number of goals while staying under the 250000 limit. To do this, sort customers by customer_pd, then select then 1 by 1 until you hit the limit.