Infeasibility due to normalization constraint?
AnsweredI am implementing Nanongkai et al.'s LP for maximum regret computation of a point p, given a set of points S. Find the LP formulation below:
max x
s.t. (p − p ′ ) · v ≥ x ∀ p ′ ∈ S
p · v = 1
v[j] ≥ 0 ∀ j ≤ d
x ≥ 0
Nanongkai et al. implemented their LP on GLPK, introducing further variables r1, r2, y_i:
max x
s.t. x − (p − pi) · v ≤ yi ∀ i ≤ t
p · v ≤ r1
−p · v ≤ r2
0 ≤ v[j] < ∞ ∀ j ≤ d
0 ≤ x ≤ ∞
−∞ < yi ≤ 0 ∀ i ≤ t
−∞ < r1 ≤ 1
−∞ < r2 ≤ −1
First of all, to ensure I understand the purpose of r1 and r2, they aim to prevent the model from becoming infeasible because perfect equality may not be achieved due to floating-point errors. Is that interpretation correct?
My current model becomes infeasible for any tuple with regret x = 0.0. I would appreciate it if someone could review my implementation below and identify what is incorrect about it. I have tried several formulations, including formulating the normalisation constraint using GRB.EQUAL, introducing slack and setting r1, r2 to fixed values 1 and -1, respectively.
I am still new to Linear Programming (continuous models in particular, the problems we solved at university were all using binary variables). Thanks in advance for any input!
public static double solveMaxRegret(GRBEnv env, Point p, List<Point> S) throws GRBException {
double[] pCoords = p.getCoordinates();
int d = pCoords.length;
int t = S.size();
GRBModel model = new GRBModel(env);
// regret variable x >= 0 (to be maximized)
GRBVar x = model.addVar(0.0, GRB.INFINITY, 1.0, GRB.CONTINUOUS, "x");
// utility vector v[j] >= 0, d >= j >= 0
GRBVar[] v = new GRBVar[d];
for (int j = 0; j < d; j++) {
v[j] = model.addVar(0.0, GRB.INFINITY, 0.0, GRB.CONTINUOUS, "v_" + j);
}
GRBVar[] y = new GRBVar[t];
// constraint: x - (p-p_i) * v <= y_i
for (int i = 0; i < t; i++) {
double[] pi = S.get(i).getCoordinates();
y[i] = model.addVar(Double.NEGATIVE_INFINITY, 0.0, 0.0, GRB.CONTINUOUS, "y_" + i);
GRBLinExpr expr = new GRBLinExpr();
// Sum over j: (p_j - p_i[j]) * v_j
for (int j = 0; j < d; j++) {
expr.addTerm(pCoords[j] - pi[j], v[j]);
}
// Add the -x term
expr.addTerm(-1.0, x);
// Add y_i term
expr.addTerm(1.0, y[i]);
// Constraint: expr >= 0 <=> x - (p-p_i).v <= y_i
model.addConstr(expr, GRB.GREATER_EQUAL, 0.0, "regret_" + i);
}
// normalization: p * v = 1
GRBLinExpr pDotV = new GRBLinExpr();
for (int j = 0; j < d; j++) {
pDotV.addTerm(pCoords[j], v[j]);
}
GRBVar r1 = model.addVar(Double.NEGATIVE_INFINITY, 1.0, 0.0, GRB.CONTINUOUS, "r1");
GRBVar r2 = model.addVar(Double.NEGATIVE_INFINITY, -1.0, 0.0, GRB.CONTINUOUS, "r2");
// constraint: p * v <= r1
model.addConstr(pDotV, GRB.LESS_EQUAL, r1, "pDotV_leq_r1");
GRBLinExpr negPDotV = new GRBLinExpr();
for (int j = 0; j < d; j++) {
negPDotV.addTerm(-pCoords[j], v[j]);
}
// constraint: -p * v <= r2
model.addConstr(negPDotV, GRB.LESS_EQUAL, r2, "neg_pDotV_leq_r2");
// objective: maximize x
GRBLinExpr obj = new GRBLinExpr();
obj.addTerm(1.0, x);
model.setObjective(obj, GRB.MAXIMIZE);
model.optimize();
double regret;
int status = model.get(GRB.IntAttr.Status);
if (status == GRB.Status.OPTIMAL) {
regret = x.get(GRB.DoubleAttr.X);
} else {
System.out.println(model.get(GRB.IntAttr.Status));
regret = Double.NaN;
}
model.dispose();
return regret;
}
-
-
Gurobi Staff
Hi Marian,
Gurobi does not have to satisfy every constraint exactly, but only up to some tolerance (see FeasibilityTol). That way, Gurobi prevents minor rounding errors from causing infeasibility. In fact, Gurobi will simplify the second formulation to the first one in its presolve step. So, it shouldn't matter to Gurobi which of those two formulations you choose.
I also checked your code, and as far as I can tell, you implemented everything correctly. I even tried it with a couple of input values, and Gurobi found optimal solutions for those examples.
However, I wonder if the formulation itself has a potential problem. Is there anything preventing (p[j] - p'[j]) from turning into a negative number? If not, then it is relatively easy to construct infeasible input values, e.g., a set S that contains a point p' s.t. p[i] < p'[j] at every position j. That's why my gut says that (p[j] - p'[j]) should be |p[j] - p'[j]|, but I don't know enough about regret computation to say for sure.
Could you post a set of input values that causes infeasibility? Then I or someone else can potentially tell you more.
The following article is helpful if you want to analyze the infeasibility yourself: How do I determine why my model is infeasible?
Best regards,
Martin
0 -
-
Hi Martin,
Thank you very much for looking into this, and thanks for the useful input regarding Gurobi's behaviour to prevent rounding errors. I have taken a course at university about linear programming using Gurobi, but we restricted ourselves mainly to binary variables and things have been easier that way… :-).
What I found out in the meantime: The LP only becomes infeasible for points which are entirely dominated by the current solution set S (i.e. their regret should be 0.0). As long as adding them would improve the solution's regret, the values are computed correctly. My workaround currently is catching the model status and returning regret = 0.0 if the model becomes infeasible.
The greedy algorithm I am implementing iteratively builds the solution set by calling the LP repeatedly (see Nanongkai et al. 2010, the same paper where the LP formulation is from), and with my workaround in place, it behaves as expected.
As an example, let D = {(125, 1000); (250, 800); (562.5, 625); (750, 250); (1000, 75)} be the set of tuples. We want to compute a regret-minimising set S \subset D of size k = 3 using GREEDY.
First of all, GREEDY adds (1000, 75) to the solution because it maximises the first dimension. The remaining two points are selected by calling the LP for each candidate in D\S and choosing the one exhibiting the largest regret. As mentioned above, the regret values are returned correctly by the LP as long as their regret is larger than zero. However, after adding point (125, 1000) in the first iteration (it has the largest regret for S = {(1000, 75)}), most remaining points will be sufficiently covered, and only (562.5, 625) has regret > 0. For the remaining points, the model becomes infeasible. By setting the regret to zero in these cases, I get the desired result ((562.5, 625) is selected as the third and last point in S) - but I guess it's just a workaround for now, and I would clearly prefer a “clean” version.
If I am not mistaken, the LP is trying to search for a (normalised) utility vector v that satisfies the constraints and maximises over x. So essentially, the LP is searching for a utility vector that would result in p outperforming the points currently in S. When p is entirely dominated, however (as is the case in the second iteration in my example), such a vector does not exist (and x would remain below zero, violating the lower bound).
I'd appreciate it if you'd take another look at the provided example.Best regards,
Marian0 -
-
-
Hi Marian,
Thanks for the explanation. I think I now understand what the expected behavior should be.
In the current formulation, the constraints \( (p − p ′ ) \cdot v \geq x \) force the variable \(x\) to take a negative value if the points in \(S\) dominate the point \(p\). This leads to infeasibility since \(x\) has a lower bound of \(0\). The correct behavior would be that the points in \(S\) are allowed to dominate the point \(p\), but that the variable \(x\) becomes \(0\) in this case.
To get the expected behavior, I'd reformulate your problem as follows:
\[\begin{align*}\text{max } & x & \\ \text{s.t. } &(p − p ′ ) \cdot v \geq x' &\forall p' \in S \\ &p \cdot v = 1 & \\ &x = \text{max}(x',0) & \\ &v[j] \geq 0 &\forall j \leq d \\ &x \geq 0 & \end{align*}\]
Note that \(x'\) has no lower bound, so \(x'\) may become negative. However, since \(x\) is the maximum of \(x'\) and \(0\), we are guaranteed that \(x\) is always non-negative.
However, it might be more efficient if you take the original LP formulation, set the lower bound of x to \(-\infty\), and handle negative regrets on your side as if they were zero.
Best regards,
Martin
0 -
Please sign in to leave a comment.
Comments
3 comments