coin-or / python-mip

Python-MIP: collection of Python tools for the modeling and solution of Mixed-Integer Linear programs
Eclipse Public License 2.0
518 stars 92 forks source link

Using python-mip library with cvxpy syntax #182

Closed Joanna-Wojdylo closed 3 years ago

Joanna-Wojdylo commented 3 years ago

I need to use CBC solver for mixed integer optimization problem, however in the target environment I cannot use CBC solver installed as an outsource software, it has to be a part of the python library. To overcome this problem I found mip library https://pypi.org/project/mip/ https://docs.python-mip.com/en/latest/index.html which comes with build-in CBC solver, which can be used only having this library imported, with no need for CBC solver installed separately. My problem is, I already have tons of code written with cvxpy (using this separate CBC solver). Now the question is, is there any possibility to use CBC built in mip library, but use it from regular cvxpy interface? Without changing the code, rewriting everything to mip sytax etc.

Sample code I would need to rewrite to mip syntax:

import numpy as np
import cvxpy as cp
import cvxopt 
import mip

def run_sample_optimization():
    demand = np.array([[100, 500, 30], [20, 200, 50], [150, 15, 35], [10, 5, 25]])
    product_supply = np.array([550, 200, 170, 40])

    allocation = cp.Variable(demand.shape, integer=True)

    objective = cp.Maximize(cp.sum(allocation/demand))

    constraints =[cp.sum(allocation, axis=1) <= product_supply,
                allocation <= demand,
                allocation >= 0]                        

    problem = cp.Problem(objective, constraints)

    optimal_value = problem.solve(solver=cp.GLPK_MI) # <-- it would be perfect to link somehow from this place to CBC implemented in mip library

    print('product supply:', product_supply)
    print('demand:\n', demand)
    print('allocation:\n', allocation.value)
    print('calculated score:', optimal_value)
    return product_supply, demand, allocation.value, optimal_value

Many thanks in advance!

tkralphs commented 3 years ago

You can call Cbc directly from cvxpy, see the list of supported solvers here: https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver.

jurasofish commented 3 years ago

fun problem. cvxpy is fortunately really modular, so you're in luck and you can do exactly what you're asking for easily.

I've modified cvxpy's cylp (thanks @tkralphs) interface to call python-mip. If you diff my code with the cvxpy code linked you'll see they're very similar.

Once you have this custom solver you can use it in cvxpy as described here and here

I made a repo for this https://github.com/jurasofish/mip_cvxpy and might flesh it out a bit more sometime.

Keep in mind the code I've written here is by no means robust, tested, or well thought out. If you do break it please let me know.

# main.py
import numpy as np
import cvxpy as cp
from mip_conif import PYTHON_MIP

def run_sample_optimization(solver):
    demand = np.array([[100, 500, 30], [20, 200, 50], [150, 15, 35], [10, 5, 25]])
    product_supply = np.array([550, 200, 170, 40])
    allocation = cp.Variable(demand.shape, integer=True)
    a = cp.Variable(shape=1)  # To help with debugging
    objective = cp.Maximize(cp.sum(allocation/demand) + cp.sum(a))
    constraints = [
        cp.sum(allocation, axis=1) <= product_supply,
        allocation <= demand,
        allocation >= 0,
        cp.sum(a) == 0,
    ]
    problem = cp.Problem(objective, constraints)

    solver_name = solver if isinstance(solver, str) else solver.name()
    print('solving with', solver_name)

    optimal_value = problem.solve(solver=solver)

    print('product supply:', product_supply)
    print('demand:\n', demand)
    print('allocation:\n', allocation.value)
    print('calculated score:', optimal_value)
    return product_supply, demand, allocation.value, optimal_value

def main():
    for solver in (cp.GLPK_MI, PYTHON_MIP()):
        run_sample_optimization(solver)

if __name__ == "__main__":
    main()
