# Solving a system of nonlinear equations with inequality constraints

Posted 3 months ago
816 Views
|
11 Replies
|
2 Total Likes
|
 Hi! I'm pretty new to Mathematica, am I am having trouble trying to solve a system of nonlinear equations. In this example I have 5 equations, with 4 inequality constraints. I've looked through the forums, and found an example (seen as the second example) using FindInstance, so I tried to use it for my problem, but unlike the one I found on the forum, when I try to execute my command (the top one), I don't get an output. I know this series of nonlinear equations should have valid solutions because I found some solutions on MATLAB, but I'm trying to double check my answers. Also, if FindInstance isn't the right command for this, please let me know. I know about Solve and Reduce, but I don't know how to incorporate inequality constraints with those commands. Thanks in advance! Answer
11 Replies
Sort By:
Posted 3 months ago
 I'm not sure how to interpret some of your latest message. See what you can make of this.I take your latest Mathematica, modify that slightly and hope I haven't made any typos entering that. sol=Solve[{-2a1x+10/(a1x+a1y-10)+g1==0, -2a1y+10/(a1x+a1y-10)+g2==0, 0.9(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3==0, 0.9(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4==0, -1/2+0.9(a2x^2/40+a2y^2/40)+g5==0,-1/2+g6==0, g1*-a1x==0,g2*-a1y==0,g3*-a2x==0,g4*-a2y==0, g5*-ix==0,g6*-iy==0,g7*(ix-24)==0,g8*(iy-24)==0}, {a1x,a1y,a2x,a2y,ix,iy,g1,g2,g3,g4,g5,g6,g7,g8}] I'm using Solve instead of Reduce. For this problem they should be about the same, but the way it formats the answer will be different.When I give that to Mathematica it issues a warning Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result That is not an error in this case. It is just saying that it looked at your problem which had a mixture of exact integer and rational coefficients with some approximate decimal numbers 0.9. The algorithm it selected to try to solve this couldn't deal with the 0.9, so it translated those into 9/10, solved the problem and then turned everything back into approximate decimals.And it gives me 91 different solutions! Two it shows me in detail and the other 89 are hidden behind a ...89... If you need to you can append a //InputForm onto the end of that Solve line and you should see about 15 pages of solutions, but InputForm changes things and is mostly just to display something, not to do further calculations with.Now I am not going to use that InputForm and I am going to try to test the first solution to see if it works. sol[] is going to extract the first solution and /. is going to substitute that into your system and show the result. {-2a1x+10/(a1x+a1y-10)+g1==0, -2a1y+10/(a1x+a1y-10)+g2==0, 0.9(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3==0, 0.9(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4==0, -1/2+0.9(a2x^2/40+a2y^2/40)+g5==0,-1/2+g6==0, g1*-a1x==0,g2*-a1y==0,g3*-a2x==0,g4*-a2y==0, g5*-ix==0,g6*-iy==0,g7*(ix-24)==0,g8*(iy-24)==0}/.sol[] That returns {True,True,False,False,False,True,True,True,True,True,True,True,True,True} and thus every one of those equations is true with that particular solution except equations 3,4,5.What is going on with those? I'm going to look a those this way {0.9(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3, 0.9(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4, -1/2+0.9(a2x^2/40+a2y^2/40)+g5}/.sol[] and that returns {-9.99201*^-17,-9.99201*^-17,1.11022*^-16} Ah, because of the inexact approximate decimal values each of those tiny floating point numbers would have been compared with zero, they were not exactly equal to zero and so that is why it reported those three equations were False. Notice the three equations it complained about were exactly the three equations with approximate decimal 0.9 in them. If I turn those three 0.9 into 9/10 and do this all over again then the solutions get bigger and more complicated and involve more square roots and powers and etc, but it looks like the warning and the False go away, but I have not checked this in every case.Learning how to understand and deal with tiny floating point discrepancies is a whole different world and we probably don't have time to even begin exploring that at this point.You can test others of your 91 solutions by doing things like {-2a1x+10/(a1x+a1y-10)+g1==0, -2a1y+10/(a1x+a1y-10)+g2==0, 0.9(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3==0, 0.9(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4==0, -1/2+0.9(a2x^2/40+a2y^2/40)+g5==0,-1/2+g6==0, g1*-a1x==0,g2*-a1y==0,g3*-a2x==0,g4*-a2y==0, g5*-ix==0,g6*-iy==0,g7*(ix-24)==0,g8*(iy-24)==0}/.sol[] which returns {True,True,True,False,False,True,True,True,True,True,True,True,True,True} So I suspect that most of the solutions are "close", but not exactly equal due to the small uncertainty involved with any floating point number.I realize that Matlab lives on approximate decimal numbers and can handle exact integers and rationals, often turning those into approximate decimals. Mathematica tends to go the other way and often turns things into exact rationals when possible.So, where are we now? I'm not certain. Is there something specific you need from this? Can you come up with a specific case or example that you need to test or understand?As I always say, check all of this very carefully before you even begin thinking you might trust something that I've done. Answer
Posted 3 months ago
 Using the Solve code in your most recent message:Each solution is enclosed in {} and all the solutions are a list enclosed in another layer of {} so I can use Length[sol] to find the length of that list and it tells me there are 91 different solutions.Next, there is nothing in the Mathematica code of your most recent message that says a2x must be greater than or equal to zero. If you want to include that condition then this sol=Solve[{-2a1x+10/(a1x+a1y-10)+g1==0, -2a1y+10/(a1x+a1y-10)+g2==0, 9/10(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3==0, 9/10(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4==0, -1/2+9/10(a2x^2/40+a2y^2/40)+g5==0,-1/2+g6==0, g1*-a1x==0,g2*-a1y==0,g3*-a2x==0,g4*-a2y==0, g5*-ix==0,g6*-iy==0,g7*(ix-24)==0,g8*(iy-24)==0, a2x>=0}, {a1x,a1y,a2x,a2y,ix,iy,g1,g2,g3,g4,g5,g6,g7,g8}]; Length[sol] reports that there are now only 63 solutions.If I go back to an earlier message of yours and use all the conditions that it did and hope that this is the right thing to do then sol=Solve[{-2a1x+10/(a1x+a1y-10)+g1==0, -2a1y+10/(a1x+a1y-10)+g2==0, 9/10(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3==0, 9/10(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4==0, -1/2+9/10(a2x^2/40+a2y^2/40)+g5==0,-1/2+g6==0, g1*-a1x==0,g2*-a1y==0,g3*-a2x==0,g4*-a2y==0, g5*-ix==0,g6*-iy==0,g7*(ix-24)==0,g8*(iy-24)==0, a1x>=0,a1y>=0,a2x>=0,a2y>=0,0<=ix<=24,0<=iy<=24, g1>=0,g2>=0,g3>=0,g4>=0,g5>=0,g6>=0,g7>=0,g8>=0}, {a1x,a1y,a2x,a2y,ix,iy,g1,g2,g3,g4,g5,g6,g7,g8}]; Length[sol] tells me that there are only two solutions and thus sol shows me {{a1x->(5 + Sqrt)/2, a1y->(5 + 5*(5 + Sqrt))/(5 + Sqrt), a2x->0, a2y->0, ix->0, iy->0, g1->0, g2->0, g3->9/10, g4->9/10, g5->1/2, g6->1/2, g7->0, g8->0}, {a1x->0, a1y->0, a2x->0, a2y->0, ix->0, iy->0, g1->1, g2->1, g3->9/10, g4->9/10, g5->1/2, g6->1/2, g7->0, g8->0}} as the only two solutions. Note: It might not be important, but if you Simplify that value for a1y in the first solution it turns out to be the same as the value for a1x, it just looks different because some functions don't take the extra time to do simplification that you may not need or want done.If on the other hand the only constraint you want to put on the solutions is that ix==24 then include that, as the only additional equation. sol=Solve[{-2a1x+10/(a1x+a1y-10)+g1==0, -2a1y+10/(a1x+a1y-10)+g2==0, 9/10(-2a2x(1-ix/40)+10/(a2x+a2y-10))+g3==0, 9/10(-2a2y(1-ix/40)+10/(a2x+a2y-10))+g4==0, -1/2+9/10(a2x^2/40+a2y^2/40)+g5==0,-1/2+g6==0, g1*-a1x==0,g2*-a1y==0,g3*-a2x==0,g4*-a2y==0, g5*-ix==0,g6*-iy==0,g7*(ix-24)==0,g8*(iy-24)==0, ix==24}, {a1x,a1y,a2x,a2y,ix,iy,g1,g2,g3,g4,g5,g6,g7,g8}]; Length[sol] and that returns zero, because sol is {} because apparently there are no solutions when ix==24.So, is what I've shown you thus far the solution to your problem? Or at least showing you the method that I use to attack problems like this using the Mathematica way of thinking?As for "simultaneously optimize two objective functions that are dependent on each other", I don't have any confidence that I know how to translate that phrase into a line or two of Mathematica code. If you can show those two functions written in Mathematica notation and you can describe what your optimization model is then perhaps I can try to understand that and see if I can translate it into a few lines of Mathematica code.Is what I am doing with this helping you at all?Is there any way you could do what you are trying to do entirely within Matlab? Is there any chance that working within one tool and language and way of thinking might be easier for you than trying to translate back and forth between two tools and languages and ways of thinking? Answer
Posted 3 months ago
 This Reduce[{-2a1x+10/(a1x+a1y-10)==0, -2a1y+10/(a1x+a1y-10)==0, 9/10(-2a2x(1-ix/40)+10/(a2x+a2y-10))==0, 9/10(-2a2y(1-ix/40)+10/(a2x+a2y-10))==0, -1/2+9/10(a2x^2/40+a2y^2/40)==0, 0<=ix<=24,0<=iy<=24},{a1x,a1y,a2x,a2y,ix,iy}] returns False indicating no solution.This NMinimize[{(-2a1x+10/(a1x+a1y-10))^2+ (-2a1y+10/(a1x+a1y-10))^2+ (9/10(-2a2x(1-ix/40)+10/(a2x+a2y-10)))^2+ (9/10(-2a2y(1-ix/40)+10/(a2x+a2y-10)))^2+ (-1/2+9/10(a2x^2/40+a2y^2/40))^2, 0<=ix<=24,0<=iy<=24},{a1x,a1y,a2x,a2y,ix,iy}] returns {0.202798,{a1x->-0.45804,a1y->-0.45804,a2x->-1.06586,a2y->-1.06586,ix->24.,iy->10.0015}} indicating it didn't find values that resulted in a zero sum of squares.If you can show one or more of the solutions you have from elsewhere then that might help track down the contradiction.Please check all this very carefully to make certain that I didn't make any mistakes. Answer
Posted 3 months ago
 Hi Bill,Thanks for the quick reply. For reference, I used the GlobalSearch tool in MATLAB to find the solutions for this system of nonlinear equations. I used a pretty different method, because I did a global optimization of my original function using upper and lower bounds. Since the notation in MATLAB is a little different, a1x=x(1), a1y=x(2), a2x=x(3), a2y=x(4), ix=x(5), iy=x(6). As you can see I did get a set of solutions here. Since part of my project involves solving strategic actions with two players and two objective functions, I can't continue using the tools that rely on just maximizing or minimizing just one objective function, I need a method that approaches the problem by solving partial derivatives. I tried again, and this time made a Lagrangian and used the Reduce command like you did. That seemed to get me somewhere, but I'm still not getting the exact same answers I did in MATLAB. And here's a couple of pictures making sense of the Lagrangian. Based on the MATLAB results, it seems like g6 and g7 should be tight constraints, since I get ix=24 and iy=6, so I should have made g6>0 and g7>0 in my string of equations, but when I did that, I got false as my output.  Thank you again! I really appreciate your help! Answer
Posted 3 months ago
 Hi Bill, I copied your format of using Solve, but this time used 9/10 instead of 0.9 and -1/2 instead of -0.5, and also got a lot of results (though where does it tell you that there are 91 solutions?). I was hoping that using this method would allow me to go through each solutions and check that I get "True" for each variable, but when I do it for the first solutions, I get "True" for all variables even though it shouldn't be right. Is this because I've changed from the decimals to the fractions? Even so, the first set of solutions can't be right because it makes a2x<0. (You'll have to scroll down all the way for the code). Is there any way to sort through all these solutions without going through them individually? For example, what if I only want to find the solutions where ix --> 24? As for what I need from this: in this specific case, I'm trying to find the variables that optimize a function by solving the partial derivatives (subject to a few constraints). I technically already know what the solutions should be (if MATLAB is to be trusted), so this is testing the reliability of solving partial derivatives. Unless you know of a mathematica command that will simultaneously optimize two objective functions that are dependent on each other, the only method I have is solving a series of partial derivatives from those two objective functions. So long story short, I need a reliable way to solve a serires of partial derivatives subject to constraints. Answer
Posted 3 months ago
 Hi Bill, it looked promising at first when the set of solutions was reduced to 2, but it seems odd that there are no solutions when we set ix==24 because MATLAb found a solution where ix==24. I know that MATLAB and Mathematica have different processing approaches, but they shouldn't have completely different answers right? I was having trouble dong this in MATLAB, which is why I was advised to switch over to Mathematica (I'm still pursuing this problem in MATLAB though).The solutions from you code seem to be getting somewhere, though I'm also surprised that I'm not getting symmetrical answers, because the equations are set up in a way that I should expect ax==a1y and a2x==a2y, with ix=/=iy. At least that's what's happened in MATLAB. So my project involves a game theoretic model involving continuous decision nodes. I've denoted the three decision variables as a1, a2 and i, with x and y denoting the two players, so 6 decision variables in total. I'm comparing two cases: a) a strategic scenario where each players have their own utility functions that they are each trying to optimize (but their decision will depend on the other player's choices as well), and b) an optimal scenario that sums up both players' utility functions to optimize total utility. In this model, all of the variables must be nonnegative, and ix and iy must <=24. Here is an example of the two objective functions I want to simultaneously optimize: subject to all the variables being nonnegative and ix and iy being <=24. funx and funy are player X and player Y's individual utility functions, respectively, when there is a strategic scenario. And funt is when I am trying to optimize their combined utility. I've changed the parameters a bit from the original question, since I thought perhaps one of the reasons why some of the solutions were pointing towards the variables approaching 0 is because all the outputs were negative numbers... I don't know, I was just hoping that changing the parameters might help a bit.Thank you again, I really appreciate your help! Answer
Posted 3 months ago
 Have you taken each of the solutions from each system and substituted it into each of the equations on each system and verified that you get 0=0 in all cases?Is this funx=-1/10a1x^2+100Log[a1x+a1y-14]-ix/2+95/100(-1/10a2x^2(1-ix/40)+100Log[a2x+a2y-14]); funy=-1/10a1y^2+100Log[a1x+a1y-14]-iy/2+95/100(-1/10a2y^2(1-ix/40)+100Log[a2x+a2y-14]); NMaximize[{funx+funy,7<=a1x&&7<=a1y&&7<=a2x&&7<=a2y&&0<=ix<=24&&0<=iy<=24},{a1x,a1y,a2x,a2y,ix,iy}] which returns {1254.92,{a1x->26.1329,a1y->26.1329,a2x->39.0282,a2y->39.0282,ix->24.,iy->1.69101*^-7}} a plausible optimum?I had difficulties with variables only being constrained to be non-negative and trying to optimize things like Log[a1x+a1y-14] being negative for small values of the variables and thus the argument to the log being negative and thus the result being complex and NMinimize complaining about this. So as a quick hack I just bumped variables up to 7 or more. If you are expecting the optimum to be symmetric then this might not be a problem, but you should probably think about this.Does this look anything like what you are trying to accomplish?You have ix in the same place in both equations, should one of those be iy? That might account for why the result is not completely symmetric. If I guess the second should be iy and I make that change then both ix and iy converge on 24.Read the documentation on NMaximize, including clicking on the orange Details and Options and see if you can understand why this is written the way it is.If anything in this isn't right then please explain. Answer
Posted 3 months ago
 Hi Bill, I am expecting the results to be symmetrical, at least in this particular example. the ix is in the right spot, by the way. This particular example relates to technology leakage from player x, which is why ix is used in both functions in the second part. Player x invests the maximum amount of 24, but player y has no incentive to invest anything since they can just freeload off player x, hence iy=0.I do have a more math-related question though. When you use Nmaximize, are you basically then just adding the two functions together and maximizing it? I'm not so sure that method would be valid in this case, because that eliminates any strategic interaction since each player isn't acting individually. Adding funx and funy together and then maximizing it is just the equivalent of maximizing funt. BUT!!!!!! The great new is that I actually got the SAME set of solutions on MATLAB!!!! (I've been very hesitant to trust the results I get from MATLAB and Mathematica because things haven't been matching up). So I think the method you used above is great for finding the solutions in the optimal (funx+funy) scenario, but not the strategic scenario. I'm wondering if I should go back to trying to solve the series of partial derivatives?  Answer
Posted 3 months ago
 If your two equations don't have a typographical error then that is great. I was just guessing why the result wasn't completely symmetric.Yes I added the two and maximized. That was because in a previous message I asked you what your optimization model was and didn't seem to get an answer. So I guessed. Exactly what do you mean by "optimize"?thanks Answer
Posted 3 months ago
 Hi Bill, I guess I wasn't exactly clear about what "optimize" means. In a strategic scenario, each player is trying to optimize their own payoffs depending on what the other player is expected to do. Typically, these types of games only involve two variables, one from each player. There are two utility functions e.g. Wx(x,y) and Wy(x,y), and then you take the partial derivative w.r.t to the player for each function (so dWx/dx and dWy/dy) and set them to zero. Then you rearrange both partial derivatives so that you get player x's best response function in terms of y and player y's best response function in terms of x, you plot them together, and that's the nash equilibrium. (sorry if that's too much detail, I'm including this because I thought it might help explain my workings). But usually, there's only two variables, and I have six. However, I looked over my functions again, and I was actually able to use this method to figure out the values for a1x and a1y. So I just need to figure out the values of a2x, a2y, ix, and iy. I also figured out that I needed to tweak the functions a little so that I had enough equations for the unknowns. I also tweaked the parameters since I'm still trying to figure out what values are "reasonable" answers. (sorry, I'm sure it's a little annoying). Here are the new objective functions, I also included the partial derivatives below: So I was able to use Reduce to get a1x in terms of a1y and a1y in terms of a1x. I then plotted them on another program and got a1x=a1y=8.475. I'm also assuming iy=0 because of its first derivative. So I guess the next step is to use Solve on the rest of the nonlinear equations, subject to constraints. However, once again, it seems like I don't have a set of feasible solutions for this second part, since Reduce seems to return False. Are these three equations really not solvable or am I missing something? Thanks again, and please let me know if I can clarify things better. I'm not great at explaining things. Answer
Posted 3 months ago
 Thank you for the explanation. That helps.In your last Reduce you were asking it to solve for iy and there was no iy in the equations. Try Reduce[{95/100(-2/10a2x(1-ix/40)+5/(a2x+a2y-14))==0, 95/100(-2/10a2y(1-ix/40)+5/(a2x+a2y-14))==0, ix^(-1/2)+95/100*1/400a2x^2==0},{a2x,a2y,ix}]//N which returns (a2x == -0.1807-7.9448*I && a2y == -0.1807-7.9448*I && ix == 44.26804+4.0376*I) || (a2x == -0.1807+7.9448*I && a2y == -0.1807+7.9448*I && ix == 44.26804-4.0376*I) If I haven't made a mistake then that is no better because it only finds three complex values.That //N turns the large complicated exact roots into approximate decimal numbers so it is possible to see the form of the results.Since in earlier problems you seemed to have plausible results that were not at the boundaries of your constraints it seems that there may be some problem in the formulation of those three equations now. Answer