Need help maximizing an expression

Question

I'm trying to use mathematica to maximize an equation subject to constraints, and it isn't working. The Maximize commmand just runs for hours without returning anything, and NMaximize does return solutions, though they aren't maxima, and giving it different neighbourhoods of the arguments gives wildly different solutions. I'm unsure if it is the equation itself, my code, or my approach, or what, so I'd appreciate any advice.

$\max\limits_{\{s,x,t,y\}} \sqrt{s} + \sqrt{t} -xs-ty-\dfrac{50st}{1+x+y} \qquad \text{st.} \qquad \{t,s\} \in [0,10^6]^2\ \land\ x,y\geq0 $

NMaximize[{Sqrt[s] + Sqrt[t] - x s - y t - 50 s t/(1 + x + y), 1000000 >= s >= 0, 1000000 >= t >= 0, x >= 0, y >= 0}, {s, t, x, y}]

returned

{19.6578, {s -> 386.427, t -> 1.27273*10^-8, x -> 2.92059*10^-17, y -> 4427.66}}

That's obviously local, at best, since setting $s=10^6; x=y=t=0$ gives a functional value of $10^3$. I ran the following, hoping that by limiting the algorithm I might get different maxima:

Do[{Print[NMaximize[{Sqrt[s] + Sqrt[t] - x s - y t - 50 s t/(1 + x + y), 1000000 >= s >= 0, x >= 0, y >= 0, 0 <= t <= 1000000}, {{s, i 1000 - 1, i 1000}, t, x, y}]]}, {i, 1000}]

Running that loop gave me hundreds of maxima, with a peak functional value of 1000. They all seem to be similar, for example {1000.,{s->1.*10^6,t->1.25033*10^-9,x->1.89242*10^-17,y->7070.2}}, with values of $t,x$ being approximately zero.

What is happening? Is mathematica rounding its computations, and is that why it returns so many different maxima? I get lots of messages like NMaximize::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.. Why doesn't the Maximize command return anything? Would increasing working precision help? Am I approaching this wrong? Theoretically there must exist a maximum over the domain I've specified, I've even limited $x,y$ to the same interval to no avail. Is there any other way to go about this?

Answer 1

So there's a lot going on here. There are many ways to get at one of the (uncountably infinitely many) solutions to your particular problem. But you asked what's going on, so here's a few tips:

What can make Maximize and NMaximize grumpy?

1. Generically: NMaximize always executes at finite precision (unlike Maximize which can work with infinite precision). So numerical error can cause it to run out of steam before finding a true maximum. Such failure is more likely when your expression is more complicated. For any of Mathematica's numeric methods, you generally have a better chance on landing on an actual solution if you set WorkingPrecision and MaxIterations higher. Of course you'll pay a price in speed / memory etc.

2. Generically: Worth considering, as KraZug suggested, comparison between multiple Method-> and seeing who resolves with minimal strife. There's a lot to learn about different methods ( http://reference.wolfram.com/language/tutorial/ConstrainedOptimizationGlobalNumerical.html ) $\leftarrow$ worth paying attention to if this is something you're going to be doing a lot of.

3. Specific to this problem: You're making Mathematica do more work than it needs to. The function Sqrt is an unnecessary headache (formally multiple-valued) that tends to complicate both numeric and analytic solutions. I assume you mean to make an unambiguous choice of Sqrt, like Abs[Sqrt[#]]&. We can resolve by either making our expression more complicated, by explicitly introducing Abs[Sqrt[_]], or we can simplify our lives by introducing a friendly change of variables {a^2->s, b^2->t}, yielding: $$a+b-a^2 x-b^2 y-50\frac{ a^2 b^2}{1+x+y}$$

Now maximized subject to constraints: $$\begin{align} 0&\le a \le 10^3\,\\ 0&\le b \le 10^3\,\\ 0&\le x \,\\ 0&\le y \end{align} $$

You can verify with the StepMonitor argument in Minimize that this converges much faster for the same WorkingPrecision.

    origExpr = Sqrt[s] + Sqrt[t] - x s - y t - 50 s t/(1 + x + y);
    simpleExpr = a + b - a^2 x - b^2 y - 50 a^2 b^2/(1 + x + y);
    (* Original with Sqrts *)
    {ansOrig, stepOrig} = NMaximize[{origExpr, 
            10^6 >= s >= 0, 10^6 >= t >= 0, x >= 0, y >= 0},
            {s, t, x, y}, StepMonitor :> Sow[{s, t, x, y, origExpr}], 
            WorkingPrecision -> 500, MaxIterations -> 400] // Reap;
    (* Simplified w/ change of variables *)        
    {ansSimple, stepSimple} = NMaximize[{simpleExpr,
            10^3 >= a >= 0, 10^3 >= b >= 0, 0 <= x,0 <= y}, 
            {a, b, x, y}, StepMonitor :> Sow[{a, b, x, y, simpleExpr}],
            WorkingPrecision -> 500, MaxIterations -> 400] // Reap;
    ListPlot[{(Last /@ Flatten[stepOrig, 1]), 
          (Last /@ Flatten[stepSimple, 1])}, 
          PlotRange -> {-10^3, 10^3}, Frame -> True, 
          PlotLegends -> {origExpr // Style, simpleExpr // Style}, 
          Joined -> True, FrameLabel -> {"steps", "Expr Val Tested"}]

Yields:

4. Specific to this problem: You can see by inspection that you're formally dealing with an uncountably infinite set of solutions, some of whose choices render other parameters completely irrelevant. Consider an extreme: {a->0}. Now we're looking at maximizing: $$b- b^2y\,.$$ This is independent of $x$, which can take on any value, and is maximized when $b$ grows and $y$ vanishes. This is symmetric about {a->b,b->a,x->y,y->x} so we have the uncountably infinite set of solutions: {a->10^3,b->0,x->0,y>=0}||{a->0,b->10^3,x>=0,y->0}. In general for solution hunting, it can be good to break degeneracy and restore symmetry later.

5. Specific to this problem: It's not hard to see that any analytic root finding in an analytic Maximize is going to be a mess for this guy. Consider what's involved in asking for something like:

   Reduce[{simpleExpr==10^3, 10^3 >= a >= 0, 10^3 >= b >= 0, 0 <= x,0 <= y},{a,b,x,y},Reals]
   

Then tracing through the result.

Answer 2

Changing the method on the NMaximize finds the solution with no problems:

NMaximize[{Sqrt[s] + Sqrt[t] - x s - y t - 50 s t/(1 + x + y), 
  1000000 >= s >= 0, 1000000 >= t >= 0, x >= 0, y >= 0}, {s, t, x, y},
  Method -> "SimulatedAnnealing"]

(* {1000., {s -> 0., t -> 1.*10^6, x -> 0.131227, y -> 0.}} *)

I don't know why the default method has a problem. The example in the question puts $t$ to (almost) zero but then doesn't manage to follow that by increasing $s$.

Answer 3

The FindMaximum command does the job:

FindMaximum[{Sqrt[s] + Sqrt[t] - x s - y t - 50 s t/(1 + x + y),  

1000000 >= s >= 0 && 1000000 >= t >= 0 && x >= 0 && y >= 0}, {s, t,  x, y}]

{1000.,{s->1.*10^6,t->2.90426*10^-11,x->0.,y->92422.5}}