3
$\begingroup$

I need to maximize the following function (the input to NMaximize below)

    NMaximize[{((1/3) + f[p]*(1 - p)*(((1/(Sqrt[2]/2))*p)^2-1))/(1+(1-p)*
    (((1/(Sqrt[2]/2))*p)^2 - 1)), p >= Sqrt[2]/2, p <= 1}, {p}]

where $f(\cdot)$ is defined as follows

    f[p_?NumericQ] := NMinimize[{-((a + a^2 + b - 2 a b + b^2 + c - 2 a c - 2 b c + c^2)
    /((-1 + a) (a + b + c))), 0 <= a, a <= b, b <= c, c <= 1, c <= a + b, (a + b + c)/3 <= p}, 
    {a, b, c}, Method -> "DifferentialEvolution"][[1]]

However, as expected, the computation does not end ,and there are several alert messages. Considering the underlying mathematical problem I am trying to solve, I could replace $f(\cdot)$ by a simple function $g(\cdot)$ that approximates it, but I need $g(p)\le f(p)$ for all $p\in[0,1]$. I tried to use InterpolatingPolynomial with a few values of $f(\cdot)$. However, $f(\cdot)$ is neither concave nor convex in $[0,1]$. I am struggling to obtain a good approximation of $f(\cdot)$ which satisfies the $g(p)\le f(p)$ in $[0,1]$.

How can we generate such approximation function $g(\cdot)$ (alternately, can we solve this optimization problem in a different way)?

Improve this question
$\endgroup$
0

4 Answers 4

Reset to default
4
$\begingroup$

The extremum of a polytope usually saturates a few of its inequalities, in this case $(a + b + c)/3 = p$ and $b=c$. With this we can find $f(p)$ exactly:

Minimize[{-((a + a^2 + b - 2 a b + b^2 + c - 2 a c - 2 b c + 
          c^2)/((-1 + a) (a + b + c))),
          0 <= a, a <= b, b == c, c <= 1, c <= a + b, (a + b + c)/3 == p},
         {a, b, c}]
(*    complicated result, the important part of which is    *)

$$ f(p) = \frac{2(p-1+\sqrt{1-p})}{p} $$ from which the exact answer follows:

Maximize[{(1/3 - (2 (-1 + p) (-1 + Sqrt[1 - p] + p) (-1 + 2 p^2))/p)/(p + 2 p^2 - 2 p^3),
         1/Sqrt[2] < p < 1}, p] // RootReduce

(*    {0.348695341134874, {p -> 0.8151520693995784}}    *)

These numbers are given as Root objects, which I've converted to numerical representation here with N.

Improve this answer
$\endgroup$
9
  • $\begingroup$ Thank you @Roman , this solution is very interesting! Could you please explain me how we can be sure that for this specific optimization problem the minimum is attained when $b=c$ (i.e., that the inequality saturated is precisely $b\le c$) ? I have the feeling there is something obvious I am missing... I would be very happy if you could provide a brief clarification about it. Thanks. $\endgroup$
    – Penelope Benenati
    Dec 10, 2020 at 14:15
  • 1
    $\begingroup$ @PenelopeBenenati I have no proof but combined two observations non-rigorously: (1) the extrema of a function in a polytope-shaped region are usually on the polytope's surface, not its interior, which means that some of your inequalities will be saturated by the solution, and (2) by inspection I noticed that the numerical solutions $f(p)$ indeed saturate $b=c$ and $a+b+c=3p$. Concerning (1): there are numerous counter-examples, so it's just a hint; for example, $x^2$ on $[-1,1]$ is minimal in the interior, and maximal on the surface. $\endgroup$
    – Roman
    Dec 10, 2020 at 15:41
  • $\begingroup$ Thank you @Roman, I see. I am trying now to rigorously prove this claim, so that I will be able to use this solution to have an analytical result (which is one of the reasons I like it). Maybe it is not difficult. $\endgroup$
    – Penelope Benenati
    Dec 10, 2020 at 15:46
  • $\begingroup$ I tried your solution, and I thought of obtaining a step-by-step solution either with Mathematica or Wolfram Pro Premium. However, it seems that the step-by-step solution is not present for minimization problems. Could you please explain me how to obtain the steps for solving this problem (in the way you showed) @Roman ? $\endgroup$
    – Penelope Benenati
    Dec 12, 2020 at 21:08
  • 1
    $\begingroup$ @PenelopeBenenati you're right; but we already know approximately where the optimum lies (from the numerical answers) and so I took a shortcut. $\endgroup$
    – Roman
    Dec 16, 2020 at 6:50
