0

I would like to solve a system of linear equations with the added constraint that the solution vector must only contain integers. I assume this can be accomplished by something like intlinprog, but it is not clear to me how a regular linear equation is supposed to be translated into a linear programming problem that this function can accept. Maybe the solution is obvious, maybe I'm just really confused about how to approach such a problem. Could somebody point me in the right direction?

As an example:

n = 7;

A = magic(n);

x = (1:n)' + (2 * rand(n, 1) - 1) / 10;

y = A * x;

clear x

How do I reconstruct x from just A and y so that x comes out as [1; 2; 3; 4; 5; 6; 7]?

Marcel Tesch
  • 184
  • 3
  • 15
  • This should be straightforward. Use f=zeros(n). – Erwin Kalvelagen Dec 03 '20 at 02:31
  • @ErwinKalvelagen How does this solve the problem? Do you mean like `intlinprog(zeros(n, 1), 1:n, A, x);`? `x` comes out as all zeros this way or am I misunderstanding something? – Marcel Tesch Dec 03 '20 at 09:19
  • You need to read the documentation a bit more carefully. You want to solve Ax=y for x (A and y are data). That means you need to pass on `Aeq=A,beq=y` (and not `A,b` which is for inequalities). Again, the documentation is quite good so it pays off to read it.. – Erwin Kalvelagen Dec 03 '20 at 09:44
  • @ErwinKalvelagen I did read the documentation and passing `A, y` as `Aeq, beq` was the first thing I tried, but this does not work. `intlinprog(zeros(n, 1), 1:n, [], [], A, x);` just gives `No feasible solution found. Intlinprog stopped because no integer points satisfy the constraints.` – Marcel Tesch Dec 03 '20 at 10:04
  • Why do you pass on x? You are solving for x. You need to pass on y. Again, to repeat myself: You want to solve Ax=y for x (A and y are data). – Erwin Kalvelagen Dec 03 '20 at 10:11
  • @ErwinKalvelagen Sorry, I mistyped. I did actually pass `y` as `beq`. Same result, no solution. – Marcel Tesch Dec 03 '20 at 10:13
  • 1
    Yes, sorry. I missed your fractional perturbation of x. I don't think there is a foolproof way to recover the original integer solution. Just do A\y and round would be a starting point. – Erwin Kalvelagen Dec 03 '20 at 10:39
  • 1
    A MIP model can look like: `min ||r||, subject to Ax=y+r` for some norm `||.||` A 1-norm allows it to keep it linear (using standard absolution value reformulations). E.g. `min sum(rplus+rmin), subject to Ax=y+rplus-rmin, rplus>=0, rmin>=0`. – Erwin Kalvelagen Dec 03 '20 at 10:49
  • 1
    A different MIP model would be: `min ||r||, subject to A(x+r)=y` Not sure which one is more appropriate. – Erwin Kalvelagen Dec 03 '20 at 11:37
  • @ErwinKalvelagen Thanks, those options look promising. Could you elaborate a little more on how something like this could get implemented in Matlab? – Marcel Tesch Dec 03 '20 at 13:34
  • 1
    You need to form Aeq=beq by assembly of submatrices. Ie. `Aeq=[A -I I]`.The vector to variables is `X = [x; vplus; vmin]` and `f = [0; 1; 1]` (with 1 a vector of ones of length n). – Erwin Kalvelagen Dec 04 '20 at 11:37
  • @ErwinKalvelagen Thanks! This means the problem would need to be formulated something like `intlinprog(f, 1:n, [], [], Aeq, y);` with `f = [zeros(n, 1); ones(2 * n, 1)];` and `[A, -eye(n), eye(n)]`? – Marcel Tesch Dec 04 '20 at 12:24
  • @ErwinKalvelagen That results in `Root LP problem is unbounded.` however. I supposed that is because there is no constraint pushing `vplus` and `vmin` towards zero, so they are free to diverge to `-Inf`? – Marcel Tesch Dec 04 '20 at 12:25
  • How can `min sum(rplus+rmin), rplus>=0, rmin>=0` be unbounded??? – Erwin Kalvelagen Dec 04 '20 at 12:30
  • @ErwinKalvelagen Ah, sorry. I forgot to add the inequality constraint that `rplus` and `rmin` need to be positive. Now it's `intlinprog(f, 1:n, Aineq, bineq, Aeq, y);` with `Aineq = [zeros(2 * n, n), -eye(2 * n)];` and `bineq = zeros(2 * n, 1);`. Does this look about right? It does now correctly reconstruct `x` as `(1:n)'`. – Marcel Tesch Dec 04 '20 at 12:45
  • 1
    These inequalities are bounds. I would use LB. – Erwin Kalvelagen Dec 04 '20 at 12:48
  • @ErwinKalvelagen Thanks, that makes it even easier. So it's `intlinprog(f, 1:n, [], [], Aeq, y, lb);` with `lb = [-Inf(n, 1); zeros(2 * n, 1)];` now. – Marcel Tesch Dec 04 '20 at 12:54
  • @ErwinKalvelagen Again, thank you very much for your help, I learned a lot. – Marcel Tesch Dec 04 '20 at 12:54

0 Answers0