PenaltyFunction Option in NMinimize

As noted in the OP, "PenaltyFunction" has the form

    pf[..., step]

where the first argument is constructed from the constraint. When minimizing an objective function of the form obj[x,..], the actual penalty function is a scaled sum of penalties constructed from pf[] and the constraints:

Max[1, Abs[obj[x,..]]] * (  (* the objective function scales *)
        (* a sum of *)
        (*   1. penalties for constraints f[x,..] == g[x,..] *)
   pf[Abs[f[x,..] - g[x,..]], step] + ...
        (*   2. penalties for constraints f[x,..] <= g[x,..] *)
   pf[f[x,..] - g[x,..], step] (1 + Sign[f[x,..] - g[x,..]])/2
  )

The penalty function is nonnegative and is zero precisely on the domain define by the constraints, unless pf[] is zero or negative outside the domain. (Normally, pf[] should be positive.)

We can get the actual penalty function from an internal variable Optimization`NMinimizeDump`penalty (the current step is given by Optimization`NMinimizeDump`itr):

myObj[x_?NumericQ, y_?NumericQ] := x + y;
myPenalty[d_?NumericQ, step_?NumericQ] := Abs[d]^step;
NMinimize[myObj[x, y], {x, y} [Element] Disk[],
 Method -> {"NelderMead", "PenaltyFunction" -> myPenalty},
 StepMonitor :> 
  Block[{x, y, Optimization`NMinimizeDump`itr = itr}, 
   Return[Optimization`NMinimizeDump`penalty, NMinimize]]
 ]

(*
  Max[1, Abs[myObj[x, y]]] *
    myPenalty[-1 + x^2 + y^2, itr] * 
     (1 + Sign[-1 + x^2 + y^2])/2
*)

Both kinds of constraints:

myObj[x_?NumericQ, y_?NumericQ, z_?NumericQ] := x + y + z;
myPenalty[d_?NumericQ, step_?NumericQ] := Abs[d]^step;
NMinimize[{myObj[x, y, z], z == x^2 + y^2}, {x, y, z} [Element] 
  Ball[],
 Method -> {"NelderMead", "PenaltyFunction" -> myPenalty},
 StepMonitor :> 
  Block[{x, y, z, Optimization`NMinimizeDump`itr = itr}, 
   Return[Optimization`NMinimizeDump`penalty, NMinimize]]
 ]

(*
  Max[1, Abs[myObj[x, y, z]]] * (
   myPenalty[Abs[-x^2 - y^2 + z], itr] +
     myPenalty[-1 + x^2 + y^2 + z^2, itr] *
      (1 + Sign[-1 + x^2 + y^2 + z^2])/2
   )
*)

The default penalty function:

myObj[x_?NumericQ, y_?NumericQ] := x + y;
NMinimize[myObj[x, y], {x, y} [Element] Disk[],
 StepMonitor :> 
  Block[{x, y, Optimization`NMinimizeDump`itr = itr}, 
   Return[Optimization`NMinimizeDump`penalty, NMinimize]]
 ]

(*
  Max[1, Abs[myObj[x, y]]]*
   2^Floor[itr/10] (-1 + x^2 + y^2)^2 *
    (1 + Sign[-1 + x^2 + y^2])/2
*)

The Google group example

Here's the actual penalty function from the Google Group example in the OP:

NMinimize[{x + y + z, (1/20)*x + y + 5*z == 100,
  (x | y | z) [Element] Integers, 0 < x < 99, 0 < y < 99, 0 < z < 99},
 {x, y, z}, 
 Method -> {"SimulatedAnnealing", 
   "PenaltyFunction" -> (100*(#1 - Floor[#1]) &), 
   "SearchPoints" -> 250}, MaxIterations -> 500,
 StepMonitor :> (Hold[{x, y, z}] /. Round[v_] :> v /. 
    Hold[v_] :> Block[v, 
      Block[{Optimization`NMinimizeDump`itr = itr}, 
       Return[Optimization`NMinimizeDump`penalty, NMinimize]]])
 ]

(*
  Max[1, Abs[Round[x] + Round[y] + Round[z]]]* (
    100 (Abs[Round[x]/20 + Round[y] + 5 (-20 + Round[z])] -
        Floor[Abs[Round[x]/20 + Round[y] + 5 (-20 + Round[z])]]) +
     50 (-Floor[-Round[x]] - Round[x]) (1 + Sign[1 - Round[x]]) +
     50 (-Floor[-Round[y]] - Round[y]) (1 + Sign[1 - Round[y]]) +
     50 (-Floor[-Round[z]] - Round[z]) (1 + Sign[1 - Round[z]])
   )
*)

Note the last three terms are always zero because Floor[Round[u]] is equal to Round[u], when u is a real number. What's worse is that the first term is sometimes zero when the variables do not satisfy the constraint, whenever x is a multiple of 20 no matter what y or z are.

This suggests that removing Floor[#] from the "PenaltyFunction" should improve the performance of NMinimize, which it happens to do:

NMinimize[{x + y + z, (1/20)*x + y + 5*z == 100,
  (x | y | z) [Element] Integers, 0 < x < 99, 0 < y < 99, 0 < z < 99},
 {x, y, z}, 
 Method -> {"SimulatedAnnealing", "PenaltyFunction" -> (100*#1 &), 
   "SearchPoints" -> 250}, MaxIterations -> 500
 ]

(*  {43., {x -> 20, y -> 4, z -> 19}}  *)