# FindMaximum returns a solution that violates the constraints

I am looking at the following optimization problem

\begin{align*} \max\ & 1000 r_1 + \frac{1}{2}r_2 + \frac{1}{3}r_3\\ \text{s.t. }& 1000^2 r_1 + \frac{1}{4}r_2 + \frac{1}{9}r_3 = \frac{1}{9},\\ & 1000^2 p_1 + \frac{1}{4}p_2 + \frac{1}{9}p_3 = \frac{1}{9},\\ & r_1 + r_2 + r_3 + r_4 = 2\\ & p_1 + p_2 + p_3 + p_4 = 2\\ & 0\leq p_i\leq 1,\quad i = 1,\dots,4\\ & r_i^2\leq p_i,\quad i = 1,\dots,4 \end{align*} with the following Mathematica code (it is clear that $$r_1=p_1=0$$, $$r_2=p_2=0$$, $$r_3=p_3=1$$, $$r_4=p_4=1$$ is a feasible solution so I picked it to be the initial solution)

FindMaximum[{1000 r1 + 1/2 r2 + 1/3 r3,
{1000^2 r1 + 1/4 r2 + 1/9 r3 == 1/9,
r1 + r2 + r3 + r4 == 2,
r1^2 <= p1, r2^2 <= p2, r3^2 <= p3, r4^2 <= p4,
1000^2 p1 + 1/4 p2 + 1/9 p3 == 1/9,
p1 + p2 + p3 + p4 == 2,
0 <= p1 <= 1, 0 <= p2 <= 1, 0 <= p3 <= 1, 0 <= p4 <= 1}
},
{{r1, 0}, {r2, 0}, {r3, 1}, {r4, 1}, {p1, 0}, {p2, 0}, {p3, 1}, {p4, 1}}, AccuracyGoal -> 10]


{0.45521, {r1 -> -6.30073*10^-8, r2 -> 0.268328, r3 -> 0.963328,
r4 -> 0.768345, p1 -> 0., p2 -> 0.0719999, p3 -> 0.928, p4 -> 1.}}


If you look at the constraint that $$1000^2 p_1 + \frac{1}{4}p_2 + \frac{1}{9}p_3 = \frac{1}{9}$$, you can see that in the solution Mathematica returns, the left-hand side evaluates to $$0.12111\cdots$$, so the constraint is not satisfied.

Is this a bug or just a numerical inaccuracy? Just the value seems too large compared with $$1/9$$ -- the additive error is not at the scale of $$10^{-4}$$ or $$10^{-5}$$.

• Yes, I can reproduce that, it fails these conditions {r1^2 <= p1, 1000000 p1 + p2/4 + p3/9 == 1/9, p1 + p2 + p3 + p4 == 2}, with the largest error on the one you mentioned. Aug 21, 2019 at 11:23
• I think r1^2 <= p1 and p1 + p2 + p3 + p4 == 2 are pretty minor violations, as the additive error is around $10^{-4}$ or $10^{-5}$, which is understandably to be usual precision problem. But the constraint I pointed out is too big to be called a precision problem. Aug 21, 2019 at 12:07

It looks like a bug. If you define the error variable like this

sol = FindMaximum[{
1000 r1 + 1/2 r2 + 1/3 r3,
{
1000^2 r1 + 1/4 r2 + 1/9 r3 == 1/9,
r1 + r2 + r3 + r4 == 2,
r1^2 <= p1,
r2^2 <= p2,
r3^2 <= p3,
r4^2 <= p4,
1000^2 p1 + 1/4 p2 + 1/9 p3 == 1/9,
error == 1000^2 p1 + 1/4 p2 + 1/9 p3 - 1/9, (* just change this *)
p1 + p2 + p3 + p4 == 2,
0 <= p1 <= 1, 0 <= p2 <= 1, 0 <= p3 <= 1, 0 <= p4 <= 1}
}, {{r1, 0}, {r2, 0}, {r3, 1}, {r4, 1}, {p1, 0}, {p2, 0}, {p3,
1}, {p4, 1}, {error, 0}}, AccuracyGoal -> 10]


you get the same solution:

{0.45521, {r1 -> -6.30073*10^-8, r2 -> 0.268328, r3 -> 0.963328,
r4 -> 0.768345, p1 -> 0., p2 -> 0.0719999, p3 -> 0.928, p4 -> 1.,
error -> 5.75352*10^-18}}


however

1000^2 p1 + 1/4 p2 + 1/9 p3 - 1/9 /. Last[sol]
error /. Last[sol]


yields

0.00999999
5.75352*10^-18


although clearly the two figures should be identical. I suggest you contact Wolfram.

• Thanks! Contacted Wolfram. Confirmed as a bug. Aug 21, 2019 at 23:59

I don't find it a bug. The feasible set consists of only one element

Reduce[{1000^2 r1 + 1/4 r2 + 1/9 r3 == 1/9, r1 + r2 + r3 + r4 == 2,
r1^2 <= p1, r2^2 <= p2, r3^2 <= p3, r4^2 <= p4, 1000^2 p1 + 1/4 p2 + 1/9 p3 == 1/9,
p1 + p2 + p3 + p4 == 2, 0 <= p1 <= 1, 0 <= p2 <= 1, 0 <= p3 <= 1, 0 <= p4 <= 1},
{r1, r2,  r3, r4, p1, p2, p3, p4}, Reals]


r1 == 0 && r2 == 0 && r3 == 1 && r4 == 1 && p1 == 0 && p2 == 0 && p3 == 1 && p4 == 2 - p3

and Mathematica does her best.