1
$\begingroup$

My objective function, $f(r,k;s,d)$, takes three different expressions, denoted by $f^1$, $f^2$, and $f^3$, respecitvely, depending on the conditions for the arguments and parameters as follows:

$f^1(r,k;s,d) = \frac{3 s \left(d^2-k^2 r\right)+k^2 r (3 d+k r)-6 d s^2+3 s^3}{6 s^2}$ when $0\leq d\leq 1,s\geq 2 d,0\leq k<d,0\leq r\leq 1$,

$f^2(r,k;s,d) = \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}$ when $0\leq d\leq 1,s\geq 2 d,d\leq k\leq s,0\leq r<\frac{d}{k}$,

$f^3(r,k;s,d) = \frac{d^3+3 d^2 k r+3 d r \left(k^2 (r-1)-2 s^2\right)+k^3 (r-1)^2 r+3 r s^3}{6 r s^2}$ when $0\leq d\leq 1,s\geq 2 d,d\leq k\leq s,\frac{d}{k}\leq r\leq 1$.

My ultimate goal is to find the optimal value of $r$ and $k$ at which the objective function $f$ gets maximized over all the three cases. It has turned out that working with analytical solutions is quite daunting. So I'm working with numerical simulations limited to the parameter values of $d \in [0,1]$ and $s \in [0,2]$.

My strategy is as follows:

First, construct a code that solves the maximization problem for each of the above three functions and computes the maximized value of the function and the corresponding optimal $r$ and $k$ (Obviously, these will be a function of $s$ and $d$.) For instance, for the first case of $f^1$, let me denote them by $f^{1*}$, $r^{1*}$, and $k^{1*}$; and the same for the cases of $f^2$ and $f^3$.

Second, using Plot3D, plot $f^{1*}$, $f^{2*}$, and $f^{3*}$ in one same diagram against $d \in [0,1]$ and $s \in [0,2]$; similarly, plot $r^{1*}$, $r^{2*}$, and $r^{3*}$ in one same diagram against $d \in [0,1]$ and $s \in [0,2]$, and plot $k^{1*}$, $k^{2*}$, and $k^{3*}$ in one same diagram against $d \in [0,1]$ and $s \in [0,2]$.

Lastly, since the first diagram allows comparing $f^{1*}$, $f^{2*}$, and $f^{3*}$, I can get the maximum value of the objective function $f$, and by looking at the second and third diagrams, I can get the corresponding optimal $r$ and $k$, which are what I'm ultimately looking for.

