NMinimize : different results when changing a little bit the constraints

Question

I want to solve multiple equations so I thought using NMinimize. The basic idea is that if you have multiple equations with multiple variables :

$$\left\{\begin{split}&f(x,y,z)=0\\ &g(x,y,z)=0 \\& h(x,y,z)=0\\\end{split}\right.$$

You can find the solution by minimizing : $|f(x,y,z)|+|g(x,y,z)|+|h(x,y,z)|$.

So this is what I'm trying to do here in this code with the functions p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z],Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z], Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z].

My problem is that when I'm changing a little bit the condition : z>... the result becomes different.

Here is the introduction to the NMinimize :

a0=0.0748294;
b0=0.629316;
h[x_, y_] = 
  a*x*Log[x] + b*(1 - x - y)*Log[1 - x - y] + c*y*Log[y] - x^2 - d*x*y;
"------------------------------"
Hx[x_, y_] = D[h[x, y], x];
Hy[x_, y_] = D[h[x, y], x];
hxx[x_, y_] = D[h[x, y], x, x];
hxy[x_, y_] = D[h[x, y], x, y];
hyy[x_, y_] = D[h[x, y], y, y];
det[x_, y_] = Det[{{hxx[x, y], hxy[x, y]}, {hxy[x, y], hyy[x, y]}}];
p[x_, y_] = h[x, y] - x*Hx[x, y] - y*Hy[x, y];
h[a_, b_, c_, d_, x_, y_] = h[x, y];
p[a_, b_, c_, d_, x_, y_] = p[x, y];
Hx[a_, b_, c_, d_, x_, y_] = Hx[x, y];
Hy[a_, b_, c_, d_, x_, y_] = Hy[x, y];
det[a_, b_, c_, d_, x_, y_] = det[x, y];
x = 0.001;
y = 0.001;
w = 0.4;

And here are the several NMinimize changing a little bit the condition. As you can see, the condition I'm changing is z>0.2,z>0.3,z>0.4. And as you can see all the results give $z>0.4$ so none of them contradicts any of the conditions.

sol = NMinimize[{Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + 
    Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + 
    Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]], 
   c > 0 && d > 0 && z > 0.2 && z < 1 - w && 
    det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
  AccuracyGoal -> 20, PrecisionGoal -> 22, MaxIterations -> 1000]

sol = NMinimize[{Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + 
    Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + 
    Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]], 
   c > 0 && d > 0 && z > 0.3 && z < 1 - w && 
    det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
  AccuracyGoal -> 20, PrecisionGoal -> 22, MaxIterations -> 1000]

sol = NMinimize[{Abs[p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z]] + 
    Abs[Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z]] + 
    Abs[Hy[a0, b0, c, d, x, y] - Hy[a0, b0, c, d, w, z]], 
   c > 0 && d > 0 && z > 0.4 && z < 1 - w && 
    det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
  AccuracyGoal -> 20, PrecisionGoal -> 22, MaxIterations -> 1000]

The results :

-> {0., {c -> 0.584662, d -> 2.09663, z -> 0.486347}}

-> {1.94289*10^-16, {c -> 0.560407, d -> 2.26232, z -> 0.507207}}

-> {1.8735*10^-16, {c -> 0.114809, d -> 3.31333, z -> 0.571615}}

As you can see they are all different. How can I do to find the exact result ?

EDIT :

What I ended up doing is to write :

sol = Table[
  NMinimize[{((p[a0, b0, c, d, x, y] - 
          p[a0, b0, c, d, w, z])^2 + (Hx[a0, b0, c, d, x, y] - 
          Hx[a0, b0, c, d, w, z])^2 + (Hy[a0, b0, c, d, x, y] - 
          Hy[a0, b0, c, d, w, z])^2)*1000000, 
    c > 0 && d > 0 && z > nn && z < 1 - w && 
     det[a0, b0, c, d, w, z] > 0}, {c, d, z}, Reals, 
   AccuracyGoal -> 20, PrecisionGoal -> 22, 
   MaxIterations -> 1000], {nn, 0.3, 0.575, 0.025}]

And I plotted the results, as one can see it converges :

ListPlot[sol[[All, 1]]]
ListPlot[Table[z /. sol[[n, 2, 3]], {n, 1, Length[sol]}]]
ListPlot[Table[d /. sol[[n, 2, 2]], {n, 1, Length[sol]}]]
ListPlot[Table[c /. sol[[n, 2, 1]], {n, 1, Length[sol]}]]

I'm not sure it's a good method, even now convergence is generally a good sign.

Answer 1

Edit Corrected error, not taking second equation into account.

In this case solution with NMinimize is not good, since you get a whole curve of solutions, not only single points.

If i got it right the difference of the p and the Hx have to be zero. Hy is the same as Hx. I rationalized to get a solution with Solve.

sol = Solve[{p[a0, b0, c, d, x, y] - p[a0, b0, c, d, w, z] == 0, 
Hx[a0, b0, c, d, x, y] - Hx[a0, b0, c, d, w, z] == 0} // 
Rationalize[#, 0] &, {c, d}]

(csol[z_] = c /. First@sol) // N

(*   (7.99877*10^-29 (-1.16357*10^37 - 3.27364*10^37 z - 
 1.91963*10^37 Log[0.6\[VeryThinSpace]- 1. z] + 
 1.92156*10^37 z Log[0.6\[VeryThinSpace]- 1. z]))/(1.68711*10^7 + 
 2.44235*10^9 z Log[z])   *)

(dsol[z_] = d /. First@sol) // N

(*   -((9.51773*10^-11 (3.68704*10^12 + 
6.61204*10^12 Log[0.6\[VeryThinSpace]- 1. z]))/(-1. + 1000. z))   *)

.

ParametricPlot3D[{csol[z], dsol[z], z}, {z, 1/5, 3/5}, 
 AspectRatio -> 1, AxesLabel -> {c, d, z}, PlotStyle ->  Thick, 
 RegionFunction -> 
 Function[{c, d, 
z}, -d^2 - 
  157968671/(138459256 (3/5 - z)) + (157329 d)/(125000 (3/5 - 
       z)) - (3625853 c)/(2000000 z) + (157329 c)/(250000 (3/5 - 
       z) z) > 0 && c > 0]]

Test with the point minima you got

dsol[0.507207]
(*   2.26232   *)

csol[0.507207]
(*   0.560406   *)