# PenaltyFunction Option in NMinimize

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.

• 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? – chris Oct 29 '14 at 11:24
• @chris Try ac = {}; NMinimize[{x^2 - y^2, x x + y y <= 1}, {x, y}, Method -> {"SimulatedAnnealing", "PenaltyFunction" -> ((AppendTo[ac, {##}]; 1) &)}] ; Print@ac – Dr. belisarius Oct 29 '14 at 12:26
• 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) &). – Michael E2 Oct 29 '14 at 16:24
• To make it thrash, give it a negative penalty: "PenaltyFunction" -> (-(Norm[{##}] - 1) &) – Michael E2 Oct 29 '14 at 16:26
• Great question! I always wondered about this, in the end giving up and turning to C++... – dr.blochwave Oct 29 '14 at 21:42

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 OptimizationNMinimizeDumppenalty (the current step is given by OptimizationNMinimizeDumpitr):

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, OptimizationNMinimizeDumpitr = itr},
Return[OptimizationNMinimizeDumppenalty, 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, OptimizationNMinimizeDumpitr = itr},
Return[OptimizationNMinimizeDumppenalty, 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, OptimizationNMinimizeDumpitr = itr},
Return[OptimizationNMinimizeDumppenalty, 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
*)


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[{OptimizationNMinimizeDumpitr = itr},
Return[OptimizationNMinimizeDumppenalty, 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}}  *)

• 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. – xzczd Dec 28 '20 at 6:15