Using solve with inequalities

I would like to use Solve to check if there are solutions to various systems of equalities/inequalities. Here's a toy example:

syms x y
S = solve(x+y==1,x<0.4,y<0.4,[x,y]);

Obviously there's no solution to that one, so Matlab gives me empty S.x and S.y. Interestingly, when I do the following, even though there are obviously a continuum of solutions, it spits out one answer:

S = solve(x+y==1,x<0.6,y<0.6,[x,y]);

It gives S.x = S.y = 1/2, which is, in some sense, the most intuitive possible solution. I like that it gives me that. Unfortunately, when I do the following, it's more concerned about that fact that there is a continuum of solutions:

S = solve(x+y==1,x<0.6,y<0.7,[x,y]);

In this case, it gives me:

Warning: Cannot find explicit solution.

What I would like to do is to have it spit out just one possible solution -- I don't really care which it chooses. I know I can get a parametrized solution doing the following:

S = solve(x+y==1,x<0.6,y<0.7,[x,y],'ReturnConditions', true)

But then the mathworks page for Solve suggests how I can then propose a solution and have Matlab check whether it meets those conditions. But it doesn't offer advice on how to have Matlab propose a solution. Is there a way to do so?

Incidentally, I've also considered doing this numerically with linprog, which seems a little faster. The issue I have is with the tolerance on the constraints. Since I have strict inequalities, I've been making them weak inequalities by saying a <= b+1e-4 is roughly equivalent to a<b. However, this causes some issues for me (where it fails to find solutions where some exist). I'd really need to be able to set the constraint tolerances to zero to do it that way I think, and I don't believe that's possible.

Any advice (on how to do this numerically using linprog or symbolically) would be greatly appreciated.

Answer 1

[Matt J]

My suggestion would be to use LCON2VERT ( Download ) if your system of in/equalities isn't too large (it couldn't be, I guess, if you were hoping to do the analysis symbolically) and if they define a bounded region.

LCON2VERT will give you the vertices of the polyhedral region corresponding to the in/equalities. If you just want an arbitrary point that strictly satisfies the inequalities, just take the mean of all the vertices.

As an example, instead of

   S = solve(x+y==1,x<0.6,y<0.7,[x,y]);

you can do this,

    >> Aeq=[1 1]; beq=1; A=eye(2); b=[0.6;0.7];

    >> vertices=lcon2vert(A,b,Aeq,beq)

    vertices =

        0.6000    0.4000
        0.3000    0.7000

    >> mean(vertices)

    ans =

        0.4500    0.5500

[Shane]

Ooh - that's interesting. I'll give that a go tomorrow and update. Thanks Matt.

[Shane]

Sorry for the delay, Matt. I got distracted on other projects. This does work for the application I'm working on. Thanks for bringing LCON2VERT to my attention.

[R. Can]

Dear Matt, Thank you for this answer. I think is a very useful function. Is it possible to modify this routine LCON2VERT If I always want to find solution for systems with strict inequality and without allowing the equality:

Aeq*x=beq

with A*x<b ?

Thank you again.

All best

[Matt J]

Dear R. Can,

If no equality constraints are present, you can follow the same approach as in my answer above with the Aeq, beq arguments to lcon2vert omitted,

vertices=lcon2vert(A,b)

Or, you can pass [] arguments,

vertices=lcon2vert(A,b,[],[])

[R. Can]

Dear Matt,

Thank you. Forgive me for miss formulating the question. I want to solve Aeq*x=beq and A*x<b at the same time. The unexisting equality I was refering to is that the program specifies that Aeq*x=beq and A*x<=b with A*x lower-equal to b. What I would want to have strick is that I'm solving "Aeq*x=beq and A*x strictly lower than b".

If I pass lcon2vert(A,b,[],[]) then I am solving only A*x<=b with the condition lower-equal to. Still I wouldn't be solving a strict lower to A*x lower-strictly-than b.

Hope I manage to transmit my question clearly this time.

Thank you again

[Matt J]

Dear R.Can,

Your question is clear now, but I think it is the same question as originally posted by Shane, and so my answer will be the same. By taking the mean of the vertices output by lcon2vert, you will obtain a solution strictly satisfying the inequalities.

[R. Can]

Thank you again, Matt.