4
$\begingroup$

I make a slight adjustment to f[p], to restrict a so as to prevent the denominator from vanishing. Also it is faster in general to use FindMinValue and will likely work just as well for present purposes.

f[p_?NumericQ] := 
 FindMinValue[{-((a + a^2 + b - 2 a b + b^2 + c - 2 a c - 2 b c + 
        c^2)/((-1 + a) (a + b + c))), 0 <= a <= .99, a <= b, b <= c, 
   c <= 1, c <= a + b, (a + b + c)/3 <= p}, {a, b, c}]

I think the idea of using an interpolation is sound and will improve speed.

gvals = Table[{p, g[p]}, {p, Sqrt[2]/2., 1., .001}];
interp = Interpolation[gvals];

Now maximize:

In[2018]:= FindMaximum[{((1/3) + 
     interp[p]*(1 - p)*(((1/(Sqrt[2]/2))*p)^2 - 1))/(1 + (1 - 
        p)*(((1/(Sqrt[2]/2))*p)^2 - 1)), p >= Sqrt[2]/2, p <= 1}, {p}]

(* Out[2018]= {0.348695, {p -> 0.815155}} *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Nice, we confirm one another. $\endgroup$
    – user64494
    Dec 9, 2020 at 19:17
  • $\begingroup$ Thank you for the proposed solution. I am sorry, I did not understand how we can be sure that interp[p]$\le$f[p] for all $p$ in the considered interval, as required in the question text. $\endgroup$
    – Penelope Benenati
    Dec 9, 2020 at 21:41
  • $\begingroup$ Furthermore, in this problem I am also interested in the case $a$ approaches $1$. Do you think there is anything we can do for studying this case too? $\endgroup$
    – Penelope Benenati
    Dec 9, 2020 at 21:45
4
$\begingroup$

Realizing your idea concerning the interpolation and correcting your syntax (&& instead of , in the constraints), I obtain

f[p_?NumericQ] :=  NMinimize[{-((a + a^2 + b - 2 a b + b^2 + c - 2 a c - 2 b c + 
      c^2)/((-1 + a) (a + b + c))), 
 0 <= a && a <= b && b <= c && c <= 1 && 
  c <= a + b && (a + b + c)/3 <= p}, {a, b, c}, 
Method -> "DifferentialEvolution", AccuracyGoal -> 4, 
PrecisionGoal -> 4][[1]];
if = Interpolation[Table[{p, f[p]}, {p, 0.7, 1, 0.05}]];
NMaximize[{((1/3) + 
 if[p]*(1 - p)*(((1/(Sqrt[2]/2))*p)^2 - 1))/(1 + (1 - 
    p)*(((1/(Sqrt[2]/2))*p)^2 - 1)), p >= Sqrt[2]/2 && p <= 1}, {p}]

(*{0.348699, {p -> 0.815213}}*)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you for the correction and for the proposed solution. I did not understand how can we be sure that if[p]$\le$f[p] for all $p$ in the considered interval, as required in the question text. $\endgroup$
    – Penelope Benenati
    Dec 9, 2020 at 21:29
  • $\begingroup$ @PenelopeBeneti: I think the condition if[p]<=f[p] is not of importance for solving the problem under consideration.It is important that the difference RealAb[if[p]-f[p]] may be done uniformly arbitrary small on the interval $[0.7,1]$. The base of that statement needs in analysis of differential properties of f[p].. $\endgroup$
    – user64494
    Dec 10, 2020 at 7:41
  • $\begingroup$ I understand what you mean @user64494 and I thank you again because I could learn from your solution. In my question there is that constraint because in the proof I am writing, I need to use a lower bound of $f()$ to provide a lower bound of the maximum I am computing with NMaximize. In other words, I cannot state that the maximum is at least $0.348699$ if the approximation function if[] is, for some values, strictly larger than than $f()$, because the real maximum could be smaller (e.g., 0.348698; BTW Roman and Akku14 obtained 0.348695...). $\endgroup$
    – Penelope Benenati
    Dec 10, 2020 at 14:09
  • $\begingroup$ Finally, if one could verify that for a certain number of points (e.g., the ones you choose : {p, 0.7, 1, 0.05}), that we have if[x]<=f[x] for all $x$ in the considered interval, than I would be sure that the result of the maximization is lower bounded by $0.348699$. To further clarify, I would content myself to state that the maximization result is at least $0.347$, if was sure about it. $\endgroup$
    – Penelope Benenati
    Dec 10, 2020 at 14:11
3
$\begingroup$

Let me use a little trick to get an anylytical solution.

Since Minimize can't find minimum depending on variable p with (a + b + c)/3 <= p , set p to a mathematical constant (here p -> 2/E , that is within alowed p values.

I know, of course, you have to test later with other methods, wether this is valid, but i made experience, it often works well.

{min, varsmin} = 
  Minimize[{-((a + a^2 + b - 2 a b + b^2 + c - 2 a c - 2 b c + 
     c^2)/((-1 + a) (a + b + c))), 0 <= a, a <= b, b <= c,   c <= 1,
 c <= a + b, (a + b + c)/3 <= 2/E}, {a, b, c}];

f[p_] = min /. E -> 2/p // FullSimplify

(*   2 (1 + Sqrt[(1 - p)/p^2] - 1/p)   *)

Plot[f[p], {p, 1/Sqrt[2], 1}]

nmax = NMaximize[{((1/3) + 
  f[p]*(1 - p)*(((1/(Sqrt[2]/2))*p)^2 - 1))/(1 + (1 - 
     p)*(((1/(Sqrt[2]/2))*p)^2 - 1)), p >= Sqrt[2]/2, 
   p <= 1}, {p}]

(*   {0.348695, {p -> 0.815152}}   *)

Analytical solution

max = Maximize[{((1/3) + 
   f[p]*(1 - p)*(((1/(Sqrt[2]/2))*p)^2 - 1))/(1 + (1 - 
      p)*(((1/(Sqrt[2]/2))*p)^2 - 1)), p >= Sqrt[2]/2, 
 p <= 1}, {p}] // FullSimplify

(*   {Root[40292242 - 521152218 #1 + 2609177160 #1^2 - 6515590038 #1^3 + 
8910881883 #1^4 - 7184393418 #1^5 + 3522825570 #1^6 - 
1014076554 #1^7 + 147469221 #1^8 - 8264916 #1^9 + 34992 #1^10 &, 
3], {p -> 
Root[-48 - 296 #1 + 269 #1^2 + 2044 #1^3 - 2940 #1^4 - 60 #1^5 + 
  1176 #1^6 + 144 #1^8 - 432 #1^9 + 144 #1^10 &, 4]}}   *)

max // N[#, 50] &

(*   {0.34869534113486204014442063948484350125554473832604, {p -> 
 0.81515206939957835833133995379472132118970645053175}}   *)
Improve this answer
$\endgroup$
5
  • $\begingroup$ Sometimes that trick works, sometimes doesn't. I don't see how to ground it. $\endgroup$
    – user64494
    Dec 9, 2020 at 20:06
  • 1
    $\begingroup$ The reason why MInimize or Solve or Reduce in simular cases fails, is, that checking all conditions for a free variable p is quite difficult. Using a mathematical constant, you circumvent the built in conditional checking routines. Then of course you have to do that condtional checking by hand with other methods. Ultimate check is whether result is right. $\endgroup$
    – Akku14
    Dec 10, 2020 at 3:46
  • $\begingroup$ I like your approach @Akku14 . Unfortunately my software never ends running the first input command. Do you think it happens because I have Wolfram|Alpha Notebook Edition (wolfram.com/wolfram-alpha-notebook-edition)? I have another question. You provided an analytical solution as Roman. Why your solution is different from his solution (starting from the 14th digit)? $\endgroup$
    – Penelope Benenati
    Dec 10, 2020 at 17:06
  • 1
    $\begingroup$ @PenelopeBenenati, may be WolframAlpha is not as mighty as Mathematica in this case. Applying // N[#, 50] & to both, i find no difference between Roman's and my result. $\endgroup$
    – Akku14
    Dec 10, 2020 at 18:34
  • $\begingroup$ Thank you @Akku14 . I am afraid I need to buy Mathematica because Notebook is not sufficient... $\endgroup$
    – Penelope Benenati
    Dec 10, 2020 at 22:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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