# Why cannot Minimize find a solution?

I am trying to minimize the following function. Mathematica gives out an error saying: NMinimize was unable to generate any initial points satisfying the inequality constraints. Mathematica eventually runs out of memory without finding any solution.

In this case, the values (a1=87, a2=8, a3=4, a4=1) should work (I don't know if it is THE solution).

f[a1_?NumericQ, a2_?NumericQ, a3_?NumericQ, a4_?NumericQ] :=
Integrate[
1 - (1 - E^(-0.79788*34/3140 *r))^a1*(1 - E^(-0.79788*27/3140 *r))^
a2*(1 - E^(-0.79788*20/3140 *r))^a3*(1 - E^(-0.79788*14/3140 *r))^
a4, {r, 475, 1200}];
Minimize[{f[a1, a2, a3, a4], {a1 >= 0,
a2 >= 0, a3 >= 0, a4 >= 0, a1 + a2 + a3 + a4 == 100,
1 - (1 - E^(-0.79788*34/3140 *475))^
a1*(1 - E^(-0.79788*27/3140 *475))^
a2*(1 - E^(-0.79788*20/3140 *475))^
a3*(1 - E^(-0.79788*14/3140 *475))^a4 >= 0.9}}, {a1, a2, a3,
a4}, Integers]


Any help is much appreciated!

The full error is: NMinimize::incst: NMinimize was unable to generate any initial points satisfying the \ inequality constraints {-0.1+0.815439^Round[a4] 0.91054^Round[a3] \ <<19>>^Round[a2] \ 0.983489^Round[100-Round[<<1>>]-Round[<<1>>]-Round[<<1>>]]<=0,-Round[<\ <1>>]<=0}. The initial region specified may not contain any feasible \ points. Changing the initial region or specifying explicit initial \ points may provide a better solution

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Commented Aug 7, 2016 at 2:37
• Replacing decimal numbers by rational numbers eliminates the error message but does not produce an answer. Also, I recommend replacing Integrate by NIntegrate. Commented Aug 7, 2016 at 4:32
• Thank you for the suggestion. If I replace 0.9 with 9/10 and 0.79788 with sqrt(2/pi), the error disappears but the output is now the same as the input. Any ideas why this might be happening? Commented Aug 7, 2016 at 4:44
• Replacing Minimize by NMinimize does a little better, but not much. Note that f[100, 0, 0, 0] is much smaller than f[87, 8, 4, 1]. I would guess that {100,0,0,0} is the minimum, but I have not proven it. Commented Aug 7, 2016 at 4:50
• With NMinimize, I again get the NMinimize:incst error. You are right that 'f[100, 0, 0, 0]' would be the minimum, but it does not satisfy the last constraint (the one with >= 0.9), whereas 'f[87, 8, 4, 1]' does. Commented Aug 7, 2016 at 4:55

f[a1_?NumericQ, a2_?NumericQ, a3_?NumericQ, a4_?NumericQ] :=
NIntegrate[
1 - (1 - E^(-79788/100000*34/3140*r))^a1*(1 - E^(-79788/100000*27/3140*r))^a2*
(1 - E^(-79788/100000*20/3140*r))^a3*(1 - E^(-79788/100000*14/3140*r))^a4,
{r, 475, 1200}];
NMinimize[{f[a1, a2, a3, a4],
a1 >= 0 && a2 >= 0 && a3 >= 0 && a4 >= 0 && a1 + a2 + a3 + a4 == 100 &&
1 - (1 - E^(-79788/100000*34/3140*475))^a1*(1 - E^(-79788/100000*27/3140*475))^a2*
(1 - E^(-79788/100000*20/3140*475))^a3*(1 - E^(-79788/100000*14/3140*475))^a4 >= 9/10},
{{a1, 70, 80}, {a2, 20, 30}, {a3, 0, 1}, {a4, 0, 1}}]

(* {186.059, {a1 -> 71.7089, a2 -> 28.2911, a3 -> 3.93472*10^-8, a4 -> 5.32211*10^-9}} *)


The nearest set of integers that satisfies the constraints is {71, 29, 0, 0}, and

f[71, 29, 0, 0]
(* 187.215 *)


is, incidentally, much smaller than

f[87, 8, 4, 1]
(* 210.032 *)


This is not a rigorous proof, because NMinimize finds only local minima for nonlinear problems. Nonetheless, the solution obtained is plausible: It sets a3 and a4 to zero and makes a2 as small as possible subject to the last constraint. ({100, 0, 0, 0} would be the minimum, were it not for the last constraint.)

• Thank you! Even without the bounds, I do get a solution in ~77 seconds now. Commented Aug 7, 2016 at 6:02

Let's attack this problem a different way.

IntegerPartitions is a convenient way of finding all possible ways that four numbers sum to 100. But IntegerPartitions is based on numbers greater than zero and you are including the possibility of some of your numbers being zero. So we need a trick. Let's look at four numbers that sum to 10.

IntegerPartitions[10 + 4, {4}] - 1


which gives us

{{10,0,0,0}, {9,1,0,0}, {8,2,0,0}, {8,1,1,0}, {7,3,0,0}, {7,2,1,0},
{7,1,1,1}, {6,4,0,0}, {6,3,1,0}, {6,2,2,0}, {6,2,1,1}, {5,5,0,0},
{5,4,1,0}, {5,3,2,0}, {5,3,1,1}, {5,2,2,1}, {4,4,2,0}, {4,4,1,1},
{4,3,3,0}, {4,3,2,1}, {4,2,2,2}, {3,3,3,1}, {3,3,2,2}}


and that gives you four numbers >=0 and which sum to 10. That doesn't take into account possible permutations of each of those sets, but we will deal with that in a moment.

Next some utility functions.

f[{a1_, a2_, a3_, a4_}] :=
If[1-(1-E^(-0.79788*34/3140*475))^a1*(1-E^(-0.79788*27/3140*475))^a2*
(1-E^(-0.79788*20/3140*475))^a3*(1-E^(-0.79788*14/3140*475))^a4>=0.9,
NIntegrate[
1-(1-E^(-0.79788*34/3140*r))^a1*(1-E^(-0.79788*27/3140*r))^a2*
(1-E^(-0.79788*20/3140*r))^a3*(1-E^(-0.79788*14/3140*r))^a4,
{r, 475, 1200}],
Infinity];
g[v_] := Map[{f[# - 1], # - 1} &, Permutations[v]];


How big is this problem going to be?

Length[IntegerPartitions[100 + 4, {4}] - 1]


gives you 8037, not counting permutations, so this may take a while. So let's sneak up on the problem by looking at smaller problems where the numbers sum to 20, to 30, to 40 and to 50 where we can quickly get an answer and see how the problem is going to scale as we look towards your 100.

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[20 + 4, {4}]], 1]]]]


