I have a moderate-sized LP with ~13k variables and ~25k constraints. When I try to solve using MOSEK, it says that the problem is primal infeasible. By construction, I know a feasible solution exists. Also, when I set the objective function to a constant and then solve the problem, it turns out to be optimal. Furthermore, both problems seems to be optimal when I use KNITRO.

How and why it is happening? Any help would be greatly appreciated.

Below is the code and attached is the GDX file. "model sigma" is infeasible and "model dummy_problem" is optimal.

Code: Select all

```
$gdxin divided.gdx
Sets
n1
n2
m2
m3
m_ineq
m_eq
m_coupeq
;
$load n1 n2 m2 m3 m_ineq m_eq m_coupeq
Parameters
c1_sigma(n1)
f1_sigma(n2)
c1_sch(n1)
f1_sch(n2)
Aineq(m_ineq,n1)
bineq1(m_ineq)
Aeq(m_eq,n1)
beq1(m_eq)
Acoupeq(m_coupeq, n1)
Bcoupeq(m_coupeq,n2)
bcoupeq1(m_coupeq)
C(m2,n2)
d1(m2)
F(m3,n2)
g1(m3)
a_i
d_i
x_feas(n1)
y_feas(n2)
$load c1_sigma f1_sigma c1_sch f1_sch C d1 F g1 a_i d_i Aineq bineq1 Aeq beq1 Acoupeq Bcoupeq bcoupeq1 x_feas y_feas
;
Variables
obj_sigma
x(n1)
;
Positive Variable
y(n2)
;
y.fx(n2)$(n2.pos<= a_i) = 0;
y.fx(n2)$(n2.pos > d_i and n2.pos <= 96) = 0;
y.fx(n2)$(n2.pos > 96 and n2.pos <= 96 + a_i) = d1('96');
Equations
const1_eq(m_eq)
const1_ineq(m_ineq)
const1_coupeq(m_coupeq)
const2(m2)
const3(m3)
objeq_sigma
objeq_dummy
;
obj_sigma.up = 0;
obj_sigma.lo = -1;
const1_eq(m_eq).. sum(n1, Aeq(m_eq,n1)*x(n1)) =e= beq1(m_eq);
const1_ineq(m_ineq) .. sum(n1, Aineq(m_ineq,n1)*x(n1)) =g= bineq1(m_ineq);
const1_coupeq(m_coupeq) .. sum(n1, Acoupeq(m_coupeq,n1)*x(n1)) + sum(n2, Bcoupeq(m_coupeq,n2)*y(n2)) =e= bcoupeq1(m_coupeq);
const2(m2).. sum(n2, C(m2,n2)*y(n2)) =e= d1(m2);
const3(m3).. sum(n2, F(m3,n2)*y(n2)) =g= g1(m3);
objeq_sigma.. obj_sigma =e= sum(n1, c1_sigma(n1)*x(n1)) + sum(n2, f1_sigma(n2)*y(n2));
objeq_dummy.. obj_sigma =e= -0.5;
option lp=mosek;
Model dummy_problem / const1_eq, const1_ineq, const1_coupeq, const2, const3, objeq_dummy /;
Model sigma / const1_eq, const1_ineq, const1_coupeq, const2, const3, objeq_sigma /;
sigma.OptFile = 1;
dummy_problem.OptFile = 1;
*Solve dummy_problem using lp minimizing obj_sigma ;
Solve sigma using lp minimizing obj_sigma ;
*display sigma.modelstat;
```