NMinimize go out of bounds on some linear constraints

Question

I'm trying to minimize a non-linear function of four variables with some linear constraints. NMinimize on Mathematica 8 seems to go out of bounds on some constraints giving complex values of the function at some point in the iteration. The function to minimize is

ff[lxw_, lwz_, c_, d_] := - J1 (lxw + lwz) - 2 J2 c +
T (-Log[2] - 1/2 (1 - lxw) Log[(1 - lxw)/4] - 
1/2 (1 + lxw) Log[(1 + lxw)/4] - 
1/2 (1 - lwz) Log[(1 - lwz)/4] - 
1/2 (1 + lwz) Log[(1 + lwz)/4] + 1/2 (1 - d) Log[(1 - d)/16] + 
1/8 (1 + 2 c + d - 2 lwz - 2 lxw) Log[
1/16 (1 + 2 c + d - 2 lwz - 2 lxw)])

where

T = 10;
J1 = 1;
J2 = -0.2;

are constant parameters. The constraints are

var = {lxw, lwz, c, d};
cons = And @@ Cases[ff @@ var, Log[x_] -> x > 0, Infinity] // Simplify

which amounts to

d < 1 && lwz < 1 && 1 + lwz > 0 && 1 + 2 c + d > 2 (lwz + lxw) && 
lxw < 1 && 1 + lxw > 0

Then

NMinimize[{ff @@ var, cons}, var, Method -> "DifferentialEvolution"]

gives

NMinimize::nrnum: The function value 7.60939\[VeryThinSpace]-3.87314 I
is not a real number at {c,d,lwz,lxw} =  {-0.267966,0.319033,0.899803,-0.0151082}. >>
{-12.5741, {c -> -0.236255, d -> -0.978425, lwz -> -0.681637, 
lxw -> -0.939595}}

Answer 1

There are two things that can be done to improve matters. One is to use constraints with weak inequalities that keep the log arguments above a minmium threshold. The other is to use an altered logarithm that allows argumens less-equal to zero and simply returns a suitable large negative.

T = 10;
J1 = 1;
J2 = -0.2;
bigValue = 10^6;
myLog[x_?NumberQ] := If[x <= 0, -bigValue, Log[x]]
vars = {lxw, lwz, c, d};
ff[lxw_, lwz_, c_, 
  d_] := -J1 (lxw + lwz) - 2 J2 c + 
   T (-Log[2] - 1/2 (1 - lxw) Log[(1 - lxw)/4] - 
      1/2 (1 + lxw) Log[(1 + lxw)/4] - 
      1/2 (1 - lwz) Log[(1 - lwz)/4] - 
      1/2 (1 + lwz) Log[(1 + lwz)/4] + 1/2 (1 - d) Log[(1 - d)/16] + 
      1/8 (1 + 2 c + d - 2 lwz - 2 lxw) Log[
        1/16 (1 + 2 c + d - 2 lwz - 2 lxw)]) /. Log -> myLog
eps = 1/1000;

So the constraints are now as follows.

cons = And @@ Cases[ff @@ var, myLog[x_] :> x >= eps, Infinity];

(* ut[57]= (1 - d)/16 >= 1/1000 && (1 - lwz)/4 >= 1/1000 && (1 + lwz)/
  4 >= 1/1000 && 
 1/16 (1 + 2 c + d - 2 lwz - 2 lxw) >= 1/1000 && (1 - lxw)/4 >= 1/
  1000 && (1 + lxw)/4 >= 1/1000 *)

Here is the optimization.

{min, vals} = 
 NMinimize[Evaluate[{ff @@ vars, cons}], vars, 
  Method -> "DifferentialEvolution"]

(* Out[55]= {-27.7261, {lxw -> -0.996, lwz -> 0.996, c -> 4.33532, 
  d -> -4.65512}} *)

Notice that the variables lxw and lwz really want to live at the boundary of the allowed region. This makes me suspect there is an issue with the forumlation, that maybe it's not solving the problem you want it to handle? Or maybe you want to let them be -1 and 1 respectively, and remove terms with logs that would then contain zero? This seems plausible since they have factors that are also zero and, in a limiting sense the values "should" be zero. When I make those adjustments, (skipping the code changes and going right to the results) I get

{min, vals} = 
 NMinimize[Evaluate[{ff @@ vars, cons}], vars, 
  Method -> "DifferentialEvolution"]

(* Out[99]= {-28.0146, {c -> 4.33553, d -> -4.65527}} *)