Decision Optimization

 View Only
  • 1.  Max weighted coverage problem: CP versus MILP formulation, CP Optimizer much slower?

    Posted Sun May 02, 2021 11:36 AM

    How do you efficiently express a max weighted set coverage problem in CP Optimizer? The conventional MILP formulation is solved very rapidly by CPLEX (via docplex.mp), but the same formulation in CP Optimizer (via docplex.cp) takes seemingly forever to solve.

    The reason that I am trying to do a max weighted set coverage problem in CP Optimizer is that this is a toy version of my real problem contains scheduling with optional intervals, with the presence of each interval being identified with the presence of a corresponding set.

    Imports, problem setup

    from docplex.cp import model as cp
    from docplex.mp import model as mp
    import numpy as np
    
    n_elements = 10000
    n_sets = 1000
    n_elements_per_set = 100
    max_sets_present = 10
    
    # Random weights for each element
    weights = np.random.uniform(size=n_elements)
    weights /= weights.sum()
    
    # Random assignments of elements to sets
    assigns = np.asarray([
        np.random.choice(n_elements, n_elements_per_set, replace=False)
        for _ in range(n_sets)])
    
    assigns_inverse = [np.nonzero(assigns == i)[0] for i in range(n_elements)]

    MILP version

    with mp.Model() as m:
        elems_present = m.binary_var_list(n_elements)
        sets_present = m.binary_var_list(n_sets)
        m.add_(m.sum(sets_present[j] for j in assign_inverse)
               >= elem_present for assign_inverse, elem_present
               in zip(assigns_inverse, elems_present))
        m.add_(m.sum(sets_present) <= max_sets_present)
        m.maximize(m.scal_prod(elems_present, weights))
    
        m.print_information()
        %time solution = m.solve()
    
    objective_value = solution.get_objective_value()
    print('objective_value =', objective_value)

    Output:

    Model: docplex_model1
    - number of variables: 11000
    - binary=11000, integer=0, continuous=0
    - number of constraints: 10001
    - linear=10001
    - parameters: defaults
    - objective: maximize
    - problem type is: MILP
    CPU times: user 1min 23s, sys: 1.69 s, total: 1min 25s
    Wall time: 25.6 s
    objective_value = 0.11137026929164126

    CP version

    with cp.CpoModel() as m:
        elems_present = m.binary_var_list(n_elements)
        sets_present = m.binary_var_list(n_sets)
        m.add(m.sum(sets_present[j] for j in assign_inverse)
              >= elem_present for assign_inverse, elem_present
              in zip(assigns_inverse, elems_present))
        m.add(m.sum(sets_present) <= max_sets_present)
        m.maximize(m.scal_prod(elems_present, weights))
    
        m.write_information()
        %time solution = m.solve(LogVerbosity='Verbose', Workers='Auto')
    
    objective_value, = solution.get_objective_values()
    print('objective_value =', objective_value)

    Output:

    Model: <ipython-input-4-9b2e732a693a>
    - source file: <ipython-input-4-9b2e732a693a>
    - modeling time: 0.41 sec
    - number of integer variables: 11000
    - number of interval variables: 0
    - number of expressions: 10003
    - number of expression nodes: 41008
    - operations: greaterOrEqual: 10000, lessOrEqual: 1, maximize: 1, scalProd: 1, sum: 10001
    (...never terminates)


    ------------------------------
    Leo Singer
    ------------------------------

    #DecisionOptimization


  • 2.  RE: Max weighted coverage problem: CP versus MILP formulation, CP Optimizer much slower?

    Posted Mon May 03, 2021 07:28 AM
    Edited by System Fri January 20, 2023 04:15 PM
    Dear Leo,

    Although the CP version of your program doesn't return from solve(), our tests shows display intermediate solutions are displayed and the solver continues its work as it is seeking to prove optimality. So for comparing the CP and MP performances you could choose to set a time limit to CP.
    The current CP model gets stuck at around 0.04 after around 20s, But you can  improve its performance by specifying a search phase that will reduce the number of decision variables for the problem.
    Simply adding

       m.add(m.search_phase(sets_present))

    will improve CPO results (the objective raises up to around 0.097).

    Also, I note that the randon seed of the random number generator isn't initialized which means you solve a different problem each time. I've done my experiments by setting the random seed to 0.
    It would be interesting for us to know what version of  COS you are using. 
    I hope this helps.

    Edit: In addition it would be interesting to test the scheduling model since it'd rely on a temporal relaxation similar to what the cplex model does.

    ------------------------------
    Renaud Dumeur
    ------------------------------



  • 3.  RE: Max weighted coverage problem: CP versus MILP formulation, CP Optimizer much slower?

    Posted Mon May 03, 2021 08:37 AM
    As a complement to Renaud's answer, it could be interesting to know if there is some sort of correlation between the interval variables and their corresponding set. I suppose that in the real problem, the sets are not random sets but are somehow related with the interval variables or other aspects of the scheduling problem. Maybe there could be a compact and efficient way to constrain the elements of the sets and count the total number of covered elements by exploiting this relation ...



    ------------------------------
    Philippe Laborie
    ------------------------------



  • 4.  RE: Max weighted coverage problem: CP versus MILP formulation, CP Optimizer much slower?

    Posted Mon May 03, 2021 09:10 AM
    The problem is planning and scheduling observations with a large field-of-view space telescope to cover a weighted region of interest. The elements represent a grid of discrete points of interest on the sky with associated weights, and the sets represent the points contained within the field of view of the telescope at a given attitude of the spacecraft.

    See https://www.youtube.com/watch?v=68voaJnBjWk&t=770s at 12m12s.

    ------------------------------
    Leo Singer
    ------------------------------



  • 5.  RE: Max weighted coverage problem: CP versus MILP formulation, CP Optimizer much slower?

    Posted Mon May 03, 2021 01:05 PM
    Very interesting ...

    There is indeed probably not much to exploit about the elements of the sets (the pixels) and the interval variables (the fields or the observations). You could of course group together the pixels belonging to exactly the same sets (but maybe this does not occur very often).

    I tried a very naive formulation using interval variables only (one optional interval variable of duration 1, a single no-overlap constraint and an horizon equal to max_sets_present that will limit the number of present intervals). The scheduling part is thus not realistic at all (no time-windows, no repetition of observations), it is here only to post an upper bound max_sets_present on the number of sets to get an idea of the solution quality.

    x = [ interval_var(size=1, optional=True, end=[0,max_sets_present]) for i in range(n_sets) ]
    z = [ max([0] + [presence_of(x[i]) for i in assigns_inverse[l]]) for l in range(n_elements) ]

    model.add(no_overlap(x))
    model.maximize(scal_prod(z, weights))

    sol = model.solve(LogVerbosity='Verbose', LogPeriod=1000000)

    CPO quickly find solutions that seem decent (around 0.1) after a few seconds. But the bound (here the upper bound) is terrible. This is not unusual for scheduling problems like this one. It is probable that you should compute upper bounds by using other techniques (like typically a MIP relaxation of the problem). It makes me think of a slightly similar satellite observation scheduling problem described in the following paper (though in this paper the objective function is more about the regularity of the observations than about the coverage):

    P. Laborie, B. Messaoudi. New Results for the GEO-CAPE Observation Scheduling Problem. Proc. 27th International Conference on Automated Planning and Scheduling, 2017.

    Finding a good solution (with CPO) and providing an optimality bound (with MIP) were done using different formulations.

    ------------------------------
    Philippe Laborie
    ------------------------------