0
$\begingroup$

I'm applying the answer (If method) suggested by 'ThatGravityGuy' in this post in a more complicated case.

Consider my objective function (pardon its formidable expression):

$f=\int _d^{d+k r}\int _k^s A dydx+ \int _{d+k r}^s\int _k^s B dydx+ \int _d^{d+k r}\int _{\frac{d+k r-x}{r}}^k D dydx+ \int _{d+k r}^s\int _0^k E dydx$

where

$A = \frac{(1-t) (r (x-k)+k r+(1-r) x -J)-(d-t)}{s^2}$ if $A > 0$, while $A=0$ if $A\leq 0$,

$B = \frac{(1-t) (r (x-k)+k r+(1-r) x-J)-(d-t)}{s^2}$ if $B > 0$, while $B=0$ if $B\leq 0$,

$D = \frac{(1-t) (r (x-k)+(1-r) x+r y-J)-(d-t)}{s^2}$ if $D > 0$, while $D=0$ if $D\leq 0$,

$E = \frac{(1-t) (r (x-k)+(1-r) x+r y-J)-(d-t)}{s^2}$ if $E > 0$, while $E=0$ if $E\leq 0$,

$J =c \left(\frac{(q-1) \left(d^2+k^2 r\right)}{s^2}+\frac{d s+\frac{k^2 r}{2}}{s^2}\right)-\frac{(1-d) \left(d s+\frac{k^2 r}{2}\right)}{s^2}$

I would like to maximize my objective function $f$ with respect to $r$ and $k$ and Plot3D them, i.e. the maximum value of $f$, and the corresponding $r$ and $k$, against $c \in [0,1]$ and $q \in [1,2]$ under the following parameter values and constraints: $s=2$, $d=0.8$, $t=0$, $0 \leq k<d$, $0 \leq r \leq 1$, $0\leq c \leq 1$, $1 \leq q \leq 2$, $cq \leq 1$.

