Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 2022-04-29 16:16:18 +0200

Max Alekseyev gravatar image

To find just one solution you can use pose the problem as MILP and specify the objective as None. See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To find just one solution you can use pose the problem as MILP and specify the objective as None. See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To find just one solution you can pose the problem as MILP and specify the objective as with None.set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides in integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides in into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='GLPK' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print(mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron().volume()

Since the volume is zero, it means is has no internal points and thus the problem has no solutions.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='GLPK' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print(mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron().volume()

Since the volume is zero, it means is the polyhedron defined by $Ax\geq 0$ has no internal points and thus the problem inquality $Ax>0$ has no solutions.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='GLPK' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print(mylp.get_values(x))
print('Solution:', mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron().volume()

Since the volume is zero, the polyhedron defined by $Ax\geq 0$ has no internal points and thus the inquality $Ax>0$ has no solutions.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

from sage.numerical.mip import MIPSolverException

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='GLPK' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print('Solution:', mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron().volume()

Since the volume is zero, the polyhedron defined by $Ax\geq 0$ has no internal points and thus the inquality $Ax>0$ has no solutions.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

from sage.numerical.mip import MIPSolverException

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='GLPK' solver='ppl' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print('Solution:', mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron().volume()
mylp.polyhedron().volume(base_ring=QQ)

Since the volume is zero, the polyhedron defined by $Ax\geq 0$ has no internal points and thus the inquality $Ax>0$ has no solutions.

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

from sage.numerical.mip import MIPSolverException

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='ppl' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print('Solution:', mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron().volume(base_ring=QQ)
mylp.polyhedron(base_ring=QQ).volume()

Since the volume is zero, the polyhedron defined by $Ax\geq 0$ has no internal points and thus the inquality $Ax>0$ has no solutions.

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