# Most efficient way to solve partly integer optimization problem

Let $$0 < u < v < w < 10$$ be natural numbers. I want to maximize $$r$$ across all choices of $$u, v, w$$ subject to the inequalities $$(42 - 100 r) - (6 - 10 r)(w + v) + vw(1 - r), -(42 - 100 r) + (6 - 10 r)(w + u) + uw(1 - r), (42 - 100 r) - (6 - 10 r)(u + v) + vu(1 - r) \ge 0.$$ For example, $$(u,v,w)=(4,5,6)$$ should give $$r \le 4/15$$ (this is the 1st $$(u,v,w)$$ I guessed, but it doesn't have to be the best), and now we need $$83$$ more values. Without the natural numbers condition, I got a working command but the resulting expression was ugly. I added the conditions, but now the command failed:

{
{Maximize[{r, (42 - 100 r) - (6 - 10 r)*(w + v) + v*w *(1 - r) >=
0 && -(42 - 100 r) + (6 - 10 r)*(w + u) + u*w *(1 - r) >=
0 && (42 - 100 r) - (6 - 10 r)*(u + v) + v*u *(1 - r) >= 0,
0 < u < v < w < 10  , u \[Element] PositiveIntegers ,
v \[Element] PositiveIntegers , w \[Element] PositiveIntegers  },
r]},
{\[Placeholder]}
}


Is there any way to do things more efficiently so that the command doesn't fail? If it's possible, I would also like to know how to return the values of $$u, v, w$$ which achieve the maximum.

I just tried

f[u_, v_, w_] :=
Maximize[{r, (42 - 100 r) - (6 - 10 r)*(w + v) + v*w *(1 - r) >=
0 && -(42 - 100 r) + (6 - 10 r)*(w + u) + u*w *(1 - r) >=
0 && (42 - 100 r) - (6 - 10 r)*(u + v) + v*u *(1 - r) >= 0}, r]


so that I could call

Maximize[{f (u, v, w), 0 < u < v < w < 10,
u \[Element] PositiveIntegers, v \[Element] PositiveIntegers,
w \[Element] PositiveIntegers}, {u, v, w}]


next, but when I defined f, I got the output $Failed. Update: I get Maximize[{{4/15, {r -> 4/15}}, True, True, True, True}, {4, 5, 6}]  after defining the function and calling the last code snippet. Certainly, Mathematica can't have read my mind and have known that I already did the case $$(4,5,6)$$ and got $$4/15,$$ so I'm tempted to believe that the first value I tried really did give the maximum. But why does Mathematica output \$Failed for every input then? The function declaration had \$Failed, and this output had \$Failed off to the side.

Final update: In retrospect, what follows would've been easier.

f[u, v, w] = Inverse[{{1, 1, 1}, {u, v, w}, {u^2, v^2, w^2}}]

NMaximize[{r,
VectorGreaterEqual[{f[u, v, w].{1 - r, 6 - 10 r, 42 - 100 r}, {0, 0,
0}}] && 0 < u < v < w < 10 && u \[Element] PositiveIntegers &&
v \[Element] PositiveIntegers &&
w \[Element] PositiveIntegers}, {r, u, v, w},
Method -> "DifferentialEvolution"]


After getting a 1/0 error and "NMaximize: The following constraints are not valid:...", manually substitute in for the vector in VectorGreaterEqual, remove the denominators from the output, and run again.

NMaximize[{r, {(42 - 100 r) (-v + w) + (6 - 10 r) (v^2 - w^2) + (1 -
r) (-v^2 w + v w^2), (42 - 100 r) (u - w) + (6 -
10 r) (-u^2 + w^2) + (1 - r) (u^2 w - u w^2), (42 -
100 r) (-u + v) + (6 - 10 r) (u^2 - v^2) + (1 - r) (-u^2 v +
u v^2)} \[VectorGreaterEqual] {0, 0, 0} &&
0 < u < v < w < 10 && u \[Element] Integers && u > 0 &&
v \[Element] Integers && v > 0 && w \[Element] Integers &&
w > 0}, {r, u, v, w}, Method -> "DifferentialEvolution"]


The usage of NMinimize instead of Minimize gives

NMaximize[{r, (42 - 100 r) - (6 - 10 r)*(w + v) + v*w*(1 - r) >=
0 && -(42 - 100 r) + (6 - 10 r)*(w + u) + u*w*(1 - r) >=
0 && (42 - 100 r) - (6 - 10 r)*(u + v) + v*u*(1 - r) >= 0 &&
0 < u < v < w < 10 && u \[Element] PositiveIntegers &&
v \[Element] PositiveIntegers && w \[Element] PositiveIntegers}, {r, u, v, w},
Method -> "DifferentialEvolution"]
(*{0.321429, {r -> 0.321429, u -> 1, v -> 2, w -> 3}}*)


whereas

N[4/15]
(*0.266667*)


NMaximize[{r, (42 - 100 r) - (6 - 10 r)*(w + v) + v*w*(1 - r) >=
0 && -(42 - 100 r) + (6 - 10 r)*(w + u) + u*w*(1 - r) >=
0 && (42 - 100 r) - (6 - 10 r)*(u + v) + v*u*(1 - r) >= 0 &&
0 < u < v < w < 10 && u \[Element] PositiveIntegers && v \[Element] PositiveIntegers &&
w\[Element]PositiveIntegers},{r, u, v, w},Method-> {"DifferentialEvolution","ScalingFactor"->1}]
(*{0.95122, {r -> 0.95122, u -> 7, v -> 8, w -> 9}}*)


Many thanks from me to @joka for the valuable notice.

• Good point, I just checked it here: wolframalpha.com/input/?i=solve+%2842-100r%29+-+%286-10r%29+*+%283%2B2%29+%2B+2*3+%281-r%29+%3E+0+and+-%2842-100r%29+%2B+%286-10r%29+*+%283%2B1%29+%2B+1*3+%281-r%29+%3E+0+and+%2842-100r%29+-+%286-10r%29+*+%281%2B2%29+%2B+2*1+%281-r%29+%3E+0+for+r I'll have to go back and spot a hidden assumption in my inequalities since I proved earlier that $r \le 3/11$ regardless of $u, v, w.$ – Display name Nov 4 '20 at 7:35
• What really needs to happen is for every entry of the vector here (wolframalpha.com/input/…*+%7B1-r%2C+10-6r%2C+100-42r%7D) to be non-negative, but I figured Mathematica wouldn't be able to eliminate all the fluff and filler, and would take too long to do the computations. – Display name Nov 4 '20 at 7:40
• @Displayname Mathematica evaluates the inequalities ...<r<… for given {u,v,w} without problems. – Ulrich Neumann Nov 4 '20 at 8:16

Get an easy analytical solution with the help of Reduce .

red = Reduce[(42 - 100 r) - (6 - 10 r)*(w + v) + v*w*(1 - r) >=
0 && -(42 - 100 r) + (6 - 10 r)*(w + u) + u*w*(1 - r) >=
0 && (42 - 100 r) - (6 - 10 r)*(u + v) + v*u*(1 - r) >= 0 &&
0 < u < v < w < 10 && u \[Element] Integers &&
v \[Element] Integers && w \[Element] Integers, r];

Maximize[{r, red}, {r, u, v, w}]

(*   {39/41, {r -> 39/41, u -> 7, v -> 8, w -> 9}}   *)


Since $$w,u,v$$ are integers and their common range (rng-see code) is not prohibitively large (1 through 9, inclusive) it's possible to iterate through all the combinations and select the best one.

The code offered below, iterates over the admissible values for the integer variables and records a configuration of variables that respect the constraints and maximize the objective function if the optimum value is strictly greater than any previous recorded value.

(* common range of u, v, w *)
rng = Range

(* formulating the constraints *)
poly1[r_] := 42 - 100 r
poly2[r_] := 6 - 10 r
poly3[r_] := 1 - r

c1[r_,u_,v_,w_] := poly1[r] - poly2[r] (w + v) + v w poly3[r]
c2[r_,u_,v_,w_] := -poly1[r] + poly2[r] (w + u) + u w poly3[r]
c3[r_,u_,v_,w_] := poly1[r] - poly2[r] (v + u) + u v poly3[r]


Preparing the loop

(* ranges for iterators *)
{i0, ie} = Through[{Min, Max}[rng]];
{j0, je} = {k0, ke} = {i0, ie};

(* store results in bag *)
vars = {u, v, w};
args = {0, 0, 0};

(* first entry is not a avalid sol *)
bag = {{0, Thread[vars -> args]}};

(* naive counter *)
c = 0;


The actual loop; Monitor is used to give an update on the current iteration (c).

Monitor[

Do[

ClearAll[c1N, c2N, c3N, cNs, objWCs, res, fMax, xOptRl, args];

(* create constraints for given integer values of u,v,w  *)
c1N[r_] := c1[r, i, j, k] // Evaluate;
c2N[r_] := c2[r, i, j, k] // Evaluate;
c3N[r_] := c3[r, i, j, k] // Evaluate;

(* prep the inequality constraints *)
cNs[r_] := {c1N[r], c2N[r], c3N[r]};

(* list of objective with all constraints *)
objWCs[r_] := Join[{r}, Thread[cNs[r] >= 0], {i < j, j < k}] // Evaluate;

(* the actual optimization under constraints *)
res = Check[NMaximize[objWCs[r], r], $Failed] // Quiet; Which[ (* there was a result *) res // FailureQ /* Not, ( {fMax, xOptRl} = res; If[ (* new opt better than current opt *) fMax > (bag // Last /* First), args = Thread[vars -> {i, j, k}]; (* update res *) bag = Flatten[{bag, {{fMax, args}}}, 1] ] ), (* no opt *) True,$Failed

];

c += 1, {i, i0, ie}, {j, j0, je}, {k, k0, ke}], c]


