# Optimizing a nested function in mathematica getting stuck at local maximum and giving wrong answer

I am trying to maximize this function over the value "p":

    g2ces[p_, theta_, I_, a_, c1_]  :=
theta*p*Min[I,
a*p^c1] + (theta*
theta)*(((I - Min[I, a*p^c1])/a)^(1/c1)*(I -
Min[I, a*p^c1]) + ((I/a)^(1/c1))*I)


At the following other parameter values (which I take as an example), we get:

FindMaximum[{g2ces[p, 0.5, 2, 2.5, -2]}, p]


p = 1.25. This is the correct solution.

And plotting this over various "p" values gives:

    Plot[Evaluate@
Table[g2ces[p, theta, 2, 2.5, -2], {theta, 0.5, 0.5,1}], {p, 0, 5}]


This solution is analytically derived and a special case. This function is equivalent to the following function, "VDPces", which is easier for me to generalize (and thus I want to work with this version of the function instead of the "g2ces" one):

Rev[p_, I_, a_, c1_] := p*Min[I, a*p^c1]
V2ces[I_, a_, c1_]  := MaxValue[{Rev[p, I, a, c1]}, p]
VDPces[p_, theta_, I_, a_, c1_] :=
theta*(Rev[p, I, a, c1] +
theta*V2ces[I - Min[I, a*p^c1], a, c1]) + (1 - theta)*(0 +
theta*V2ces[I, a, c1])


When I maxmimize this function, I get a different solution:

FindMaximum[{VDPces[p, 0.5, 2, 2.5, -2]}, p]


p = 1.11803 (incorrect)

and a different plot

Plot[Evaluate@
Table[VDPces[p, theta, 2, 2.5, -2], {theta, 0.5, 0.5, 1}], {p, 0,
5}]


that looks like:

Obviously something is going wrong with the optimization.

• is I the imaginary unit? – vi pa Jan 19 '20 at 23:56
• No I is just another variable, no imaginary numbers. Set it equal to 2. – wolfsatthedoor Jan 20 '20 at 1:15
• Nothing wrong with the maximization. The problem is the function definitions which are different. Show[Plot[{Evaluate[VDPces[p, 0.5, 2, 2.5, -2]], Evaluate[g2ces[p, 0.5, 2, 2.5, -2]]}, {p, 0, 2}]] – Cesareo Jan 20 '20 at 11:41
• I pretty confident that theoretically these should be identical though is my point. In the g2ces, I just find a closed form solution for that particular function (a*p^c1). – wolfsatthedoor Jan 29 '20 at 4:44