# 3D-plot of a maximization of a definite integral

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, $$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"}]}

• 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!! Jul 15, 2019 at 1:29
• 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?
– ppp
Jul 15, 2019 at 1:53
• 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. Jul 15, 2019 at 3:22
• @ppp Integrals are not calculated, so the optimization process does not end. Jul 15, 2019 at 5:14
• @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. Jul 15, 2019 at 15:11