My Mathematica Code is as follows (Since C, D, and E are protected letters, I'm starting with unprotecting them):

s = 2; d = 0.8; t = 0; 
J = -(((1 - d) ((k^2 r)/2 + d s))/s^2) + c (((-1 + q) (d^2 + k^2 r))/s^2 + ((k^2 r)/2 + d s)/s^2);  
A = If[(((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2 > 0, (((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2, 0]; 
B = If[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2 > 0, (((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, 0]; 
D1 = If[(((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2 > 0, (((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2, 0]; 
E1 = If[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2 > 0, (((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, 0]; 
max = Flatten[Table[{C1, Q, MaxValue[{Integrate[A /. {c -> C1, q -> Q}, {x, d, d + r*k}, {y, (d + r*k - x)/r, k}] + Integrate[B /. {c -> C1, q -> Q}, {x, d, d + r*k}, {y, k, s}] + Integrate[D1 /. {c -> C1, q -> Q}, {x, d + r*k, s}, {y, 0, k}] + Integrate[E1 /. {c -> C1, q -> Q}, {x, d + r*k, s}, {y, k, s}], 0 <= C1 <= 1, 1 <= Q <= 2, C1*Q <= 1, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {C1, 0, 1, .1}, {Q, 1, 2, .1}], 1]; 
maxk = Flatten[Table[{C1, Q, k /. Last@Maximize[{Integrate[A /. {c -> C1, q -> Q}, {x, d, d + r*k}, {y, (d + r*k - x)/r, k}] + Integrate[B /. {c -> C1, q -> Q}, {x, d, d + r*k}, {y, k, s}] + Integrate[D1 /. {c -> C1, q -> Q}, {x, d + r*k, s}, {y, 0, k}] + Integrate[E1 /. {c -> C1, q -> Q}, {x, d + r*k, s}, {y, k, s}], 0 <= C1 <= 1, 1 <= Q <= 2, C1*Q <= 1, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {C1, 0, 1, .1}, {Q, 1, 2, .1}], 1]; 
maxr = Flatten[Table[{C1, Q, r /. Last@Maximize[{Integrate[A /. {c -> C1, q -> Q}, {x, d, d + r*k}, {y, (d + r*k - x)/r, k}] + Integrate[B /. {c -> C1, q -> Q}, {x, d, d + r*k}, {y, k, s}] + Integrate[D1 /. {c -> C1, q -> Q}, {x, d + r*k, s}, {y, 0, k}] + Integrate[E1 /. {c -> C1, q -> Q}, {x, d + r*k, s}, {y, k, s}], 0 <= C1 <= 1, 1 <= Q <= 2, C1*Q <= 1, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {C1, 0, 1, .1}, {Q, 1, 2, .1}], 1]; 
{ListPlot3D[max, AxesLabel -> {"c", "q", "f"}], ListPlot3D[maxk, PlotRange -> {0, 2}, AxesLabel -> {"c", "q", "k"}], ListPlot3D[maxr, PlotRange -> {0, 1}, AxesLabel -> {"c", "q", "r"}]}

This code is running forever. Any help will be greatly appreciated!

Alex Trounev's suggested approach is working! Here is the code:

s = 2; d = 0.8; t = 0;
J = -(((1 - d) ((k^2 r)/2 + d s))/s^2) + c (((-1 + q) (d^2 + k^2 r))/s^2 + ((k^2 r)/2 + d s)/s^2);
regA = ImplicitRegion[1/800 (-16 c (3 + 2 q) + 25 c k^2 (1 - 2 q) r + 8 (-18 + 25 x) + 5 r ((-40 + k) k + 40 y)) > 0 &&     d <= x <= d + r*k && (d + r*k - x)/r <= y <= k, {x, y}];
regB = ImplicitRegion[1/800 (-144 - 16 c (3 + 2 q) + 5 k^2 r + 25 c k^2 (1 - 2 q) r + 200 x) > 0 && d <= x <= d + r*k && k <= y <= s, {x, y}];
regD = ImplicitRegion[1/800 (-16 c (3 + 2 q) + 25 c k^2 (1 - 2 q) r + 8 (-18 + 25 x) + 5 r ((-40 + k) k + 40 y)) > 0 && d + r*k <= x <= s && 0 <= y <= k, {x, y}];
regE = ImplicitRegion[1/800 (-144 - 16 c (3 + 2 q) + 5 k^2 r + 25 c k^2 (1 - 2 q) r + 200 x) > 0 && d + r*k <= x <= s && k <= y <= s, {x, y}];
max = Flatten[Table[{c, q, MaxValue[{Integrate[(((1 - r) x + r (x - k) + r y - J) (1 - t) - (d - t))/s^2, {x, y} \[Element] regA] + Integrate[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regB] Integrate[(((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regD] + Integrate[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regE], 0 <= c <= 1, 1 <= q <= 2, c*q <= 1, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {c, 0, 1, .1}, {q, 1, 2, .1}], 1];
maxk = Flatten[Table[{c, q, k /. Last@Maximize[{Integrate[(((1 - r) x + r (x - k) + r y - J) (1 - t) - (d - t))/s^2, {x, y} \[Element] regA] + Integrate[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regB] + Integrate[(((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regD] + Integrate[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regE], 0 <= c <= 1, 1 <= q <= 2, c*q <= 1, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {c, 0, 1, .1}, {q, 1, 2, .1}], 1];
maxr = Flatten[Table[{c, q, r /. Last@Maximize[{Integrate[(((1 - r) x + r (x - k) + r y - J) (1 - t) - (d - t))/s^2, {x, y} \[Element] regA] + Integrate[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regB] + Integrate[(((1 - r)*x + r*(x - k) + r*y - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regD] + Integrate[(((1 - r)*x + r*(x - k) + r*k - J)*(1 - t) - (d - t))/s^2, {x, y} \[Element] regE], 0 <= c <= 1, 1 <= q <= 2, c*q <= 1, 0 <= k < d, 0 <= r <= 1}, {k, r}]}, {c, 0, 1, .1}, {q, 1, 2, .1}], 1];
{ListPlot3D[max, AxesLabel -> {"c", "q", "f"}], ListPlot3D[maxk, PlotRange -> {0, 2}, AxesLabel -> {"c", "q", "k"}], ListPlot3D[maxr, PlotRange -> {0, 1}, AxesLabel -> {"c", "q", "r"}]}
$\endgroup$
  • $\begingroup$ I would recommend using esc+dsD+esc double sided symbols—it is generally recommended to not unprotect and redefine the symbols you used, especially D and E. Hope this helps!! $\endgroup$ – CA Trevillian Jul 15 at 1:29
  • $\begingroup$ Thanks, CA Trevillian! I tried, for example, both esc+dsD+esc and \[DoubleStruckCapitalD] to format, but didn't work. (They did work in Mathematica.) Do you know how to do that to format? $\endgroup$ – ppp Jul 15 at 1:53
  • $\begingroup$ I'm not following, they work in Mathematica, and don't work where? You are doing this in Mathematica? The shortcuts you listed should work in Mathematica, as you stated. When you redefine those symbols, functions that rely on them can have issues. On the surface, this is what I gather is happening or what is causing issues. $\endgroup$ – CA Trevillian Jul 15 at 3:22
  • 2
    $\begingroup$ @ppp Integrals are not calculated, so the optimization process does not end. $\endgroup$ – Alex Trounev Jul 15 at 5:14
  • 3
    $\begingroup$ @ppp (Replace first of all C->C1, D->D1,E->E1 ) It is necessary for each integral to determine the region of integration, for example reg1 = ImplicitRegion[1/4 (-(4/5) + k r + 1/20 (8/5 + (k^2 r)/2) - c (1/4 (8/5 + (k^2 r)/2) + 1/4 (-1 + q) (16/25 + k^2 r)) + (1 - r) x + r (-k + x)) > 0 && d <= x <= d + r*k && (d + r*k - x)/r <= y <= k, {x, y}];. Then instead of Integrate[A... there will be Integrate[(((1 - r) x + r (x - k) + r k - J) (1 - t) - (d - t))/s^2, {x, y} \[Element] reg1]. Then the integral is calculated. $\endgroup$ – Alex Trounev Jul 15 at 15:11

Your Answer

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

Browse other questions tagged or ask your own question.