Ask Your Question
1

A linear program returns infeasible but I think to know a feasible solution.

asked 2021-12-10 09:52:22 +0100

CyrilleP gravatar image

updated 2023-05-19 22:07:02 +0100

FrédéricC gravatar image

I have the following linear programming program which is solved without problem by SageMath.

   nv=6 #nombre de contraintes
mv=4 #nombre de variables
F=10
a=1.5
Av= matrix(RR,nv,mv,[1,0,0,0,
                     0,1,0,0,
                     0,0,1,0,
                     0,0,0,1,
                     1,1,1,1,
                     1,-a,0,-a])     #les coefficients
bvmin=vector(RR,[3/20*F,3/20*F,11/40*F,3/20*F,17/20*F,0]) #les bornes inférieures pour les contraintes
bvmax=vector(RR,[9/20*F,9/20*F,9/20*F,9/20*F,17/20*F,9/20*F]) #les bornes supérieures pour les contraintes
show(LatexExpr('A='),Av)
show(LatexExpr('bmin='),bvmin)
show(LatexExpr('bmax='),bvmax)
#Création du programme
v=MixedIntegerLinearProgram(maximization=True, solver='GLPK') 
#On pose les nouvelles variables allant de x_0 à x_5
ind=[0..mv-1]
x=v.new_variable(integer=False,indices=ind) 
#On utilise la fonction linéaire pour les contraintes
Bv=Av*x
#On fixe l'objectif
v.set_objective(x[0])
for i in range(mv) :
    v.set_min(x[i],0)
    v.set_max(x[i],F)
for i in range(0,nv-1):
    v.add_constraint(Bv[i],min=bvmin[i], max=bvmax[i])    
v.show()
v.solve()
xx=v.get_values(x)
show(xx)

Now I want to change the objective to $maxmin(x_0,x_1,x_2,x_3)$. So I follow the usual approach which leads to the addition of a variable says $x_4$ and four new constraints $x_0\geq x_4$,$x_1\geq x_4$,$x_2\geq x_4$ and $x_3 \geq x_4$ as in the following code.

nv=10 #nombre de contraintes
mv=5 #nombre de variables
F=10
a=1.5
Av= matrix(RR,nv,mv,[1,0,0,0,0,
                     0,1,0,0,0,
                     0,0,1,0,0,
                     0,0,0,1,0,
                     1,1,1,1,0,
                     1,-a,0,-a,0,
                     1,0,0,0,-1,
                     0,1,0,0,-1,
                     0,0,1,0,-1,
                     0,0,0,1,-1])     #les coefficients
bvmin=vector(RR,[3/20*F,3/20*F,11/40*F,3/20*F,17/20*F,0,0,0,0,0]) #les bornes inférieures pour les contraintes
bvmax=vector(RR,[9/20*F,9/20*F,9/20*F,9/20*F,17/20*F,9/20*F,F,F,F,F]) #les bornes supérieures pour les contraintes
show(LatexExpr('A='),Av)
show(LatexExpr('bmin='),bvmin)
show(LatexExpr('bmax='),bvmax)
#Création du programme
v=MixedIntegerLinearProgram(maximization=True, solver='GLPK') 
#On pose les nouvelles variables allant de x_0 à x_5
ind=[0..mv-1]
x=v.new_variable(integer=False,indices=ind) 
#On utilise la fonction linéaire pour les contraintes
Bv=Av*x
#On fixe l'objectif
v.set_objective(x[4])
for i in range(mv) :
    v.set_min(x[i],0)
    v.set_max(x[i],F)
for i in range(0,nv):
    v.add_constraint(Bv[i],min=bvmin[i], max=bvmax[i])    
v.show()
v.solve()
xx=v.get_values(x)
show(xx)

But this time Sagemath return an infeasability of this program. What bother me is that if I set $x_4$ to the minimum value of the former program and use the optimal solution in this second program, I have a perfectly feasible solution.

