ListPlot3D for optimization with varying parameter values

Question

This is a followup question from Alex Trounev's answer to my question in the following post:

3D-Plot optimization results for varying parameter values

Consider my objective function, objF:

1/(12 s^2) (6 k s (s (-1 + t) - 2 t) + (3 (-1 + c) k^4 (-1 + r) (-1 + t))/s^2 - (6 d^4 (-1 + t))/(r^2 s) - (2 k^3 (3 + 3 c (-1 + r) - 3 r + s - 2 r s) (-1 + t))/s + (3 d^3 (k^2 r - 6 k r s + s (2 + 2 c - 4 c q + s + 2 r s)) (-1 + t))/(r^2 s^2) + (2 s^2 (-(1 + 3 r) s (-1 + t) + 3 t + 6 r t))/r + (3 k^2 (-1 + c (-1 + r) (1 + 2 r) (-1 + t) + t + r (-1 - 2 r (1 + s) (-1 + t) + 3 t)))/r + 1/(r s^2) 3 d (2 k s^2 (1 + 2 r (1 + s - t) + c (1 + 2 r) (-1 + t) - t) + k^4 (-1 + r) r (-1 + t) + 2 k^3 r (-1 + c + s - r s) (-1 + t) - 2 s^3 (s + 2 r s + 2 t) + k^2 s (2 + s - r (6 + s + 2 r s) - 2 t + 6 r t + (-1 + r) (1 + 2 r) s t + 2 c (-1 + 3 r + t - 3 r t))) + 1/(r^2 s^2) 3 d^2 (2 k^3 r^2 (-1 + t) + 2 k r s (3 - c (1 + 2 q) + s + 2 r s) (-1 + t) - k^2 r (1 + c - 2 c q - 2 s + 6 r s) (-1 + t) + s^2 (1 + c (-1 + 2 q) (1 + 2 r) (-1 + t) - t + 2 r (1 + s + (-1 + s) t))))

with parameter values: $t=0.2$, $s=2$, $d=0.8$ and $d\leqslant k\leqslant 1$, $\frac{d}{k} \leqslant r \leqslant 1$, $0 \leqslant c \leqslant 1$, and $1 \leqslant q\leqslant 2$.

I'm trying to maximize the above objective function with respect to r and k.

Eventually, I would like to Plot3D each of the optimal values of $objF$, $r$, and $k$ against $c$ and $q$.

For this, I followed Alex Trounev's answer in the above link. The main idea therein is that due to the singularity $\frac{1}{r}$ in the objective function, which generates the infinity with $r \rightarrow 0$, we have to trim $r$, for instance, start with $r=r_0$ where $r_0 = 10^{-3}$.

In all, my mathematica codes are as follows.

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

When run, I get the following result.

Can anyone help?

Answer 1

Here we must add a small value to the upper limit of r and divide the inequality d/k <= r <= 1 into two

Block[{d = .8, s = 2, t = .2, r0 = 10^-3, r1 = 1 + 10^-10}, 
 objF = 1/(12 s^2) (6 k s (s (-1 + t) - 
        2 t) + (3 (-1 + c) k^4 (-1 + r) (-1 + t))/
      s^2 - (6 d^4 (-1 + t))/(r^2 s) - (2 k^3 (3 + 3 c (-1 + r) - 
          3 r + s - 2 r s) (-1 + t))/
      s + (3 d^3 (k^2 r - 6 k r s + 
          s (2 + 2 c - 4 c q + s + 2 r s)) (-1 + 
          t))/(r^2 s^2) + (2 s^2 (-(1 + 3 r) s (-1 + t) + 3 t + 
          6 r t))/r + (3 k^2 (-1 + c (-1 + r) (1 + 2 r) (-1 + t) + t +
           r (-1 - 2 r (1 + s) (-1 + t) + 3 t)))/r + 
     1/(r s^2) 3 d (2 k s^2 (1 + 2 r (1 + s - t) + 
           c (1 + 2 r) (-1 + t) - t) + k^4 (-1 + r) r (-1 + t) + 
        2 k^3 r (-1 + c + s - r s) (-1 + t) - 
        2 s^3 (s + 2 r s + 2 t) + 
        k^2 s (2 + s - r (6 + s + 2 r s) - 2 t + 
           6 r t + (-1 + r) (1 + 2 r) s t + 
           2 c (-1 + 3 r + t - 3 r t))) + 
     1/(r^2 s^2) 3 d^2 (2 k^3 r^2 (-1 + t) + 
        2 k r s (3 - c (1 + 2 q) + s + 2 r s) (-1 + t) - 
        k^2 r (1 + c - 2 c q - 2 s + 6 r s) (-1 + t) + 
        s^2 (1 + c (-1 + 2 q) (1 + 2 r) (-1 + t) - t + 
           2 r (1 + s + (-1 + s) t)))); 
 max = Flatten[
   Table[{c, q, 
     MaxValue[{objF, d <= k <= 1, r <= r1, k*r >= d}, {k, r}]}, {c, 0,
      1, .1}, {q, 1, 2, .1}], 1]; 
 maxk = Flatten[
   Table[{c, q, 
     k /. Last@
       Maximize[{objF, d <= k <= 1, r <= r1, k*r >= d}, {k, r}]}, {c, 
     0, 1, .1}, {q, 1, 2, .1}], 1]; 
 maxr = Flatten[
   Table[{c, q, 
     r /. Last@
       Maximize[{objF, d <= k <= 1, r <= r1, k*r >= d}, {k, r}]}, {c, 
     0, 1, .1}, {q, 1, 2, .1}], 1];] 

In the options ListPlot3D also need to add PlotRange -> {0, 1 + 10^-10}

{ListPlot3D[max, AxesLabel -> {"c", "q", "max"}], 
 ListPlot3D[maxk, PlotRange -> {0, 1}, 
  AxesLabel -> {"c", "q", "maxK"}], 
 ListPlot3D[maxr, AxesLabel -> {"c", "q", "maxR"}, 
  PlotRange -> {0, 1 + 10^-10}]}