MILP modelling on minimal disturbance of right-hand-side to make a linear system infeasible

Question

I try to model the following problem: given $z\in\{0,1\}^m$ and a linear system $Ax\le b(z), x\in\Bbb R^n, A\in\Bbb R^{d\times n}, b(z)\in\Bbb R^{d}$, where $b(z)$ means that some entries of $b$ are functions of $z$. Specifically, the entries in $b$ that depend on $z$ are in the form of $c + 1-z$, where $c$ indicates some constant.

The statement I want to model is: minimize $\|z\|_0$ such that the linear system $Ax\le b(z)$ is infeasible. It is already known that this linear system is feasible when $z = 0$ and there exists a $z$ making the system infeasible. Is there any way to model it into MILP?

Answer 1

By Farkas lemma, infeasibility of $Ax\leq b$ is equivalent to feasibility of $A^Ty = 0, y^Tb < 0, y\geq 0$, or more practically useful $A^Ty=0, y^Tb \leq -1, y\geq 0$.

Unfortunately, this will lead to a bilinear model when you parameterize $b(z)$. It is fairly similar to an application I worked on a decade ago Oops! I cannot do it again: Testing for recursive feasibility in MPC where we started with bilinear stuff coming from Farkas lemma, and it was possible to trade that bilinearity with a complementarity condition (Sec 4.2) which (at the time) was more efficiently solved.

EDIT: Missed the fact that $z$ is binary. In that case, you can linearize $y^Tb(z)$ using standard big-M methods, and you end up with a MILP.

Here is a small proof-of-concept implementation in the MATLAB Toolbox YALMIP (disclaimer, developed by me). I was too lazy to manually linearize so I used a built-in command to derive the big-M model

% Feasible model
m = 50;n = 5;
A = randn(m,n);
b0 = A*randn(n,1) + 0.99

% Perturb
z = binvar(m,1);
b = b0 - z;

% Solve
y = sdpvar(m,1);
[linear_b_y,cut]=binmodel(b'*y,[-100 <= y <= 100]);
Model = [A'*y == 0, y >= 0, linear_b_y == -1, cut];
optimize(Model,sum(z))
nnz(z)

% Check that it is infeasible
x = sdpvar(n,1);
optimize(A*x<=value(b))

Answer 2

I think this can be approached using a constraint generation technique (variant of Benders decomposition), although I have no idea if it would efficient. By reordering the rows of $A$, we can assume that $$b(z)=\left[\begin{array}{c} \hat{b}\\ c+e-z \end{array}\right]$$where $\hat{b}\in\mathbb{R}^{d-m}$, $c\in \mathbb{R}^m$, $e=(1,\dots,1)^\prime \in \mathbb{R}^m$ and $z\in \lbrace 0,1 \rbrace^m$.

You start out minimizing $\parallel z\parallel_{0}$ subject only to $z\in \lbrace 0,1 \rbrace^m$. Each time the solver finds what it thinks is a new incumbent $\tilde{z}$ (satisfying all constraints added to that point), you solve the inequality $Ax\le b(\tilde{z})$. If the inequality is infeasible, accept the new incumbent and keep chugging along. If a feasible solution $\tilde{x}$ is found, let $$A\tilde{x}=\left[\begin{array}{c} \tilde{u}\\ \tilde{v} \end{array}\right],$$noting that $\tilde{u}\le \hat{b}$ (about which you can do nothing) and $\tilde{v}\le c + e - \tilde{z}$. Let $J=\lbrace j : \tilde{z}_j = 0 \land \tilde{v}_j > c_j\rbrace$. If $J=\emptyset$ you are screwed: there is no way to make $\tilde{x}$ infeasible. Otherwise, you can add the constraint $$\sum_{j\in J} z_j \ge 1$$to the IP model, making $\tilde{z}$ infeasible, and resume solving.