# 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

## 1 Answer

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


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[{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
{"ft_lang_label":"__label__en","ft_lang_prob":0.6951083,"math_prob":0.98979396,"size":3402,"snap":"2021-21-2021-25","text_gpt3_token_len":1084,"char_repetition_ratio":0.1433196,"word_repetition_ratio":0.25,"special_character_ratio":0.3659612,"punctuation_ratio":0.18167701,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9973019,"pos_list":[0,1,2,3,4,5,6],"im_url_duplicate_count":[null,5,null,5,null,5,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2021-06-15T16:52:10Z\",\"WARC-Record-ID\":\"<urn:uuid:b4f9b1c9-d62e-4b27-a232-695564ffec4d>\",\"Content-Length\":\"178529\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:0bc2a3ff-5284-4e50-a8ff-52dcd33eaf62>\",\"WARC-Concurrent-To\":\"<urn:uuid:d144f654-e52a-4d2a-b390-f697bd0ff0ab>\",\"WARC-IP-Address\":\"151.101.129.69\",\"WARC-Target-URI\":\"https://mathematica.stackexchange.com/questions/64388/penaltyfunction-option-in-nminimize\",\"WARC-Payload-Digest\":\"sha1:THWDZ7KJ5T2HW432JI5USSJHAL3MTZB4\",\"WARC-Block-Digest\":\"sha1:QKXZ4JV2PFVJDSZMZKBRV6DP345KRKAV\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2021/CC-MAIN-2021-25/CC-MAIN-2021-25_segments_1623487621450.29_warc_CC-MAIN-20210615145601-20210615175601-00553.warc.gz\"}"}