After evaluating the loop above, the following line

bag // Rest


evaluates to

{
{0.321429, {u -> 1, v -> 2, w -> 3}},
{0.333333, {u -> 5, v -> 7, w -> 8}},
{0.4, {u -> 5, v -> 8, w -> 9}},
{0.5, {u -> 6, v -> 7, w -> 8}},
{0.75, {u -> 6,v -> 8,w -> 9}},
{0.95122, {u -> 7, v -> 8, w -> 9}}
}


Obviously, the optimum value for r under the constraints is 0.95122 and this can be achieved when u -> 7, v -> 8 and w -> 9.

## Edit in order to accommodate a minor consideration raised in the comments

In order to address the issue raised by @user64494 in the comments, it is possible-but not necessary-to replace NMaximize with Maximize inside the loop and obtain the optimum value of r in rational form.

Performing the replacement and evaluating, produces

 {
{9/28, {u -> 1, v -> 2, w -> 3}},
{1/3, {u -> 5, v -> 7, w -> 8}},
{2/5, {u -> 5, v -> 8, w -> 9}},
{1/2, {u -> 6, v -> 7, w -> 8}},
{3/4, {u -> 6, v -> 8, w -> 9}},
{39/41, {u -> 7, v -> 8, w -> 9}}
}


Now, evaluating the following lines of code,