My Mathematica code for doing this is as follows (which is an extension of Alex Trounev's answer in 3D-Plot optimization results for varying parameter values):

Block[{t = 0}, f1 = (k^2 r (3 d + k r) + 3 (d^2 - k^2 r) s - 6 d s^2 + 3 s^3)/(6 s^2); max1 = Flatten[Table[{d, s, MaxValue[{f1, 0 <= d <= 1, s >= 2 d, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxk1 = Flatten[Table[{d, s, k /. Last@Maximize[{f1, 0 <= d <= 1, s >= 2 d, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxr1 = Flatten[Table[{d, s, r /. Last@Maximize[{f1, 0 <= d <= 1, s >= 2 d, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1];] Block[{t = 0}, f2 = 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)); max2 = Flatten[Table[{d, s, MaxValue[{f2, 0 <= d <= 1, s >= 2 d, d <= k <= s, 0 <= r < d/k}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxk2 = Flatten[Table[{d, s, k /. Last@Maximize[{f2, 0 <= d <= 1, s >= 2 d, d <= k <= s, 0 <= r < d/k}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxr2 = Flatten[Table[{d, s, r /. Last@Maximize[{f2, 0 <= d <= 1, s >= 2 d, d <= k <= s, 0 <= r < d/k}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1];] Block[{t = 0}, f3 = (d^3 + 3 d^2 k r + k^3 (-1 + r)^2 r + 3 r s^3 + 3 d r (k^2 (-1 + r) - 2 s^2))/(6 r s^2); max3 = Flatten[Table[{d, s, MaxValue[{f3, 0 <= d <= 1, s >= 2 d, d <= k <= s, d/k <= r <= 1}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxk3 = Flatten[Table[{d, s, k /. Last@Maximize[{f3, 0 <= d <= 1, s >= 2 d, d <= k <= s, d/k <= r <= 1}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1]; maxr3 = Flatten[Table[{d, s, r /. Last@Maximize[{f3, 0 <= d <= 1, s >= 2 d, d <= k <= s, d/k <= r <= 1}, {k, r}]}, {d, 0, 1, .1}, {s, 0, 2, .1}], 1];] {ListPlot3D[max1, max2, max3, AxesLabel -> {"d", "s", "V"}], ListPlot3D[maxk1, maxk2, maxk3, PlotRange -> {0, 2}, AxesLabel -> {"d", "s", "k"}], ListPlot3D[maxr1, maxr2, maxr3, PlotRange -> {0, 1}, AxesLabel -> {"d", "s", "r"}]} 

(In the code, you can simply by-pass the meaning of $t=0$.)

When this code is run, I get a list of error messages. Can anyone help? Thanks so much, in advance!

$\endgroup$
  • 1
    $\begingroup$ When s = 0 you will divide by zero. You have some semicolons missing between statements (mathematica will interpret this to mean you want to multiply). ListPlot3D takes all the data in a single argument, not a sequence, i.e. ListPlot3D[{a,b,c}] not ListPlot3D[a, b, c]. $\endgroup$ – wxffles Jul 9 at 23:25
  • $\begingroup$ Thanks, wxffles! $\endgroup$ – ppp Jul 11 at 2:35
1
$\begingroup$

Both variables d, s must be separated from 0. Then there are no messages.

Block[{t = 0, d0 = 10^-6, s0 = 10^-6}, 
 f1 = (k^2 r (3 d + k r) + 3 (d^2 - k^2 r) s - 6 d s^2 + 
     3 s^3)/(6 s^2); 
 max1 = Flatten[
   Table[{d, s, 
     MaxValue[{f1, 0 <= d <= 1, s >= 2 d, 0 <= k < d, 
       0 <= r <= 1}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1];
 maxk1 = Flatten[
   Table[{d, s, 
     k /. Last@
       Maximize[{f1, 0 <= d <= 1, s >= 2 d, 0 <= k < d, 
         0 <= r <= 1}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1]; 
 maxr1 = Flatten[
   Table[{d, s, 
     r /. Last@
       Maximize[{f1, 0 <= d <= 1, s >= 2 d, 0 <= k < d, 
         0 <= r <= 1}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1];]

Block[{t = 0, d0 = 10^-6, s0 = 10^-6}, 
 f2 = 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)); 
 max2 = Flatten[
   Table[{d, s, 
     MaxValue[{f2, 0 <= d <= 1, s >= 2 d, d <= k <= s, 
       0 <= r < d/k}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1]; 
 maxk2 = Flatten[
   Table[{d, s, 
     k /. Last@
       Maximize[{f2, 0 <= d <= 1, s >= 2 d, d <= k <= s, 
         0 <= r < d/k}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1];
  maxr2 = Flatten[
   Table[{d, s, 
     r /. Last@
       Maximize[{f2, 0 <= d <= 1, s >= 2 d, d <= k <= s, 
         0 <= r < d/k}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1];]

Block[{t = 0, d0 = 10^-6, s0 = 10^-6}, 
 f3 = (d^3 + 3 d^2 k r + k^3 (-1 + r)^2 r + 3 r s^3 + 
     3 d r (k^2 (-1 + r) - 2 s^2))/(6 r s^2); 
 max3 = Flatten[
   Table[{d, s, 
     MaxValue[{f3, 0 <= d <= 1, s >= 2 d, d <= k <= s, 
       d/k <= r <= 1}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 1]; 
 maxk3 = Flatten[
   Table[{d, s, 
     k /. Last@
       Maximize[{f3, 0 <= d <= 1, s >= 2 d, d <= k <= s, 
         d/k <= r <= 1}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 
   1]; maxr3 = 
  Flatten[Table[{d, s, 
     r /. Last@
       Maximize[{f3, 0 <= d <= 1, s >= 2 d, d <= k <= s, 
         d/k <= r <= 1}, {k, r}]}, {d, d0, 1, .1}, {s, s0, 2, .1}], 
   1];] 

{ListPlot3D[{max1, max2, max3}, AxesLabel -> {"d", "s", "V"}], 
 ListPlot3D[{maxk1, maxk2, maxk3}, PlotRange -> {0, 2}, 
  AxesLabel -> {"d", "s", "k"}], 
 ListPlot3D[{maxr1, maxr2, maxr3}, PlotRange -> {0, 1}, 
  AxesLabel -> {"d", "s", "r"}]}

Figure 1

$\endgroup$
  • $\begingroup$ Thanks, Alex, as always! It was very helpful, this time as well!! $\endgroup$ – ppp Jul 11 at 2:33
  • $\begingroup$ @ppp you're welcome! $\endgroup$ – Alex Trounev Jul 11 at 2:40
  • $\begingroup$ May I ask you a followup question? I'm trying to add one more complication. Regarding the conditions in the second Block[] for $f_2$, I would like to make a slight change as follows: $k=d$ if $r=0$, and $d <k\leq s$ if $0 < r < \frac{d}{k}$. For this, I replaced d <= k <= s, 0 <= r < d/k by If[r == 0, k = d, d < k <= s]. And I don't get a result. I also tried k==d instead of k=d. It gives some result, which however doesn't seem right. Can you please help me once more? I really appreciate! $\endgroup$ – ppp Jul 11 at 2:45
  • $\begingroup$ @ppp A special case would be better to allocate a separate block, and then combine the results. $\endgroup$ – Alex Trounev Jul 11 at 9:21
  • $\begingroup$ Yes, it is working! Thanks so much, Alex!! You are so great! $\endgroup$ – ppp Jul 11 at 17:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.