PenaltyFunction Option in NMinimize

Question

I am curious if anybody has some experience in using the built-in Option "PenaltyFunction" in any of the numerical optimization functions like NMinimize.

According to the documentation the user should be able to provide NMinimize with a custom penalty function to control the result of the minimization. The default option value is Automatic. I tried to find some examples of how to use this Option, but there is hardly any information available. The official documentation gives no example.

In the this Google.group post I found the following example, which is not working any more in Version 10.0.1.

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

(* {55., {x -> 20, y -> 19, z -> 16}} *)

now with "PenaltyFunction"

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

This gives a warning of the solution not meeting the constraints and the result (which I don't understand):

(* {22., {x -> 20, y -> 1, z -> 1}}*)

I tried also to play with the setting myself:

From the documentation I took:

 NMinimize[x + y, {x, y} ∈ Disk[]]
 Show[ContourPlot[x + y, {x, y} ∈ Disk[]], 
       Graphics[{Red, PointSize[Large], Point[{x, y} /. Last[%]]}]]

Now with some PenaltyFunction

NMinimize[x + y, {x, y} ∈ Disk[], Method -> {"NelderMead", 
           "PenaltyFunction" -> (Min[#, 0] &)}]
Show[ContourPlot[x + y, {x, y} ∈ Disk[]], 
      Graphics[{Red, PointSize[Large], Point[{x, y} /. Last[%]]}]]

or

NMinimize[x + y, {x, y} ∈ Disk[],  Method -> {"NelderMead", 
           "PenaltyFunction" -> ((# - Round[#])^3 &)}]
Show[ContourPlot[x + y, {x, y} ∈ Disk[]], 
      Graphics[{Red, PointSize[Large], Point[{x, y} /. Last[%]]}]]

I can't decide whether I am too stupid or just too lazy to figure out, what "PenaltyFunction" actually does. Here my questions:

• What arguments does "PenaltyFunction" use/accept?
• How is it possible to penalize individual fitting parameters?
• Do you have any example use cases that shed light on the whole issue?

EDIT

I found another Little Piece of Information here:

The author states that the Default Setting for "DifferentialEvolution" is

"PenaltyFunction"-> Function[{d,i},d*10^(4*i)]

Function applied to penalize invalid Parameter values outside constraints (d =distance from allowed value, i =number of iteration)

To me it is unclear if this applies to all constraints and how this is usefull.

Answer 1

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}}  *)