(* make a list of rules for the values of r, u, v, w at the optimum *)
sol = bag // Last /* (MapAt[Rule[r, #] &, #, 1] &) /* Flatten;

(* make a list of the non-negativity constraint functions *)
cs = Map[RightComposition[Apply[Sequence], #] &, {c1, c2, c3}];

(* evaluate the non-negativity constraints at the optimum solution *)
Through[cs[{r, u, v, w} /. sol]] >= 0 // Thread /* Apply[And]


produces True.

This means that all the non-negativity constraints are satisfied at the optimum solution obtained above.

Obviously, the u < v < w constraint is also trivially accommodated at the obtained solution.

Simple inspection of the values of $$u,v,w$$ at the optimum reveals that indeed 7 < 8 < 9 is also True.

• Unfortunately, (42 - 100 r) - (6 - 10 r)*(w + v) + v*w*(1 - r) >= 0 && -(42 - 100 r) + (6 - 10 r)*(w + u) + u*w*(1 - r) >= 0 && (42 - 100 r) - (6 - 10 r)*(u + v) + v*u*(1 - r)>=0 /. {r -> 0.95122, u -> 7, v -> 8, w -> 9} performs False. – user64494 Nov 4 '20 at 12:03
• @user64494 try r->39/41 instead of r->0.95122; the expression you evaluate returns False because 39/41->0.9512195121951219... – joka Nov 4 '20 at 12:25