You might get better results by using LogicalExpand on the result of Reduce on the constraints and taking the best result.

In[1]:= n = 100;
Y = 1000000;

In[15]:= r1 = 
  Reduce[{p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && 
     p2 >= 0 && ap10 >= 0 && 
     p100 >= 0 && (1/100)*Y*(-10*ap10 + 30*p1 + 40*p2 + 100*p100) == 
      1000}, {p0, p1, p2, p100, ap10}];

In[16]:= Quiet[
 NMaximize[{n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2 , #}, {p0, p1, p2, 
     p100, ap10}, Method -> "DifferentialEvolution"] & /@ 
  LogicalExpand[r1]]

Out[16]= {0., {p0 -> 0.996973, p1 -> 0.00289558, p2 -> 0, 
   p100 -> 0.000131325, ap10 -> 0}} || {0., {p0 -> 0.998052, 
   p1 -> 0.00135491, p2 -> 0, p100 -> 0.000593527, 
   ap10 -> 0}} || {0., {p0 -> 0.630183, p1 -> 0.0949542, p2 -> 0, 
   p100 -> 0, ap10 -> 0.274863}} || {2.87524*10^-10, {p0 -> 0.997462, 
   p1 -> 0.0011485, p2 -> 0.000577472, p100 -> 0.000459685, 
   ap10 -> 0.000352236}} || {3.76256*10^-8, {p0 -> 0.996667, 
   p1 -> 0.000555751, p2 -> 0.000648745, p100 -> 0.000715146, 
   ap10 -> 0.00141369}} || {2.87269*10^-9, {p0 -> 0.9975, 
   p1 -> 0.000357294, p2 -> 0.000416864, p100 -> 0.000816955, 
   ap10 -> 0.000908887}} || {0.0251849, {p0 -> 0.940801, 
   p1 -> 0.00797022, p2 -> 0.00746356, p100 -> 1.02811*10^-12, 
   ap10 -> 0.0437649}}

If I symbolically solve the Kuhn-Tucker equations for the problem with the second constraint = 1000, I get 51 solutions, the first of which has an objective function value of 0

In[1]:= KTEqs[obj_, cons_List, vars_] :=
 Module[{consconvrule = {GreaterEqual[x_, y_] -> LessEqual[y - x, 0], 
     Equal[x_, y_] -> Equal[x - y, 0], 
     LessEqual[x_, y_] -> LessEqual[x - y, 0],
     LessEqual[lb_, x_, ub_] -> LessEqual[(x - lb) (x - ub), 0], 
     GreaterEqual[ub_, x_, lb_] -> LessEqual[(x - lb) (x - ub), 0]} 
    , stdcons, eqcons, ineqcons, lambdas, mus, lagrangian, eqs1, eqs2,
    eqs3, alleqns, allvars },
  stdcons = cons /. consconvrule;
  eqcons = Cases[stdcons, Equal[_, 0]][[All, 1]];
  ineqcons = Cases[stdcons, LessEqual[_, 0]][[All, 1]];
  lambdas = Array[\[Lambda], Length[eqcons]];
  mus = Array[\[Mu], Length[ineqcons]];
  lagrangian = obj + lambdas.eqcons + mus.ineqcons;
  eqs1 = Thread[ D[lagrangian, {vars}] == 0];
  eqs2 = Thread[mus >=   0];
  eqs3 = Thread[mus*ineqcons == 0];
  alleqns = Join[eqs1, eqs2, eqs3, stdcons];
  allvars = Join[vars, lambdas, mus];
  {alleqns, allvars}
  ]

In[2]:= n = 100;
Y = 1000000;

In[4]:= kts1 = 
  KTEqs[-n!/((94)!*4!)*p0^(94)*ap10^4*p1*
    p2, {p0 + p1 + p100 + p2 + ap10 == 
     1, (1/100)*Y*(-10*ap10 + 30*p1 + 40*p2 + 100*p100) == 1000, 
    p0 >= 0, p1 >= 0, p2 >= 0, p100 >= 0, ap10 >= 0}, {p0, p1, p2, 
    p100, ap10}];

In[5]:= res1 = 
  Reduce[Sequence @@ kts1, Reals, Backsubstitution -> True];

In[6]:= Length[res1]

Out[6]= 51

In[7]:= res2 = 
  Reduce[{obj == -n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, res1}, {obj, p0,
     p1, p2, p100, ap10}];

In[8]:= Length[res2]

Out[8]= 51

In[9]:= res2[[1]]

Out[9]= \[Lambda][1] == 0 && \[Lambda][2] == 0 && \[Mu][1] == 
  0 && \[Mu][2] == 0 && \[Mu][3] == 0 && \[Mu][4] == 0 && \[Mu][5] == 
  0 && obj == 0 && p0 == 0 && p1 == 101/800 && p2 == 101/1000 && 
 p100 == 0 && ap10 == 3091/4000

In[10]:= N[res2[[1]]]

Out[10]= \[Lambda][1.] == 0. && \[Lambda][2.] == 0. && \[Mu][1.] == 
  0. && \[Mu][2.] == 0. && \[Mu][3.] == 0. && \[Mu][4.] == 
  0. && \[Mu][5.] == 0. && obj == 0. && p0 == 0. && p1 == 0.12625 && 
 p2 == 0.101 && p100 == 0. && ap10 == 0.77275

I wrote up solving the Kuhn-Tucker equations with Reduce some years ago in the Mathematica Journal.

http://www.mathematica-journal.com/issue/v9i4/contents/Tricks9-4/Tricks9-4_2.html

There are some typos in the article but you can use it for a reference.

Thank you for your answer. I've made some experiments and found that it gives good results with adding scalingfactor->1. But if equations becomes more complicated it gives strange results again. So what does scalefactor do and which paramenter should I use or maybe I should add some more methodes? According to iMin with global it gives better results giving the zero but the right answer is 0.025. It would be amazing if NMaximize function will be the single method of finding this naximum. But even NMaximize cant cope with more complicated equations((

n = 100; 
Y = 1000000; 
st = 1000;
o1 = 1;
o2 = 1;

q = Table[
   NMaximize[{n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, 
     p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && 
      p2 >= 0 && ap10 >= 0 && 
      p100 >=  
       0 && (1/100)*Y*(-10*ap10 + 30*p1 + 40*p2  + 100*p100) == 
       st*o}, {p0, p1, p2, p100, ap10}, 
    Method -> {"DifferentialEvolution", "ScalingFactor" -> 1}], {o, 
    o1, o2}];
{{0.0251849, {p0 -> 0.940801, p1 -> 0.00797022, p2 -> 0.00746355, 
   p100 -> 4.39477*10^-9, ap10 -> 0.0437649}}}

If you dont mind please take a look at a bit more complicated example.

Attachments:

Maximums.nb