I do not know where is the problem : a misunderstanding, Glpk (normaly it's a sure program.....). I am not certain that is a true Sagemath question...

I will add that I have changed the objective of the first program from x[0] to x[1], x[2] and x[3]. In all cases I find a solution. Of course in all those solutions the min of x[i] is always the same. But I am not sure that can explain anything.

edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
1

answered 2021-12-10 19:20:01 +0100

Your program is infeasible.

  • you have constraint 0.0 <= x_0 - 1.5 x_1 - 1.5 x_3 <= 4.5. Since x_1 and x_3 are at least 1.5 and x_0 at most 4.5, you need to set x_0=4.5, x_1=x_3=1.5 to satisfy the lower bound.
  • you also have constraint x_0 + x_1 + x_2 + x_3 == 8.5 which can be satisfied only if x_2 = 8.5 - 4.5 - 1.5 - 1.5 = 1. But since you have 2.75 <= x_2 <= 4.5, this is impossible.

Without the lower bound on constraint 0.0 <= x_0 - 1.5 x_1 - 1.5 x_3 <= 4.5, you can get a solution.

Let me rewrite the code in a way I'm more familiar with:

def foo():
    F = 10
    a = 1.5
    p = MixedIntegerLinearProgram(maximization=True)
    x = p.new_variable(nonnegative=True)
    p.add_constraint(x[0], min=3/20*F, max=9/20*F)
    p.add_constraint(x[1], min=3/20*F, max=9/20*F)
    p.add_constraint(x[2], min=11/40*F, max=9/20*F)
    p.add_constraint(x[3], min=3/20*F, max=9/20*F)
    p.add_constraint(x[0] + x[1] + x[2] + x[3] == 17/20*F)
    p.add_constraint(x[0] - a*x[1] - a*x[3] <= 9/20*F)  # upper bound only
    for i in range(4):
        p.add_constraint(x[4] <= x[i])
    p.set_objective(x[4])
    p.show()
    p.solve()
    x_val = p.get_values(x)
    show(x_val)

Then you get:

sage: foo()
Maximization:
  x_4 

Constraints:
  1.5 <= x_0 <= 4.5
  1.5 <= x_1 <= 4.5
  2.75 <= x_2 <= 4.5
  1.5 <= x_3 <= 4.5
  8.5 <= x_0 + x_1 + x_2 + x_3 <= 8.5
  x_0 - 1.5 x_1 - 1.5 x_3 <= 4.5
  - x_0 + x_4 <= 0.0
  - x_1 + x_4 <= 0.0
  - x_2 + x_4 <= 0.0
  - x_3 + x_4 <= 0.0
Variables:
  x_0 is a continuous variable (min=0.0, max=+oo)
  x_1 is a continuous variable (min=0.0, max=+oo)
  x_2 is a continuous variable (min=0.0, max=+oo)
  x_3 is a continuous variable (min=0.0, max=+oo)
  x_4 is a continuous variable (min=0.0, max=+oo)
{0: 1.916666666666667,
 1: 1.916666666666667,
 2: 2.749999999999999,
 3: 1.916666666666667,
 4: 1.916666666666667}
edit flag offensive delete link more

Comments

David thanks for your answer. Unfortunatelly I cannot bypass the lower bound. In one way, this is not a true problem. Many problem are infeasible and must be recast. Starting from your observation my true problem is why is there a solution of the ismple maximisation of one x[i]. But what bother me is that I find a solution for each other problems.

I have also a remarks. Infeasibility is not a computational error It's a property of the system. So the expected return should not be an error but simply "The probelm is infeasible".

CyrilleP gravatar imageCyrilleP ( 2021-12-12 09:51:38 +0100 )edit

You can at least use try except statements to catch infeasibility (may be not the best solution).

from sage.numerical.mip import MIPSolverException
try:
    p.solve()
except MIPSolverException as msg:
    msg = str(msg)
    if "no feasible solution" in msg:
        # do something
    else:
        # do something else
David Coudert gravatar imageDavid Coudert ( 2021-12-12 19:14:55 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2021-12-10 09:52:22 +0100

Seen: 338 times

Last updated: Dec 10 '21