solve: Cannot find explicit solution but it have solutions, matlab's bug?

this is my code:

syms x0 x1 x2 x3 x4;

cond1 = (1.00028*x0+1.00045*x1)>=1.00035*(x0+x1);
cond2 = ((1.00028*x0+1.00045*x1) + 1.00063*x2)>=1.00054*(x0+x1+x2);
cond3 = (((1.00028*x0+1.00045*x1) + 1.00063*x2)+1.00082*x3)>=1.00076*(x0+x1+x2+x3);
cond4 = ((((1.00028*x0+1.00045*x1) + 1.00063*x2)+1.00082*x3)+1.00101*x4)>=1.00091*(x0+x1+x2+x3+x4);

cond5 = x0+x1+x2+x3+x4>=0;

condconds = [cond1 cond2 cond3 cond4 cond5];
sol = solve(condconds, [x0 x1 x2 x3 x4])

Does my code have any error or matlab have bugs???

My solution :

x0 1

x1 1

x2 3.44

x3 19.94

x4 38.07

by the way, i just need integer x

[Sara Boznik]

Your code looks fine. On my computer either does not work.

[John D'Errico]

Actually, no. The code is not fine, at least not in terms of what MATLAB allows solve to handle. Solve is not designed to solve that class of probem. The problem has infinitely many solutions, and solve cannot solve problems with infinitely many solutions happily.

[andrew wu]

So can you help me which function would solve my problem? I just need a few solutions

[John D'Errico]

I showed exactly how to solve that problem in my answer! I guess I thought it would be clear that since your problem is linear and the system effectively homogeneous, then you can trivially scale that solution by any constant. I'll concede that may not have been obvious.

X = linprog(ones(5,1),-DA,zeros(5,1),[],[],ones(5,1))
Optimal solution found.
X =
                         1
                         1
          3.88888888889285
          21.5925925926071
          41.2222222222452

Thus, if X is this vector, than ANY multiple of X is also a solution. This must be true, since your equations form a homogeneous system.

syms t
vpa(subs(condconds,[x0,x1,x2,x3,x4],t*X.'),10)
ans =
[ 2.0007*t <= 2.00073*t, 5.892068889*t <= 5.892068889*t, 27.50236741*t <= 27.50236741*t, 68.76622407*t <= 68.76622407*t, 0.0 <= 68.7037037*t]

The vector X*t is clearly a solution for any positive value of t.

[andrew wu]

I dont understand your code, it dont include my data, I changed my question that just need integer solutions

[andrew wu]

I run your code but it shows undefined variable DA

Answer 1

[John D'Errico]

Look at what you have. A linear system of inequalities. Solve is designed to solve equalities. The solution that you have may be a solution. However, there are infinitely many solutions. One of those solutions must be the zero solution, since your problem is a homogeneous one, with no constant terms.

But, is your "solution truly a solution? I'm afraid not.

vpa(subs(condconds,[x0,x1,x2,x3,x4],[1,1,3.44,19.94,38.07]),10)
ans =
[ 2.0007 <= 2.00073, 5.4429376 <= 5.4428972, 25.3992888 <= 25.399248, 63.5077395 <= 63.5076987, 0.0 <= 63.45]

So the 2nd, 3rd and 4th constraints are all invalidated by your solution, though the errors are small. So I'll ssume you have merely rounded off what may be a valid solution.

Now let me look at the equality form of your constraint system that Stefan found. Is that actually rank deficient?

svd(double(A))
ans =
       2.2361
   0.00074952
   0.00018996
   0.00011237
   9.3138e-05

NO. Not at all.

Regardless, if we return to the inequality constrained form of your problem, can we find solutions? Probably. A common trick is to use linear programming.

linprog wants <= inequalities, so I'll swap the signs on A.

DA = double(A);
linprog(ones(5,1),-DA,zeros(5,1))
Optimal solution found.
ans =
     0
     0
     0
     0
     0
     
linprog(-ones(5,1),-DA,zeros(5,1))
Problem is unbounded.
ans =
     []

So we do find a solution as the zero vector, as we must, since A is actually non-singular. But by looking in another direction, linprog tells us the problem is unbounded. So there are infinitely many solutions in some direction.

Is there a solution where all of the knowns are no smaller than 1?

format long g
>> X = linprog(rand(5,1),-DA,zeros(5,1),[],[],ones(5,1))

Optimal solution found.

X =
                         1
                         1
          3.88888888889285
          21.5925925926071
          41.2222222222452

Interesting, in that linprog finds one near what you called the solution. As I said though, there will be infinitely many solutions, but perhaps this will be good enough for you.

[andrew wu]

Although you don't completely solve my problem, you give me the function "linprog()" that bring me a lot about thoughts, then I use "intlinprog()" to resolve the problem. your answer is best among the answers. Thank you for your help.

Answer 2

[Stephan]

Congratulations, you have 1 solution of an infinite number of possible solutions - look to the borderline case of equal to:

syms x0 x1 x2 x3 x4;
assume([x0 x1 x2 x3 x4], {'real', 'positive'})
cond1 = (1.00028*x0+1.00045*x1)==1.00035*(x0+x1);
cond2 = ((1.00028*x0+1.00045*x1) + 1.00063*x2)==1.00054*(x0+x1+x2);
cond3 = (((1.00028*x0+1.00045*x1) + 1.00063*x2)+1.00082*x3)==1.00076*(x0+x1+x2+x3);
cond4 = ((((1.00028*x0+1.00045*x1) + 1.00063*x2)+1.00082*x3)+1.00101*x4)==1.00091*(x0+x1+x2+x3+x4);
cond5 = x0+x1+x2+x3+x4==0;
condconds = [cond1 cond2 cond3 cond4 cond5];


[A,~] = equationsToMatrix(condconds);
det_A = vpa(det(A))

leads to

det_A =
 
A =
 
[                              -7/100000,                               1/10000,                                      0,                                      0,                             0]
[                              -13/50000,                             -9/100000,  1266637395196663/14073748835532800000,                                      0,                             0]
[                                -3/6250,                            -31/100000, -1829587348620553/14073748835532800000,                                3/50000,                             0]
[ -8866461766385191/14073748835532800000, -1294784892868923/2814749767106560000,           -78812993479/281474976710656, -1266637395197479/14073748835532800000, 450359962737/4503599627370496]
[                                      1,                                     1,                                      1,                                      1,                             1]
 
 
b =
 
 0
 0
 0
 0
 0
 
 
det_A =
 
0.0000000000000033319999999983153682592275926204

which means, that the determinant of your system is zero in fact. Since solve is looking for a explicit solution it will not find any i guess.

[John D'Errico]

I'm sorry, but this is incorrect.

svd(double(A))
ans =
       2.2361
   0.00074952
   0.00018996
   0.00011237
   9.3138e-05

A is not even remotely close to being a singular matrix. The condition number is not even that terribly large.

cond(double(A))
ans =
        24008

Using the determinant to test for singularity is about as poor a solution as any I would ever recommend. Would you tell me that

B = eye(5)*0.001;

is a singular matrix? det seems to think it is, yet you should agree it is not.

det(B)
ans =
        1e-15

NEVER use det to test a matrix for singularity, unless you are solving a homework assignment as directed. And even then, you should tell your teacher why it is such a bad idea to do so. Or let me tell them why.

[Stephan]

My first idea was to delete my answer, but I (and may be other people) can learn from your comment. Instead of deleting my incorrect answer, i vote for yours.

;-)

[andrew wu]

Thank you