coin-or / Cbc

COIN-OR Branch-and-Cut solver
Other
815 stars 115 forks source link

Clipping Function #412

Open FlorianStraubIng opened 3 years ago

FlorianStraubIng commented 3 years ago

Dear CBC-team,

thank you for providing CBC. I'm currently running CBC using the Python-MIP interface.

I’m trying to integrate a clipping function into a mip. I do this by formulating constraints that clip input values to the interval [clip_lower, clip_upper] using the BigM notation. I suspect that CBC may not solve the problem correctly for all input values. I've already described the issue in the Python-MIP Github (https://github.com/coin-or/python-mip/issues/194), but I think that no one knows why this occurs.

"""
Demonstration of mip implementation of a clipping function. Constrains clip input values to the interval
[clip_lower, clip_upper]
"""
from typing import List
import mip
import numpy as np
import matplotlib.pyplot as plt

def solve_mip(x, upper_clip=1, lower_clip=0):
    """
    Solves a mip that clips the input value between 0 and the defined upper border.

    Args:
        x (float): input value
        upper_clip (float): Max value where the input value is to be clipped
        lower_clip (float): Min value where the input value is to be clipped

    Returns:
        upper_bound and lower_bound of the input_value. (Note that we want this mip to have a defined output, so
        lower_bound and upper_bound are expected to coincide.)
    """
    m = mip.Model(solver_name='CBC', sense=mip.MAXIMIZE)
    y = m.add_var()
    d = m.add_var(var_type=mip.BINARY)
    d2 = m.add_var(var_type=mip.BINARY)
    M = 1000

    # d == 1 -> y = upper_clip; d == 0 & d2 == 1 -> y = x
    m.add_constr(y <= upper_clip)
    m.add_constr(y <= x + (1 - d2) * M)  # d2 == 1 & d == 0 -> y = x
    m.add_constr(y >= x - M * d)
    m.add_constr(y >= upper_clip - M * (1 - d))

    # lower boundary
    m.add_constr(y >= lower_clip)
    m.add_constr(y <= lower_clip + d2 * M)  # d2 == 0 -> y = lower_clip

    # To make sure that the solution that is found is the only solution, the mip is solved minimizing and maximizing the
    # output values -> if lower == upper we can guarantee that the mip's solution is well defined.
    m.objective = y
    m.optimize()
    upper = y.x

    m.objective = -y
    m.optimize()
    lower = y.x
    return lower, upper

def plot_values(plotdata: List[tuple]):
    """
    Plots multiple datasets in the same figure.

    Args:
        plotdata (List[tuple]): Each tuple contains (x_data, y_data, color_definition). color_definition can be 'r',
            'b', ...
    """
    fig = plt.figure()

    ax = fig.add_subplot(111)

    for x, y, color in plotdata:
        ax.plot(x, y, color)
        ax.spines['left'].set_position('zero')
        ax.spines['right'].set_color('none')
        ax.spines['bottom'].set_position('zero')
        ax.spines['top'].set_color('none')

        # remove the ticks from the top and right edges
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')
    plt.show()

if __name__ == '__main__':
    x_in = np.linspace(-2, 2, 500)
    y_lower = []
    y_upper = []

    # run the mip for all values on x-axis
    for x in x_in:
        lower, upper = solve_mip(x)
        y_lower.append(lower)
        y_upper.append(upper)

    # If the constraints are strict enough, we expect the data from the lower and upper boundary to be identical
    plot_data = [
        (x_in, y_lower, 'blue'),
        (x_in, y_upper, 'red')
    ]

    print('lower and upper boundaries are identical: ', y_lower == y_upper)
    plot_values(plot_data)

clip_values

As you can see in the plot, this mip can be solved for most of the input values. However there are some values that return None/Infeasible. Since these values are only a small subset of points in the linear domain of the problem, I assume that the mip is actually feasible in these regions. Might this be a solver related issue? Or a numerical issue? Are there any CBC properties that I could change so this problem would be feasible for all input values? Or do you think that this kind of problem can't be solved reliably?

Any hints on what might cause this issue would be much appreciated.

Thank you in advance!

Best,

Florian

tosttost commented 3 years ago

I'm not sure what you are trying to archive but it looks like a 'bad' idea (big-M with generic M is nearly always a bad idea). I suspect numerical issues of some kind. If cbc can be made more stable for this kind of problem is another question.

Generally the chances that someone solves your issue are higher if it is easy to reproduce (most if not all contributors work on CBC in their spare time). In this case:

I suggest that you identify one of the problematic values and write the MIP to a .lp / .mps file using https://docs.python-mip.com/en/latest/classes.html#mip.Model.write Then write down a feasible / expected solution and make sure it is in fact feasible. This way you can rule out modelling issues. If feasible post the file and the expected solution here.

FlorianStraubIng commented 3 years ago

Hi tostost,

thank you for your reply and the hints!

I don't expect anyone to debug the script, but I didn't want to provide the plots without a source;)

Regarding the CBC version I got indeed different results when I used the latest CBC release. But I still got infeasible results where I expected a solution to be found. I extracted some example lp files that unexpectedly returned infeasible: lpfiles_infeasible.zip, whereas this returns feasible lp_file_feasible.zip

