18
$\begingroup$

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[%]]}]]

enter image description here

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[%]]}]]

enter image description here

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[%]]}]]

enter image description here

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.

$\endgroup$
9
  • $\begingroup$ Interesting question! I got it to work with 2 arguments : NMinimize[x^2 - y^2, {x, y} \[Element] Disk[], Method -> {"SimulatedAnnealing", "PenaltyFunction" -> ((#1^4 + #2^4) &)}] Show[ContourPlot[x + y, {x, y} \[Element] Disk[]], Graphics[{Red, PointSize[Large], Point[{x, y} /. Last[%]]}]] which would suggest it can take this as a penalty for points not satisfying the constraint? $\endgroup$ – chris Oct 29 '14 at 11:24
  • 1
    $\begingroup$ @chris Try ac = {}; NMinimize[{x^2 - y^2, x x + y y <= 1}, {x, y}, Method -> {"SimulatedAnnealing", "PenaltyFunction" -> ((AppendTo[ac, {##}]; 1) &)}] ; Print@ac $\endgroup$ – Dr. belisarius Oct 29 '14 at 12:26
  • 2
    $\begingroup$ Related, but modestly helpful: library.wolfram.com/infocenter/MathSource/7623/…. I think the penalty function is applied (added to the objective function) only when the point goes outside the constraint. I don't know how Nelder-Mead works, but I suspect it enforces constraints in a way that is clobbered by using your kinds of "PenaltyFunction". For the unit disk, try "PenaltyFunction" -> (2(Norm[{##}] - 1) &). $\endgroup$ – Michael E2 Oct 29 '14 at 16:24
  • 1
    $\begingroup$ To make it thrash, give it a negative penalty: "PenaltyFunction" -> (-(Norm[{##}] - 1) &) $\endgroup$ – Michael E2 Oct 29 '14 at 16:26
  • 1
    $\begingroup$ Great question! I always wondered about this, in the end giving up and turning to C++... $\endgroup$ – dr.blochwave Oct 29 '14 at 21:42
9
$\begingroup$

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}}  *)
$\endgroup$
1
  • 1
    $\begingroup$ Interesting. Given the solution is posted in 2004, I tested the code in v5.0, the correct result is given, and the corresponding penalty function is: i.stack.imgur.com/sWCQz.png Seems that z has been eliminated using the constraints in some way, this might explain why the original "PenaltyFunction works in v5.0. $\endgroup$ – xzczd Dec 28 '20 at 6:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.