mathematical optimization – Possible bug in Maximize

I have a maximization problem with several inequality constraints. (code below) I expected one of the constraints to be binding. Maximize gave a solution where it wasn’t binding so I was suspicious. Moreover, the maximum value Maximize finds is negative, while it is very easy to get 0 using a corner solution. I tried Nmaximize with the same objective and the same constraints, it basically gave the numbers I expected. The maximum value it returned is positive and argmax is clearly in the feasible region.

Here is the problem in its simplest form I can reduce it to. You can skip to the code below.

Suppose $F$ is a CDF over $(0,1)$. In fact, suppose it is the uniform distribution $F(x)=x$. Then, the objective function is:

$$ pp=-F(underline{theta}_s) dfrac{F(overline{theta}_b)-F(underline{theta}_b)}{2(F(overline{theta}_b)-F(underline{theta}_b) + F(overline{theta}_s)-F(underline{theta}_s))} left( mathbb{E}(theta| thetain (underline{theta}_b,overline{theta}_b) ) – underline{theta}_s right) – (1-F(overline{theta}_s)) dfrac{F(overline{theta}_b)-F(underline{theta}_b)}{2(F(overline{theta}_b)-F(underline{theta}_b) + F(overline{theta}_s)-F(underline{theta}_s))} dfrac{F(overline{theta}_b) -F(overline{theta}_s)}{F(overline{theta}_b) -F(underline{theta}_b)} left( mathbb{E}(theta| thetain (overline{theta}_s,overline{theta}_b) ) – overline{theta}_s right) – F(underline{theta}_b) dfrac{F(overline{theta}_s)-F(underline{theta}_s)}{2(F(overline{theta}_b)-F(underline{theta}_b) + F(overline{theta}_s)-F(underline{theta}_s))} dfrac{F(underline{theta}_b) -F(underline{theta}_s)}{F(overline{theta}_s) -F(underline{theta}_s)} left( underline{theta}_b – mathbb{E}(theta| thetain (underline{theta}_s,underline{theta}_b) ) right) – (1-F(overline{theta}_b)) dfrac{F(overline{theta}_s)-F(underline{theta}_s)}{2(F(overline{theta}_b)-F(underline{theta}_b) + F(overline{theta}_s)-F(underline{theta}_s))} left( underline{theta}_b – mathbb{E}(theta| thetain (underline{theta}_s,overline{theta}_s) ) right) – underline{theta}_sF(underline{theta}_s) + overline{theta}_b(1-F(overline{theta}_b)) $$

Definitions in Mathematica (I gave the problem in a slightly different form because otherwise the conditioning probability in the expectation can return $frac{1}{0}$ error when the bounds of the integration are equal):

F(x_):= x

pp(ls_, lb_, us_, 
  ub_) := -F(ls)*((F(ub) - F(lb))/(
    F(ub) - F(lb) + F(us) - F(ls)))*((!(
*SubsuperscriptBox(((Integral)), (lb), (ub))(
*FractionBox((x), ((F(ub) - F(lb)))) (F')(
        x) (DifferentialD)x)) - ls)/
    2) - (1 - F(us))*((F(ub) - F(lb))/(
    F(ub) - F(lb) + F(us) - F(ls)))*(1/(F(ub) - F(lb)))*((!(
*SubsuperscriptBox(((Integral)), (us), (ub))(x*(F')(
        x) (DifferentialD)x)) - us*(F(ub) - F(us)))/2) - 
  F(lb)*((F(us) - F(ls))/(F(ub) - F(lb) + F(us) - F(ls)))*(1/(
    F(us) - F(ls)))*((lb*(F(lb) - F(ls)) - !(
*SubsuperscriptBox(((Integral)), (ls), (lb))(x*(F')(
        x) (DifferentialD)x)))/2) - (1 - F(ub))*((F(us) - F(ls))/(
    F(ub) - F(lb) + F(us) - F(ls)))*((ub - !(
*SubsuperscriptBox(((Integral)), (ls), (us))(
*FractionBox((x), ((F(us) - F(ls)))) (F')(
        x) (DifferentialD)x)))/2) - ls*F(ls) + ub*(1 - F(ub))

Given these, I tried the following maximization problem:

Maximize({pp, 0 <= ls <= lb <= ub <= 1, ls <= us <= ub, 
  ls >= 1 - ub}, {ls, us, lb, ub})

This is the output:

{pp, {ls -> 1/4, us -> 9/16, lb -> 9/16, ub -> 7/8}}

However,

pp(1/4, 9/16, 9/16, 7/8)

returns

-(37/1024).

For example,

pp(0, 0, 1, 1)

returns 0.

Moreover,

NMaximize({-F(ls)*((F(ub) - F(lb))/(
     F(ub) - F(lb) + F(us) - F(ls)))*((!(
*SubsuperscriptBox(((Integral)), (lb), (ub))(
*FractionBox((x), ((F(ub) - F(lb)))) (F')(
         x) (DifferentialD)x)) - ls)/
     2) - (1 - F(us))*((F(ub) - F(lb))/(
     F(ub) - F(lb) + F(us) - F(ls)))*(1/(F(ub) - F(lb)))*((!(
*SubsuperscriptBox(((Integral)), (us), (ub))(x*(F')(
         x) (DifferentialD)x)) - us*(F(ub) - F(us)))/2) - 
   F(lb)*((F(us) - F(ls))/(F(ub) - F(lb) + F(us) - F(ls)))*(1/(
     F(us) - F(ls)))*((lb*(F(lb) - F(ls)) - !(
*SubsuperscriptBox(((Integral)), (ls), (lb))(x*(F')(
         x) (DifferentialD)x)))/2) - (1 - F(ub))*((F(us) - F(ls))/(
     F(ub) - F(lb) + F(us) - F(ls)))*((ub - !(
*SubsuperscriptBox(((Integral)), (ls), (us))(
*FractionBox((x), ((F(us) - F(ls)))) (F')(
         x) (DifferentialD)x)))/2) - ls*F(ls) + ub*(1 - F(ub)), 
  0 <= ls <= lb <= ub <= 1, ls <= us <= ub, ls >= 1 - ub}, {ls, us, 
  lb, ub})

returns

{0.09375, {ls -> 0.249987, us -> 0.750013, lb -> 0.249987, 
  ub -> 0.750013}}

which has a positive value. I expected the argmax to be $frac{1}{4},frac{1}{4},frac{3}{4},frac{3}{4}$ so this was what I was hoping to get.

So am I missing something critical here or is this really a bug? I am on version 12.0 student edition, if it matters.