In all of the lp files above the correct variable assignments would be: var(0)=x var(1)=0 var(2)=1

I used the default CBC settings for the solving process.

As you suspected I discovered that providing different BigM values indeed changed the outcome of the optimization, but didn't lead to consistent solver behavior. However, I found that the input x-values that I provided seemed to be crucial to the outcome.

The constraints that represent the clipping function are supposed to implement

Please note that I provided input values in the interval [-2,2] with 500 samples. So the step size is 4/499. (4 is the range of the interval, the first value is the lower bound (-2) and there are 499 steps left to reach the upper bound)

Those input values seem to cause the problems

M=1000, 500 equally spaced values from [-2,2] -> Stepsize is 4/499: grafik

Only for very small M the problem can be solved with the 500 equally spaced values. M=3, 500 equally spaced values from [-2,2] -> Stepsize is 4/499: grafik (no infeasible values)

For M=4 CBC returns infeasible for many input values. M=4, 500 equally spaced values from [-2,2] -> Stepsize is 4/499 grafik

M=1000, 501 equally spaced values from [-2,2] -> Stepsize is 4*0.002 grafik

M=1000, 513 equally spaced values from [-2,2] -> Stepsize is 4/(2^9) grafik (no infeasible values)

M=1000, 64 equally spaced values from [-2,2] -> Stepsize is 4/63 grafik

M=1000, 65 equally spaced values from [-2,2] -> Stepsize is 4/(2^6) grafik (no infeasible values)

The experiments especially interesting to me are those where I sampled 65 and 513 values from the interval. Since the first x-value is the lower bound we have left 64 and 512 x-values to reach the upper bound. Since those values – 64 and 512 – are our denominators when calculating x and are powers of 2, they lead to x-values that can be represented exactly in the binary system. All x-values that cause infeasibility are within the interval (0,1). So, it seems like the limited accuracy of floating-point numbers lead to those numerical instabilities.. Is this expected? Can this be solved for cases where I can't provide input values that are easy to represent as floats in binary? Any hints would be appreciated :)

Tank you in advance,

Florian

tosttost commented 3 years ago

That looks like a problem with preprocessing:

$ cbc clipping_demo_m1000_x28e-2.lp
Welcome to the CBC MILP Solver 
Version: 2.10.5 
Build Date: Apr  8 2020 

command line - cbc clipping_demo_m1000_x28e-2.lp -solve (default strategy 1)
Continuous objective value is 0 - 0.00 seconds
Cgl0000I Cut generators found to be infeasible! (or unbounded)
Pre-processing says infeasible or unbounded
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

If I turn preprocessing off (-preprocess off for command line Cbc), all 3 files are solved correctly. I don't use python-MIP, so I don't know if/how to turn off preprocessing.

I don't know enough about the internals of preprocessing in order to improve it - maybe @jjhforrest can help.

jjhforrest commented 3 years ago

Will try and have it fixed by Monday.

John Forrest

On 26/06/2021 15:51, tosttost wrote:

That looks like a problem with preprocessing:

|$ cbc clipping_demo_m1000_x28e-2.lp -preprocess off -solve Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Apr 8 2020 command line

  • cbc clipping_demo_m1000_x28e-2.lp -solve (default strategy 1) Continuous objective value is 0 - 0.00 seconds Cgl0000I Cut generators found to be infeasible! (or unbounded) Pre-processing says infeasible or unbounded Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00 |

If I turn preprocessing off (-preprocess off for command line Cbc), all 3 files are solved correctly. I don't use python-MIP, so I don't know if/how to turn off preprocessing.

I don't know enough about the internals of preprocessing in order to improve it - maybe @jjhforrest https://github.com/jjhforrest can help.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/coin-or/Cbc/issues/412#issuecomment-869013940, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABWJYHGBNROUE7R5C4GA7TLTUXSP5ANCNFSM46ZAPPEQ.

jjhforrest commented 3 years ago

Should be fixed in master.

On 26/06/2021 15:51, tosttost wrote:

That looks like a problem with preprocessing:

|$ cbc clipping_demo_m1000_x28e-2.lp -preprocess off -solve Welcome to the CBC MILP Solver Version: 2.10.5 Build Date: Apr 8 2020 command line

  • cbc clipping_demo_m1000_x28e-2.lp -solve (default strategy 1) Continuous objective value is 0 - 0.00 seconds Cgl0000I Cut generators found to be infeasible! (or unbounded) Pre-processing says infeasible or unbounded Total time (CPU seconds): 0.00 (Wallclock seconds): 0.00 |

If I turn preprocessing off (-preprocess off for command line Cbc), all 3 files are solved correctly. I don't use python-MIP, so I don't know if/how to turn off preprocessing.

I don't know enough about the internals of preprocessing in order to improve it - maybe @jjhforrest https://github.com/jjhforrest can help.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/coin-or/Cbc/issues/412#issuecomment-869013940, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABWJYHGBNROUE7R5C4GA7TLTUXSP5ANCNFSM46ZAPPEQ.

FlorianStraubIng commented 3 years ago

Wow that was quick!

Thank you very much guys, much appreciated!