gives you {3.4375, {298.141, {0, 0, 16, 4}}}

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[30 + 4, {4}]], 1]]]]


gives you {28.5156, {260.326, {0, 9, 21, 0}}}

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[40 + 4, {4}]], 1]]]]


gives you {88.5781, {242.616, {1, 25, 14, 0}}}

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[50 + 4, {4}]], 1]]]]


gives you {191.063, {225.491, {1, 42, 7, 0}}}

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[60 + 4, {4}]], 1]]]]


gives you {334.719, {206.077, {2, 58, 0, 0}}}

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[70 + 4, {4}]], 1]]]]


gives you {539.109, {201.847, {19, 51, 0, 0}}}

So all this gives you what looks like enough information to extrapolate the amount of time it might take to give you a solution for 100. It doesn't look like it will take as long as I had feared it might.

This incorporates the constraint about >=0.9.

And it appears that the solution is

Timing[First[Sort[Flatten[Map[g, IntegerPartitions[100 + 4, {4}]], 1]]]]


which gives you {1520.27, {187.215, {71, 29, 0, 0}}}

That matches bbgodfrey's solution of looking for minima while ignoring the constraint that a1,a2,a3,a4 be integers and then looking for the nearest point which does have integers.

• Thank you for the detailed explanation Bill. You are right that the minimum would be for f[100, 0, 0, 0] but it does not satisfy my last constraint (the one with the >=0.9). For f[87, 8, 4, 1] it does satisfy. Does the minimize operation have to extrapolate every possible combination? I thought it would be able to optimize and return a solution without having to enumerate all possible combinations. Commented Aug 7, 2016 at 5:19
• I realize everyone wants a fast and easy solution to their problem. But your requirement that a1,a2,a3,a4 be integers and you have an assortment of floating point parameters and you have a combination of exponential expressions, all those might make it difficult to find a simple analytic solution to your problem. If you can show a simple way to conclude what the minimum is without looking at all the cases then I'd be interested in seeing it.
– Bill
Commented Aug 7, 2016 at 5:29
• Thank you again! It is true that I was foolishly expecting a faster answer for a difficult question. Commented Aug 7, 2016 at 5:57
• Just out of curiosity, may I ask what are the basic configurations of your machine? I am able to get results for different values of r for a1+a2+a3+a4=100, in the range of 75 - 100 seconds. I am running an i7 with 8GB RAM. Commented Aug 7, 2016 at 6:06
• IF you could justify that the "surface" in 4 variables is "nice" then you might be able to do hill climbing across that surface until you find the minimum you are looking for. I certainly don't have the intuition to be able to see what that surface is going to look like, particularly with it hidden inside your integration. And in the time I spent trying to see if I could find an explicit solution for the result of the integration without having substituted in values for a1,a2,a3,a4 I was not able to do this.
– Bill
Commented Aug 7, 2016 at 6:07