Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 2022-06-03 16:30:03 +0100

Max Alekseyev gravatar image

There are two essential problems to consider here: (i) finding ONE solution or establishing that there are NONE; (ii) finding ALL solutions. In my experience, problem (i) is best solved with (M)ILP solvers while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only if one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = milp = MixedIntegerLinearProgram()
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

Here one solution is found momentarily (at least with Gurobi MILP solver), but I was not patient enough to get all solutions.

There are two essential problems to consider here: (i) finding ONE solution or establishing that there are NONE; (ii) finding ALL solutions. In my experience, problem (i) is best solved with (M)ILP solvers while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only if one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = milp = MixedIntegerLinearProgram()
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

Here one solution is found momentarily (at least with Gurobi MILP solver), but I was not patient enough to get all solutions.

There are two essential problems to consider here: (i) finding ONE solution or establishing that there are NONE; (ii) finding ALL solutions. In my experience, problem (i) is best solved with (M)ILP solvers while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only if one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = MixedIntegerLinearProgram()
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

Here In the given example, one solution is found momentarily (at least with Gurobi MILP solver), but I was not patient enough to get all solutions.them all.

There are two essential problems to consider here:

  • (i) finding ONE solution or establishing that there are NONE;
  • (ii) finding ALL solutions.

In my experience, problem (i) is best solved with (M)ILP solvers solvers, while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only if when one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = MixedIntegerLinearProgram()
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

In the given example, one solution is found momentarily (at least with Gurobi MILP solver), but I was not patient enough to get them all.

PS. Notice that for a polyhedron is created from MILP, it may be tricky to map the coordinates of its integral points back to the variables - see https://ask.sagemath.org/question/55520/

There are two essential problems to consider here:

  • (i) finding ONE solution or establishing that there are NONE;
  • (ii) finding ALL solutions.

In my experience, problem (i) is best solved with (M)ILP solvers, while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only when one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = MixedIntegerLinearProgram()
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

In the given example, one solution is found momentarily (at least with Gurobi MILP solver), but I was not patient enough to get them all.

PS. Notice that for a polyhedron is created extracted from MILP, it may be tricky to map the coordinates of its integral points back to the variables - see https://ask.sagemath.org/question/55520/

There are two essential problems to consider here:

  • (i) finding ONE solution or establishing that there are NONE;
  • (ii) finding ALL solutions.

In my experience, problem (i) is best solved with (M)ILP solvers, while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only when one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = MixedIntegerLinearProgram()
MixedIntegerLinearProgram(solver='ppl')
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

In the given example, one solution is found momentarily (at least with Gurobi MILP solver), but I was not patient enough to get them all.

PS. Notice that for a polyhedron extracted from MILP, it may be tricky to map the coordinates of its integral points back to the variables - see https://ask.sagemath.org/question/55520/

There are two essential problems to consider here:

  • (i) finding ONE solution or establishing that there are NONE;
  • (ii) finding ALL solutions.

In my experience, problem (i) is best solved with (M)ILP solvers, while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only when one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = MixedIntegerLinearProgram(solver='ppl')
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed with finding all')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

In the given example, one solution is found momentarily (at least (tested with Gurobi MILP solver), ppl and Gurobi solvers), but I was not patient enough to get them all.

PS. Notice that for a polyhedron extracted from MILP, it may be tricky to map the coordinates of its integral points back to the variables - see https://ask.sagemath.org/question/55520/

There are two essential problems to consider here:

  • (i) finding ONE solution or establishing that there are NONE;
  • (ii) finding ALL solutions.

In my experience, problem (i) is best solved with (M)ILP solvers, while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only when one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
    n = len(A)      # number of equations
    assert n == len(V)
    m = len(A[0])   # number of variables

    milp = MixedIntegerLinearProgram(solver='ppl')
    milp.set_objective( None )
    x = milp.new_variable(integer=True, nonnegative=True)
    for i in range(n):
        milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

    try:
        milp.solve()
    except MIPSolverException:
        return []       # no solutions

    X = milp.get_values(x)
    print(f'one solution is found: {list(milp.get_values(x).values())}\nproceed {[X[i] for i in range(m)]}')
    print('proceed with finding all')
all solutions')

    # --------------------------------------------------------------

    # we can extract polyhedron from MILP, or construct it explicitly
    # P = milp.polyhedron( backend='normaliz', base_ring=QQ ) 

    # explicit polyhedron construction:
    equalities = [ [-V[i]] + A[i] for i in range(n) ]
    inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
    P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
    print('Polyhedron:', P)

    return P.integral_points()

In the given example, one solution is found momentarily (tested with ppl and Gurobi solvers), but I was not patient enough to get them all.

PS. Notice that for a polyhedron extracted from MILP, it may be tricky to map the coordinates of its integral points back to the variables - see https://ask.sagemath.org/question/55520/

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59