# mip_conif.py
"""

Based on
    https://github.com/cvxgrp/cvxpy/
    blob/5d5c7d606e39b3ea4b54391f772c7e3dc38ede20/cvxpy/reductions/
    solvers/conic_solvers/cbc_conif.py

Copyright 2016 Sascha-Dominic Schnug
Copyright 2021 Michael Jurasovic

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""

import cvxpy.settings as s
from cvxpy.reductions.solvers.conic_solvers.cbc_conif import (
    CBC,
    dims_to_solver_dict,
)
from cvxpy.reductions.solvers.conic_solvers.conic_solver import ConicSolver
from cvxpy.reductions.solution import Solution, failure_solution
import numpy as np

class PYTHON_MIP(CBC):  # uppercase consistent with cvxopt
    """ An interface to the python-mip solver
    """

    # Solver capabilities.
    MIP_CAPABLE = True
    SUPPORTED_CONSTRAINTS = ConicSolver.SUPPORTED_CONSTRAINTS

    def name(self):
        """The name of the solver.
        """
        return "PYTHON_MIP"

    def import_solver(self):
        """Imports the solver.
        """
        import mip
        _ = mip  # For flake8

    def accepts(self, problem):
        """Can python-mip solve the problem?
        """
        # TODO check if is matrix stuffed.
        if not problem.objective.args[0].is_affine():
            return False
        for constr in problem.constraints:
            if type(constr) not in self.SUPPORTED_CONSTRAINTS:
                return False
            for arg in constr.args:
                if not arg.is_affine():
                    return False
        return True

    def apply(self, problem):
        """Returns a new problem and data for inverting the new solution.

        Returns
        -------
        tuple
            (dict of arguments needed for the solver, inverse data)
        """
        data, inv_data = super(PYTHON_MIP, self).apply(problem)
        variables = problem.x
        data[s.BOOL_IDX] = [int(t[0]) for t in variables.boolean_idx]
        data[s.INT_IDX] = [int(t[0]) for t in variables.integer_idx]

        return data, inv_data

    def invert(self, solution, inverse_data):
        """Returns the solution to the original problem given the inverse_data.
        """
        status = solution['status']

        if status in s.SOLUTION_PRESENT:
            opt_val = solution['value'] + inverse_data[s.OFFSET]
            primal_vars = {inverse_data[self.VAR_ID]: solution['primal']}
            return Solution(status, opt_val, primal_vars, None, {})
        else:
            return failure_solution(status)

    def solve_via_data(self, data, warm_start, verbose, solver_opts, solver_cache=None):
        # Import basic modelling tools of cylp
        import mip

        c = data[s.C]
        b = data[s.B]
        A = data[s.A]
        dims = dims_to_solver_dict(data[s.DIMS])

        n = c.shape[0]

        # Problem
        model = mip.Model()

        # Variables
        x = []
        bool_idxs = set(data[s.BOOL_IDX])
        int_idxs = set(data[s.INT_IDX])
        for i in range(n):
            if i in bool_idxs:
                x.append(model.add_var(var_type=mip.BINARY))
            elif i in int_idxs:
                x.append(model.add_var(var_type=mip.INTEGER))
            else:
                x.append(model.add_var())

        # Constraints
        # eq
        def add_eq_constraints(_model):
            coeffs = A[0:dims[s.EQ_DIM], :]
            vals = b[0:dims[s.EQ_DIM]]
            for i in range(coeffs.shape[0]):
                coeff_list = np.squeeze(np.array(coeffs[i].todense())).tolist()
                expr = mip.LinExpr(variables=x, coeffs=coeff_list)
                _model += expr == vals[i]
        add_eq_constraints(model)

        # leq
        def add_leq_constraints(_model):
            leq_start = dims[s.EQ_DIM]
            leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM]
            coeffs = A[leq_start:leq_end, :]
            vals = b[leq_start:leq_end]
            for i in range(coeffs.shape[0]):
                coeff_list = np.squeeze(np.array(coeffs[i].todense())).tolist()
                expr = mip.LinExpr(variables=x, coeffs=coeff_list)
                _model += expr <= vals[i]
        add_leq_constraints(model)

        # Objective
        model.objective = mip.minimize(mip.LinExpr(variables=x, coeffs=c.tolist()))

        model.verbose = verbose
        status = model.optimize()

        status_map = {
            mip.OptimizationStatus.OPTIMAL: s.OPTIMAL,
            mip.OptimizationStatus.INFEASIBLE: s.INFEASIBLE,
            mip.OptimizationStatus.INT_INFEASIBLE: s.INFEASIBLE,
            mip.OptimizationStatus.NO_SOLUTION_FOUND: s.INFEASIBLE,
            mip.OptimizationStatus.ERROR: s.SOLVER_ERROR,
            mip.OptimizationStatus.UNBOUNDED: s.UNBOUNDED,
            mip.OptimizationStatus.CUTOFF: s.INFEASIBLE,
            mip.OptimizationStatus.FEASIBLE: s.OPTIMAL_INACCURATE,
            mip.OptimizationStatus.LOADED: s.SOLVER_ERROR,  # No match really
        }

        solution = {
            "status": status_map[status],
            "primal": [var.x for var in x],
            "value": model.objective_value,
        }

        return solution
Joanna-Wojdylo commented 3 years ago

Thank you very much @jurasofish! Your answer was really helpful. I have one more question, because on the PYTHON-MIP website they are saying CBC is the default solver used in mip Model. However, I need to be 100% sure that it really is CBC and ideally force mip Model to use solver_name=CBC. I tried to simply add it in mip_conif.py while calling mip.Model():

    def solve_via_data(self, data, warm_start, verbose, solver_opts, solver_cache=None):
        # Import basic modelling tools of cylp
        import mip

        c = data[s.C]
        b = data[s.B]
        A = data[s.A]
        dims = dims_to_solver_dict(data[s.DIMS])

        n = c.shape[0]

        # Problem
        model = mip.Model(solver_name=CBC)# <--- right here

but this is probably not the brightest idea. It leads to further errors:

 File "C:\ \Anaconda3\envs\cbc_testing_env\lib\site-packages\cvxpy\problems\problem.py", line 436, in solve
    return solve_func(self, *args, **kwargs)
  File "C:\ \Anaconda3\envs\cbc_testing_env\lib\site-packages\cvxpy\problems\problem.py", line 927, in _solve
    solution = solving_chain.solve_via_data(
  File "C:\ \Anaconda3\envs\cbc_testing_env\lib\site-packages\cvxpy\reductions\solvers\solving_chain.py", line 340, in solve_via_data
    return self.solver.solve_via_data(data, warm_start, verbose,
  File "C:\ \Anaconda3\envs\cbc_testing_env\lib\site-packages\cvxpy\reductions\solvers\conic_solvers\mip_conif.py", line 107, in solve_via_data
    model = mip.Model(solver_name=CBC)
  File "C:\ \Anaconda3\envs\cbc_testing_env\lib\site-packages\mip\model.py", line 82, in __init__
    if self.solver_name.upper() in ["GUROBI", "GRB"]:

I saw, that it is possible to run PYTHON_MIP("CBC") but again, I am not 100% sure if this really calls CBC solver. Would you have any further recommendations for me, what can be done to force mip to use CBC solver and be 100% sure?

jurasofish commented 3 years ago

almost, should work if you do model = mip.Model(solver_name=mip.CBC)

tkralphs commented 3 years ago

If I'm not mistaken, the code that is failing in the comment above is here:

https://github.com/coin-or/python-mip/blob/161bd783994ca157fa7cc4b03a18211fdcd27bd4/mip/model.py#L82-L101

It looks like it is expecting a string, so you would quote CBC.

model = mip.Model(solver_name="CBC")
Joanna-Wojdylo commented 3 years ago

Yes, both ways suggested by you guys work, solver_name=mip.CBC and solver_name = "CBC", so big thanks to you guys! However, as you saw on https://github.com/coin-or/Cbc/issues/371 such prepared mip-python interface works perfectly on a simple optimization problem, but freezes on the real optimization problem I need to solve (and which can be easily solved with external CBC in few minutes).

Joanna-Wojdylo commented 3 years ago

I think this issue can be closed now, the solution and an additional discussion can be found here https://github.com/cvxgrp/cvxpy/issues/1265.

Thank you very much guys!