3D Plot maximization results with varying parameter values

I'm trying to maximize the following objective function with respect to $$r$$ and $$k$$

$$objF = \frac{d^3+3 d^2 (k (r-1)-s)+3 d \left(k^2 ((r-3) r+1)+2 k r s-2 (r-1) s^2\right)+k^3 (r ((r-1) r+3)-1)-3 k^2 r^2 s+3 (r-1) s^3}{6 (r-1) s^2}$$

under the conditions of $$0\leq d\leq 1$$, $$s\geq 2 d$$, $$d\leq k\leq s$$, and $$0\leq r<\frac{d}{k}$$. And plot the maximum of the objective function, $$r^*$$, and $$k^*$$ against the varying parameter values of $$d \in [0,1]$$ and $$s\in [0,2]$$.

My Mathematica code is:

Block[{t = 0}, objF = 1/( 6 (-1 + r) s^2) (d^3 + k^3 (-1 + r (3 + (-1 + r) r)) + 3 d^2 (k (-1 + r) - s) - 3 k^2 r^2 s + 3 (-1 + r) s^3 + 3 d (k^2 (1 + (-3 + r) r) + 2 k r s - 2 (-1 + r) s^2)); max = Flatten[Table[{d, s, MaxValue[{objF, 0 <= d <= 1, s >= 2 d, d <= k <= s, 0 <= r < d/k}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxk = Flatten[Table[{d, s, k /. Last@Maximize[{objF, 0 <= d <= 1, s >= 2 d, d <= k <= s, 0 <= r < d/k}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxr = Flatten[Table[{d, s, r /. Last@Maximize[{objF, 0 <= d <= 1, s >= 2 d, d <= k <= s, 0 <= r < d/k}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1];] {ListPlot3D[max, AxesLabel -> {"d", "s", "objF"}], ListPlot3D[maxk, PlotRange -> {0, 2}, AxesLabel -> {"d", "s", "k"}], ListPlot3D[maxr, PlotRange -> {0, 1}, AxesLabel -> {"d", "s", "r"}]}

And I got the following result. I have two issues (or questions) with this result.

First, I'm curious what the gray areas in the last diagram mean.

Second, in order to check the concrete value of the objective function with specific parameter values, e.g., $$s=1.8$$ and $$d=0.7$$, I used the following code:

Block[{s = 1.8, d = 0.7}, objF = 1/( 6 (-1 + r) s^2) (d^3 + k^3 (-1 + r (3 + (-1 + r) r)) + 3 d^2 (k (-1 + r) - s) - 3 k^2 r^2 s + 3 (-1 + r) s^3 + 3 d (k^2 (1 + (-3 + r) r) + 2 k r s - 2 (-1 + r) s^2)); max = MaxValue[{objF, d <= k <= s, 0 <= r < d/k}, {k, r}]; maxk = k /. Last@Maximize[{objF, d <= k <= s, 0 <= r < d/k}, {k, r}]; maxr = r /. Last@Maximize[{objF, d <= k <= s, 0 <= r < d/k}, {k, r}]] {max, maxk, maxr}

And the result is: That is, the maximum $$objF=0$$, $$k^* =0$$, and $$r^* =0$$. This result is very strange since a quick inspection of the above diagrams show that at $$s=1.8$$ and $$d=0.7$$, we have both the maximum $$objF$$ and $$k^*$$ positive.

• How is this different from mathematica.stackexchange.com/a/201660/1783 and what is wrong with the answer there? Your code here does not maximize the objective function. – bill s Jul 7 at 20:00
• Since the analytical solution turns out to be very complicated with a lot of conditional expression with Root[#], I am trying a numerical solution by simulation. Thanks for your help, Bill! – ppp Jul 7 at 20:08
• Hi Bill, shall I merge the two posts into one or delete the previous post? I would really appreciate if anyone could help me on this post. Thanks! – ppp Jul 8 at 0:56
• My suggestion would be to use the analytical solution and plot that for specific values (rather than trying to redo the maximization). You'll need to learn to use root objects (see help for Root). They're not really that hard! – bill s Jul 8 at 1:02
• Thanks for your comments, bill s. As you already know, the analytical solution turns out to be very complicated (its expression covers more than one entire page). Do you think it is still something manageable? If I can get an analytical solution as you suggest, then there is actually no need to do the numerical simulation. – ppp Jul 8 at 2:42

The answer of your first question is quite easy: You try to plot a value which is near zero and cut the negativ part (PlotRange).

I'm wondering why you didn't use NMaximize? Here my approach(firstly without restriction 0<d<1, d <= s/2)

max[d_?NumericQ, s_?NumericQ] := {d, s,NMaximize[{1/(6 (-1 + r) s^2) (d^3 +k^3 (-1 + r (3 + (-1 + r) r)) + 3 d^2 (k (-1 + r) - s) - 3 k^2 r^2 s + 3 (-1 + r) s^3 + 3 d (k^2 (1+ (-3 + r) r) + 2 k r s - 2 (-1 + r) s^2)),
d <= k <= s, 0 <= r k < d}, {k, r}]}

That's it.

Plot the results (requires some time...)

ListPlot3D[DeleteCases[Table[{d, s, r /. max[d, s][[3, 2]]}, {d, 0.05, 1, .05}, {s, 2 d, 1, .05}], {}] // Flatten[#, 1] &, PlotRange -> All, AxesLabel -> {d, s, r}] ListPlot3D[DeleteCases[Table[{d, s, k /. max[d, s][[3, 2]]}, {d, 0.05, 1, .05}, {s, 2 d, 1, .05}], {}] // Flatten[#, 1] &, AxesLabel -> {d, s, k}] ListPlot3D[DeleteCases[Table[{d, s, max[d, s][[3, 1]]}, {d, 0.05, 1, .05}, {s, 2 d, 1, .05}], {}] // Flatten[#, 1] &, AxesLabel -> {d, s, max}] • Thanks, Ulrich! May I ask you a question? I tried the simulation without the condition $s > 2d$. And also, tried to do the Plot3D with the range of both $d$ and $s$ starting from 0. Accordingly, the parameter value range in the code would be: {d, 0, 1, .05}, {s, 0, 2, .05}. When run, I get a list of error messages; one of them is "Infinite expression 1/0. encountered". Do you know why? Can you please help me once more? I really appreciate! – ppp Jul 8 at 18:36
• As for the range of parameter values in Plot3D, I also tried many different variations. For example, {d, 0.01, 1, .05}, {s, 0.01, 2, .05} or {d, 0.05, 1, .05}, {s, d, 2, .05}. All these generate the similar error messages. – ppp Jul 8 at 18:40
• With d=0 I also got problems, that's why I started with small d. Plot3D simulated very long time. – Ulrich Neumann Jul 9 at 7:54
• Hi Ulich, may I ask you one more question? When I use PlotRange -> {0,1} for $r$ in the first diagram, I get basically the same result, but the area is in gray, not in orange. Do you know why? – ppp Jul 9 at 17:26
• Try a small lower limit in the plotrange. – Ulrich Neumann Jul 9 at 20:27