# NMinimize: How to avoid solutions that do not satisfy constraints within a certain tolerance?

I just started to use Mathematica a few weeks ago. Using NMinimze, I would like to avoid solutions that do not satisfy certain constraints (although they "almost" satisfy them). Do you know how to change the following command to find a solution satisfying "completely" all the constraints, solving the same minimization problem?

NMinimize[{((e*(1 - Sqrt[(g - e)^2 + (f - h)^2]) + (g - e)*(1 -
Sqrt[f^2 + e^2])) + (h*(1 -
Sqrt[(g - e)^2 + (f - h)^2]) + (f - h)*(1 -
Sqrt[g^2 + h^2])))/((g + f)*
Max[1 - Sqrt[(g - e)^2 + (f - h)^2], 1 - Sqrt[g^2 + h^2]]), 0 <= e <= 1, 0 <= f <= 1, e^2 + f^2 == 1, e <= g <= 1, 0 <= h <= f, Sqrt[(g - e)^2 + (f - h)^2] <= 1, g^2 + h^2 <= 1}, {e, f, g, h}]

• You can try to increase PrecisionGoal, but you have to live with the fact that NMinimize is never exact. However the error message that you get indicates that the problem might not have any feasible poins at all. At least, Mathematica cannot find any. See also reference.wolfram.com/language/ref/message/NMinimize/nosat.html – Henrik Schumacher Aug 2 '20 at 16:15
• Apparently, Mathematica has a hard time to find a feasible point at all. If you know one you can try to use FindRoot with this point as starting point. – Henrik Schumacher Aug 2 '20 at 16:27
• If you define your objective= and constraints={...} you can generate exact candidates for the constraints with FindInstance[And @@ constraints, {e, f, g, h}, Reals, 5000] then do MinimalBy[sats, Quiet@Check[objective /. #, Infinity] &] - this gets me e -> 0, f -> 1, g -> 255505/577352, h -> 1/2. This is not optimal of course, but it leads me to believe that varying g and leaving others fixed will find a good minimum. I then did NMinimize[ Prepend[constraints, objective] /. {e -> 0, f -> 1, h -> 1/2}, g] and got {0.535898, {g -> 0.866025}}. – flinty Aug 2 '20 at 16:54
• If you want to 'almost' satisfy your constraints then you could add slack (tolerance) terms to your constraints and try to minimize the total absolute slack (or squared slack) plus your objective. – flinty Aug 2 '20 at 16:57
• @PenelopeBenenati I've provided an answer using slack terms. You basically convert hard constraints like e^2+f^2 == 1 into softer constraints that allow a bit of flexibility like -sef1 < e^2 + f^2 - 1 < sef1 but then you make sure in your objective that the slack is penalized strongly so it can only be a very small tolerance. – flinty Aug 2 '20 at 17:28

This NMinimize::nosat documentation page explains the there are no solutions found if this warning message is displayed.

NMinimize[{((e*(1 - Sqrt[(g - e)^2 + (f - h)^2]) + (g - e)*(1 -
Sqrt[f^2 + e^2])) + (h*(1 -
Sqrt[(g - e)^2 + (f - h)^2]) + (f - h)*(1 -
Sqrt[g^2 + h^2])))/((g + f)*
Max[1 - Sqrt[(g - e)^2 + (f - h)^2], 1 - Sqrt[g^2 + h^2]]),
0 <= e <= 1, 0 <= f <= 1, e^2 + f^2 == 1, e <= g <= 1, 0 <= h <= f,
Sqrt[(g - e)^2 + (f - h)^2] <= 1, g^2 + h^2 <= 1}, {e, f, g, h}]

{0.404445, {e -> 0.00756254, f -> 0.999323, g -> 0.868352,
h -> 0.490688}}


0.00756254^2 + .999323^2 <1

and all the rest is not very valid in a probe too.

Most probable this is ill-conditioned or overly constrained. The is a chance that the 4-dimensionality causes some strong divergence and NMinimize to finish the search for a minimum somewhere near or at the borders, so on the circle on in the circles.

Since {e,f}, and {g,h} are unit circles around zero in 2 dimensions there is a chance to use the constraints for visual control of the solution.

objective = {((e*(1 - Sqrt[(g - e)^2 + (f - h)^2]) + (g - e)*(1 -
Sqrt[f^2 + e^2])) + (h*(1 -
Sqrt[(g - e)^2 + (f - h)^2]) + (f - h)*(1 -
Sqrt[g^2 + h^2])))/((g + f)*
Max[1 - Sqrt[(g - e)^2 + (f - h)^2], 1 - Sqrt[g^2 + h^2]])}

constraints = {0 <= e <= 1, 0 <= f <= 1, e^2 + f^2 == 1, e <= g <= 1,
0 <= h <= f, Sqrt[(g - e)^2 + (f - h)^2] <= 1, g^2 + h^2 <= 1}

Solve[constraints, {e, f, g, h}]

{{f -> ConditionalExpression[Sqrt[
1 - e^2], (0 < g <= Sqrt[1 - h^2] && 0 <= e <= g &&
1/2 < h <= 1) || (0 < g <= Sqrt[2 h - h^2] && 0 <= e <= g &&
0 <= h <= 1/2) || (0 <= h <= 1/2 &&
Sqrt[2 h - h^2] < g <= Sqrt[1 - h^2] &&
g/2 - 1/2 Sqrt[(4 h^2 - g^2 h^2 - h^4)/(g^2 + h^2)] <= e <=
g)]}, {e -> ConditionalExpression[0, 0 <= h <= 1],
f -> ConditionalExpression[1, 0 <= h <= 1],
g -> ConditionalExpression[0, 0 <= h <= 1]}}


The solutions on the border of the unit circle are irrelevant.

This can be used in a

RegionPlot3D[(0 < g <= Sqrt[1 - h^2] && 0 <= e <= g &&
1/2 < h <= 1) || (0 < g <= Sqrt[2 h - h^2] && 0 <= e <= g &&
0 <= h <= 1/2) || (0 <= h <= 1/2 &&
Sqrt[2 h - h^2] < g <= Sqrt[1 - h^2] &&
g/2 - 1/2 Sqrt[(4 h^2 - g^2 h^2 - h^4)/(g^2 + h^2)] <= e <=
g), {e, 0, 1}, {g, 0, 1}, {h, 0, 1}]


NMinimize[{((e*(1 - Sqrt[(g - e)^2 + (f - h)^2]) + (g - e)*(1 -
Sqrt[f^2 + e^2])) + (h*(1 -
Sqrt[(g - e)^2 + (f - h)^2]) + (f - h)*(1 -
Sqrt[g^2 + h^2])))/((g + f)*
Max[1 - Sqrt[(g - e)^2 + (f - h)^2], 1 - Sqrt[g^2 + h^2]]),
f == Sqrt[
1 - e^2], (0 < g <= Sqrt[1 - h^2] && 0 <= e <= g &&
1/2 < h <= 1) || (0 < g <= Sqrt[2 h - h^2] && 0 <= e <= g &&
0 <= h <= 1/2) || (0 <= h <= 1/2 &&
Sqrt[2 h - h^2] < g <= Sqrt[1 - h^2] &&
g/2 - 1/2 Sqrt[(4 h^2 - g^2 h^2 - h^4)/(g^2 + h^2)] <= e <=
g)}, {e, g, h}]


{-1.38043*10^13, {e -> 9.30374*10^-6, g -> 0.887035, h -> 0.53828}}


This again is very marginal. Might be the complete plane for e near zero is degenerate a solution.

That the problem is ill-posed might be due to an error in the question.

Or it is an situation like this

Plot[{1/(1 + g), Sqrt[1/4 + g^2], 1/4 + g^2}, {g, 0, 1}]


• Mathematica v12 (windows 64) doesn't show an error NMinimize::nosat – Ulrich Neumann Aug 3 '20 at 6:44

As indicated by Henrik Schumacher the constraints are fullfilled numerically

constraint = {0 <= e <= 1, 0 <= f <= 1, e^2 + f^2 == 1, e <= g <= 1,0 <= h <= f, Sqrt[(g - e)^2 + (f - h)^2] <= 1, g^2 + h^2 <= 1};
mini = NMinimize[{((e*(1 - Sqrt[(g - e)^2 + (f - h)^2]) + (g - e)*(1 -Sqrt[f^2 + e^2])) + (h*(1 -Sqrt[(g - e)^2 + (f - h)^2]) + (f - h)*(1 -Sqrt[g^2 + h^2])))/((g + f)*Max[1 - Sqrt[(g - e)^2 + (f - h)^2], 1 - Sqrt[g^2 + h^2]])
,constraint}, {e, f, g, h}]


The constraints #3 and #6 seem to be violated constraint /. mini[[2]] ({True, True, False, True, True, False, True})

Further inspection shows

constraint[[3]] /. Equal -> Subtract /. mini[[2]]
(*5.33813*10^-11*)

constraint[[6]] /. LessEqual -> Subtract /. mini[[2]]
(*1.66635*10^-9*)


that both constraints are fullfilled numerically quite well!

Here's how you can do it by adding some slack into the constraints and punishing slack in the objective:

SeedRandom[1];

(* the function you're trying to minimize *)
objective = ((e*(1 - Sqrt[(g - e)^2 + (f - h)^2]) + (g - e)*(1 -
Sqrt[f^2 + e^2])) + (h*(1 -
Sqrt[(g - e)^2 + (f - h)^2]) + (f - h)*(1 -
Sqrt[g^2 + h^2])))/((g + f)*
Max[1 - Sqrt[(g - e)^2 + (f - h)^2], 1 - Sqrt[g^2 + h^2]]);

(* these are the hard constraints *)
constraints = {
0 <= e <= 1,
0 <= f <= 1,
e^2 + f^2 == 1,
e <= g <= 1,
0 <= h <= f,
Sqrt[(g - e)^2 + (f - h)^2] <= 1,
g^2 + h^2 <= 1
};

(* these constraints are softer and allow for a bit of slack *)
slackedConstraints = {
0 - se <= e <= 1 + se,
0 - sf <= f <= 1 - sf,
-sef1 < e^2 + f^2 - 1 < sef1,
e <= g <= 1,
0 - sh <= h <= f + sh,
Sqrt[(g - e)^2 + (f - h)^2] <= 1,
g^2 + h^2 - 1 <= 0
};
variables = {e, f, g, h};
slackterms = {se, sf, sh, sef1};

(* solve it and harshly punish too much total squared slack *)
sol = Last[
NMinimize[{objective + 10^10*Total[slackterms^2],
slackedConstraints}, Join[variables, slackterms]]]

(* RESULT:
{e -> 0.25283, f -> 0.967511, g -> 0.944242, h -> 0.329154,
se -> 4.51664*10^-14, sf -> -2.52757*10^-13, sh -> 3.93093*10^-14,
sef1 -> 1.92914*10^-7} *)

objective /. sol
(* result: 0.304607 *)

(* Substitute back into the hard constraints to check if any violated *)
constraints /. sol
(* {True, True, False, True, True, True, True} *)

(* hard constraint #3 is violated, but only by a tiny amount: *)
e^2 + f^2 /. sol
(* result 1. *)

• Thank you a lot @flinty! This is very helpful for me! By the way, in this problem I understood that in minimizing the expression, we have something similar to a $0/0$, which in this specific case, if all the contraints are satisfied, should be equal exactly to $1$. I understood it by inspecting now more carefully the problem (by hand), but it seems difficult to see this fact with Mathematica. Thank you once again! – Penelope Benenati Aug 2 '20 at 17:31
• Do you think could be useful to combine your approach with the use of limits for the slack variables (which should be set to go to $0$)? – Penelope Benenati Aug 2 '20 at 17:36
• If you know particular variables are going to zero, then you should change your objective and remove the redundant constraints. For example objective /. {e -> 0, g -> 0} gives ((1 - Sqrt[(f - h)^2]) h + (f - h) (1 - Sqrt[h^2]))/(f Max[1 - Sqrt[(f - h)^2], 1 - Sqrt[h^2]]) – flinty Aug 2 '20 at 17:43
• Thank you @flinty! – Penelope Benenati Aug 2 '20 at 17:49