Message Boards Message Boards

0
|
8572 Views
|
34 Replies
|
7 Total Likes
View groups...
Share
Share this post:

Maximize a function with multinomial conditions?

Posted 6 years ago

Consider the following code:

Maximize[{1000000/100*(0*a + 25*b + 75*c + 100*d),  Reduce[a + b + c + d == 1 && a >= 0 && b >= 0 && c >= 0 && d >= 0 &&
     a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 50*100*a^99*b^2 >= 0.05, {a, b, c, d}, Reals]},
     {a, b, c, d}]

Hi. I'm trying to maximize this excpresssion but after several hours it is still runing. Without "reduce" it give me the wrong answer. i need to get the output somthing like this a=0.9538 b=0.0045 c=0.0280 d=0.0137. Maybe I need to use kinda gradient seearch technique?

POSTED BY: Alex Graham
34 Replies

the random starting points are generated using the RandomReal function

POSTED BY: Frank Kampas
Posted 5 years ago

Thank you for your reply. Correct me please if I'm wrong. So for each random search it sets up the random starting points to all p and searches the minimum up to the 1 for all p variables and objective function. But I have another question, what is the seed in this case? Is it an extra starting point from 0 to 1 for all p? Why is it set up on zero? Regards, Alex

POSTED BY: Alex Graham

check to see if both solutions satisfy the constraints. if so, use the one with the better result for the objective function value

POSTED BY: Frank Kampas
Posted 5 years ago

1)I did the calculation with multiple random searches, using ParametricIPOPTMinimize, using the starting values as the parameters. I got the result I posted for 10, 100, and 1000 random searches. 2)Thank you for your reply. I have one specializing question. [{{p0, p1, p2, p3, p100}, 0, 1}], 10, 0] 0,1 means it is searching within the 0:1 interval, 10 means number of serches. What does 0] mean? 3)the final input, 0, is the seed for the 10 random searches.


Seems some time has passed away since I was on the thread. I've got what I wanted, but still have some questions about minimization process itself. If I get it right, the seed is the point for all p0,p1,p2,p3,p100. And it finds the minimum checking the function from the 0 to 1 for all p. So my question are 1) Do the algorithm start to find minimum for all p0=p1=p2=p3=p4=p100=0 at the same time or take only p0=0 and finds some minimum, then take only p1=0 finds some minimum and so on? 2) If the seed is the same for all random searches, so searches are the same or not? If they are not the same what is the difference between them? Why we need several of them but not only one? Thanks in advance.

POSTED BY: Alex Graham
Posted 6 years ago

Hello again! I modified my part of the code, extended its possibilities and make it convinient. Now it generates all necessary things. I noticed when I use more than 4 error the results from imin and FindMaximum are a bit different. For example for some six errors I have these results. Starting points are automatically transferred from imim to FindMaximum

{135.292, {-617734., {Subscript[p, 0] -> 0.879188, 
   Subscript[p, 1] -> 0., Subscript[p, 2] -> 0.051358, 
   Subscript[p, 3] -> 0.00122958, Subscript[p, 4] -> 0., 
   Subscript[p, 5] -> 0.0205016, Subscript[p, 6] -> 0., 
   Subscript[p, 100] -> 0.0477225}}}

{616555., {Subscript[p, 0] -> 0.874134, Subscript[p, 1] -> 0., 
  Subscript[p, 2] -> 0.041653, Subscript[p, 3] -> 0.00561729, 
  Subscript[p, 4] -> 0.0301089, Subscript[p, 5] -> 0., 
  Subscript[p, 6] -> 0., Subscript[p, 100] -> 0.048487}}

Then which result should I choose, iMin or FindMaximum?

POSTED BY: Alex Graham

rerunning with several different random seeds gave the same result

POSTED BY: Frank Kampas
Posted 6 years ago

I compared my results with theirs. For 1,2,3,4 errors with n=101 they got 29.3 and n=10 259.9 With my set of outcomes and the code n=10

{-260546., {p0 -> 0.642418, p1 -> 0., p2 -> 0.0769742, 
  p3 -> 0.0222692, p4 -> 0., p100 -> 0.258339}}

with n 101

{-29404., {p0 -> 0.960175, p1 -> 0., p2 -> 0.00805374, 
  p3 -> 0.00260652, p4 -> 0., p100 -> 0.0291647}}

So it definately works as it should be, perfectly or with very low inaccuracy. Or they used modified set of outcomes... Thanks for your replies!

POSTED BY: Alex Graham
Posted 6 years ago

Hello again. I've made a program to generate the set of outcomes depending on errors. So I'm sure that 99% that it's correct. So I've tried to find probabilities for 4 mistakes. And got these results. This means that the probabilities for having a error with d1,d2,d4 are zero, that it's not ok in theory or not(i'm not sure). What do you think on it? Anyway, thanks for all.

POSTED BY: Alex Graham

iMin didn't run because t is a List

In[52]:= iMin[-Y/100*(0*p0 + d1*p1 + d2*p2 + d3*p3 + 100*p100), 
 List @@ (p0 + p1 + p2 + p3 + p100 == 1 && p0 >= 0 && p1 >= 0 && 
    p2 >= 0 && p3 >= 0 && p100 >= 0 && (Sequence @@ t) == 1/20), 
 Thread[{{p0, p1, p2, p3, p100}, 0, 1}], 10, 0]

During evaluation of In[52]:= {{Infeasible_Problem_Detected,8},{Solve_Succeeded,2}}

Out[52]= {-539209., {p0 -> 0.925934, p1 -> 0.0127001, p2 -> 0.0099051,
   p3 -> 0.000352082, p100 -> 0.0511088}}
POSTED BY: Frank Kampas
Posted 6 years ago

Seems my modified code dosent work with FindMaximum and Imin without reduction, because the equation output is turning out with brakets.

POSTED BY: Alex Graham

the final input, 0, is the seed for the 10 random searches.
I don't understand your problem well enough to answer your second question However, Mathematica is very powerful so I would think that the answer is yes.

POSTED BY: Frank Kampas
Posted 6 years ago

Thank you for the explanation.

POSTED BY: Alex Graham
Posted 6 years ago

Thank you for your reply. I have one specializing question. [{{p0, p1, p2, p3, p100}, 0, 1}], 10, 0] 0,1 means it is searching within the 0:1 interval, 10 means number of serches. What does 0] mean?

POSTED BY: Alex Graham

using the code at http://community.wolfram.com/groups/-/m/t/1164680 and changing the sign of the objective function, since the code does minimization:

In[11]:= n = 56;
d1 = 10;
d2 = 15;
d3 = 16;
Y = 1000000;

In[18]:= iMin[-Y/100*(0*p0 + d1*p1 + d2*p2 + d3*p3 + 100*p100), 
 List @@ (p0 + p1 + p2 + p3 + p100 == 1 && p0 >= 0 && p1 >= 0 && 
    p2 >= 0 && p3 >= 0 && p100 >= 0 && 
    p0^n + n*p0^(n - 1)*p1 + n*p0^(n - 1)*p2 + n*p0^(n - 1)*p3 + 
      n*(n - 1)/2*p0^(n - 2)*p1^2 + n*(n - 1)*p0^(n - 2)*p1*p2 + 
      n*(n - 1)/2*p0^(n - 2)*p2^2 + n*(n - 1) p0^(n - 2)*p1*p3 + 
      n (n - 1) p0^(n - 2)*p2*p3 + 
      n (n - 1) (n - 2)/6*p0^(n - 3) p1^3 + 
      n (n - 1) (n - 2)/2*p0^(n - 3)*p1^2*p2 + 
      n (n - 1) (n - 2)/2*p0^(n - 3)*p1*p2^2 + 
      n (n - 1) (n - 2)/2*p0^(n - 3)*p1^2*p3 + 
      n (n - 1) (n - 2) p0^(n - 3)*p1*p2*p3
          == 1/20), Thread[{{p0, p1, p2, p3, p100}, 0, 1}], 10, 0]

During evaluation of In[18]:= {{Infeasible_Problem_Detected,8},{Solve_Succeeded,2}}

Out[18]= {-53920.9, {p0 -> 0.925934, p1 -> 0.0127001, p2 -> 0.0099051,
   p3 -> 0.000352082, p100 -> 0.0511088}}
POSTED BY: Frank Kampas
Posted 6 years ago

Could you please show me how to get the starting points with ParametricIPOPTMinimize for FindMaximum function?. I've of course tried to do it on my own, read your promts on http://community.wolfram.com/groups/-/m/t/1164680 ..., manuals but unfortunately I havent succeed.

POSTED BY: Alex Graham
In[3]:= n = 56;
d1 = 10;
d2 = 15;
d3 = 16;
Y = 1000000;

In[13]:= FindMaximum[{Y/100*(0*p0 + d1*p1 + d2*p2 + d3*p3 + 100*p100),
   p0 + p1 + p2 + p3 + p100 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 && 
   p3 >= 0 && p100 >= 0 && 
   p0^n + n*p0^(n - 1)*p1 + n*p0^(n - 1)*p2 + n*p0^(n - 1)*p3 + 
     n*(n - 1)/2*p0^(n - 2)*p1^2 + n*(n - 1)*p0^(n - 2)*p1*p2 + 
     n*(n - 1)/2*p0^(n - 2)*p2^2 + n*(n - 1) p0^(n - 2)*p1*p3 + 
     n (n - 1) p0^(n - 2)*p2*p3 + 
     n (n - 1) (n - 2)/6*p0^(n - 3) p1^3 + 
     n (n - 1) (n - 2)/2*p0^(n - 3)*p1^2*p2 + 
     n (n - 1) (n - 2)/2*p0^(n - 3)*p1*p2^2 + 
     n (n - 1) (n - 2)/2*p0^(n - 3)*p1^2*p3 + 
     n (n - 1) (n - 2) p0^(n - 3)*p1*p2*p3
         == 1/20}, {{p0, .925}, {p1, .013}, {p2, .01}, {p3, 
   0}, {p100, .05}}]

Out[13]= {53920.9, {p0 -> 0.925934, p1 -> 0.0127001, p2 -> 0.0099051, 
  p3 -> 0.000352085, p100 -> 0.0511088}}
POSTED BY: Frank Kampas

I did the calculation with multiple random searches, using ParametricIPOPTMinimize, using the starting values as the parameters. I got the result I posted for 10, 100, and 1000 random searches.

POSTED BY: Frank Kampas
Posted 5 years ago

1)I did the calculation with multiple random searches, using ParametricIPOPTMinimize, using the starting values as the parameters. I got the result I posted for 10, 100, and 1000 random searches. 2)Thank you for your reply. I have one specializing question. [{{p0, p1, p2, p3, p100}, 0, 1}], 10, 0] 0,1 means it is searching within the 0:1 interval, 10 means number of serches. What does 0] mean? 3)the final input, 0, is the seed for the 10 random searches. Seems some time has passed away since I was on the thread. I've got what I wanted, but still have some questions about minimization process itself. If I get it right, the seed is the point for all p0,p1,p2,p3,p100. And it finds the minimum checking the function from the 0 to 1 for all p. So my question are 1) Do the algorithm start to find minimum for all p0=p1=p2=p3=p4=p100=0 at the same time or take only p0=0 and finds some minimum, then take only p1=0 finds some minimum and so on? 2) If the seed the same for all random searches, so searches are the same or not? If they are not the same what is the difference between them? Why we need several of them but not only one? Thanks in advance.

POSTED BY: Alex Graham
Posted 6 years ago

it definately makes sense as it is significantly tighter than stringer bound 64490 as it should be.

POSTED BY: Alex Graham

I tried doing the problem with some experimental code I've been writing and got

{53920.9, {p0 -> 0.925934, p1 -> 0.0127001, p2 -> 0.00990509, 
  p3 -> 0.000352097, p100 -> 0.0511088}

Does that solution make sense?

POSTED BY: Frank Kampas

why are there quotation marks around the 4 and 11? That's probably why Reduce doesn't work.

POSTED BY: Frank Kampas
Posted 6 years ago

Year. sorry this was my mistake. I dont think that these comments will brake reduce function. Now, I have NMaximizing function dosent work for a large n.

POSTED BY: Alex Graham
Posted 6 years ago

Hello again. Suppose, now a have three mistakes in the sample and the system becomes more complicated. And the problem is that the Reduce function failures to help.(((

In[4]:= 
  n = 50;
  d1 = 10;
  d2 = 15;
  d3 = 16;
  Y = 1000000;
  r = Reduce[
     p0 + p1 + p2 + p3 + p100 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 &&
       p3 >= 0 && p100 >= 0 &&
      p0^n + n*p0^(n - 1)*p1 + n*p0^(n - 1)*p2 + n*p0^(n - 1)*p3 "4" +
         n*(n - 1)/2*p0^(n - 2)*p1^2 + n*(n - 1)*p0^(n - 2)*p1*p2 + 
        n*(n - 1)/2*p0^(n - 2)*p2^2 +
        n*(n - 1) p0^(n - 2)*p1*p3 "8" + n (n - 1) p0^(n - 2)*p2*p3 + 
        n (n - 1) (n - 2)/6*p0^(n - 3) p1^3 + 
        n (n - 1) (n - 2)/2*p0^(n - 3)*p1^2*p2 "11" +
        n (n - 1) (n - 2)/2*p0^(n - 3)*p1*p2^2 + 
        n (n - 1) (n - 2)/2*p0^(n - 3)*p1^2*p3 + 
        n (n - 1) (n - 2) p0^(n - 3)*p1*p2*p3
       == 1/20, {p0, p1, p2, p3, p100}, Reals, 
     Backsubstitution -> True];
  q = NMaximize[{Y/100*(0*p0 + d1*p1 + d2*p2 + d3*p3 + 100*p100), 
     r}, {p0, p1, p2, p3, p100}]

  },
 {\[Placeholder]}
}

During evaluation of In[4]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce.

During evaluation of In[4]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce.

During evaluation of In[4]:= Reduce::ivar: 0.6565924971797801` is not a valid variable.

During evaluation of In[4]:= NMaximize::bcons: The following constraints are not valid: {Reduce[{p0+p1+p100+p2+p3==1,p0>=0,p1>=0,p2>=0,p3>=0,p100>=0,p0^50+50 p0^49 p1+1225 p0^48 p1^2+19600 p0^47 p1^3+50 p0^49 p2+2450 p0^48 p1 p2+58800 11 p0^47 p1^2 p2+1225 p0^48 p2^2+58800 p0^47 p1 p2^2+50 4 p0^49 p3+2450 8 p0^48 p1 p3+58800 p0^47 p1^2 p3+2450 p0^48 p2 p3+117600 p0^47 p1 p2 p3==1/20},{p0,p1,p2,p3,p100},\[DoubleStruckCapitalR],Backsubstitution->True]}. Constraints should be equalities, inequalities, or domain specifications involving the variables.

Out[4]= {{"n=101;\[IndentingNewLine]d1=25;\[IndentingNewLine]d2=75;\
\[IndentingNewLine]s=0;\[IndentingNewLine]Y=1000000;\
\[IndentingNewLine]r=Reduce[a+b+c+d\[Equal]1&& a\[GreaterEqual]0&&b\
\[GreaterEqual]0&&c\[GreaterEqual]0&&d\[GreaterEqual]0\
&&\[IndentingNewLine]a^n+n*a^(n-1)*b+n*a^(n-1)*c+n*(n-1)*a^(n-2)*b*c+\
n*(n-1)/2*a^(n-2)*b^2==1/20,{a,b,c,d},Reals,Backsubstitution\[Rule]\
True];\[IndentingNewLine]q=NMaximize[{Y/100*(0*a+d1*b+d2*c+100*d),r},{\
a,b,c,d}];\[IndentingNewLine]l=q+s" Null^6 NMaximize[{10000 (10 p1 + 
        100 p100 + 15 p2 + 16 p3), 
     Reduce[p0 + p1 + p100 + p2 + p3 == 1 && p0 >= 0 && p1 >= 0 && 
       p2 >= 0 && p3 >= 0 && p100 >= 0 && 
       p0^50 + 50 p0^49 p1 + 1225 p0^48 p1^2 + 19600 p0^47 p1^3 + 
         50 p0^49 p2 + 2450 p0^48 p1 p2 + 58800 "11" p0^47 p1^2 p2 + 
         1225 p0^48 p2^2 + 58800 p0^47 p1 p2^2 + 50 "4" p0^49 p3 + 
         2450 "8" p0^48 p1 p3 + 58800 p0^47 p1^2 p3 + 
         2450 p0^48 p2 p3 + 117600 p0^47 p1 p2 p3 == 1/20, {p0, p1, 
       p2, p3, p100}, Reals, Backsubstitution -> True]}, {p0, p1, p2, 
     p3, p100}]}, {\[Placeholder]}}
POSTED BY: Alex Graham

I defined a function that could help;

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, cons];
  allvars = Join[vars, lambdas, mus];
  {alleqns, allvars}
  ]

KTMinimize[obj_, cons_List, vars_List] :=
 Block[{kteqs, r, rls, objvals, minobj, objrls, res},
  kteqs = KTEqs[obj, cons, vars];
  r = LogicalExpand @  
    Reduce[Sequence @@ kteqs, Backsubstitution -> True, 
     Cubics -> True, Quartics -> True];

  If[Head[r] === And,
   rls = List @@ ToRules[r];
   objvals = obj /. rls;
   res = {objvals, rls},
   (* Else *)
   rls = List @@ (ToRules /@ r);
   objvals = obj /. rls;
   minobj = Min[objvals];
   objrls = Thread[{objvals, rls}];
   res = Select[objrls, #[[1]] == minobj &];
   If[Length[res] == 1, res = res[[1]]];
   ];

  res
  ]

In[3]:= AbsoluteTiming[
 KTMinimize[-1000000/100*(25*b + 75*c + 100*(1 - a - b - c)), {a >= 0,
     b >= 0, c >= 0, 
    a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 
      50*100*a^99*b^2 >= 1/20}, {a, b, c}] // N]

Out[3]= {3.48344, {-35843.9, {a -> 0.953819, b -> 0.00444532, 
   c -> 0.028011, \[Mu][1.] -> 
    3.42892*10^-10, \[Mu][2.] -> -5.08376*10^-10, \[Mu][3.] -> 
    5.66923*10^-10, \[Mu][4.] -> 190922.}}}
POSTED BY: Frank Kampas

to troubleshoot your code, you need to separate the various evaluations to see where the problem is occurring.

POSTED BY: Frank Kampas
Posted 6 years ago

Thanks for your reply. Just one techical quistion, I think this should be last question for now. Is it possible to get more cleaner output without nulls and placeholder?

In[9]:= {
 {n = 101;
  d1 = 25;
  d2 = 75;
  s = 20000;
  Y = 1000000;
  r = Reduce[
     a + b + c + d == 1 && a >= 0 && b >= 0 && c >= 0 && d >= 0 &&
      a^n + n*a^(n - 1)*b + n*a^(n - 1)*c + n*(n - 1)*a^(n - 2)*b*c + 
        n*(n - 1)/2*a^(n - 2)*b^2 == 1/20, {a, b, c, d}, Reals, 
     Backsubstitution -> True];
  q = NMaximize[{Y/100*(0*a + d1*b + d2*c + 100*d), r}, {a, b, c, d}];
  l = q + s


  },
 {\[Placeholder]}
}

During evaluation of In[9]:= NMaximize::incst: NMaximize was unable to generate any initial points satisfying the inequality constraints {(101 a-Sqrt[101] Sqrt[Power[<<2>>] Plus[<<2>>]])/10100+b<=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.

Out[9]= {{{55845.7 Null^7, {Null^7 (20000 + (a -> 0.953808)), 
    Null^7 (20000 + (b -> 0.00446185)), 
    Null^7 (20000 + (c -> 0.0279999)), 
    Null^7 (20000 + (d -> 0.0137303))}}}, {\[Placeholder]}}
POSTED BY: Alex Graham
In[1]:= sln = 
 FindMaximum[{1000000/100*(25*b + 75*c + 100*(1 - a - b - c)), 
   a >= 0 && b >= 0 && c >= 0 && 
    a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 
      50*100*a^99*b^2 >= 1/20}, {{a, 0.95}, {b, 0.004}, {c, 0.03}}]

Out[1]= {35843.9, {a -> 0.953819, b -> 0.00444532, c -> 0.028011}}

In[2]:= {a, b, c} = {a, b, c} /. sln[[2]]

Out[2]= {0.953819, 0.00444532, 0.028011}
POSTED BY: Frank Kampas
Posted 6 years ago

Thanks guys for your replies! Could you give me some more advice on how to use the results for a,b,c... from Nmaximize or findmaximum in one input? I've tried to make like this

r = FindMaximum[{1000000/100*(25*b + 75*c + 100*(1 - a - b - c)), 
   a >= 0 && b >= 0 && c >= 0 && 
    a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 
      50*100*a^99*b^2 == 1/20}, {{a, 0.95}, {b, 0.004}, {c, 0.03}}]
s = 50 + a /. [[2]] + b /. [[3]]

But it seems to be working only with solve equation and I also tried to use the other way

{a, b, c} = 
  With[{r = 
     FindMaximum[{1000000/100*(25*b + 75*c + 100*(1 - a - b - c)), 
       a >= 0 && b >= 0 && c >= 0 && 
        a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 
          50*100*a^99*b^2 == 1/20}, {{a, 0.95}, {b, 0.004}, {c, 
        0.03}}]}, Extract[r, Position[r, _?NumericQ]]];
s = 50 + a + b

But both without success. Thanks in advance!

POSTED BY: Alex Graham

The problem can be solved exactly using the Kuhn-Tucker equations as a gigantic root object. I changed it to a minimization, as that is what my KTEqs function is set up for. I only show the numerical result here, which agrees with Jim's result.

In[8]:= 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, cons];
  allvars = Join[vars, lambdas, mus];
  {alleqns, allvars}
  ]

In[16]:= eqs = 
  KTEqs[-1000000/100*(25*b + 75*c + 100*(1 - a - b - c)), {a >= 0, 
    b >= 0, c >= 0, 
    a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 
      50*100*a^99*b^2 >= 1/20}, {a, b, c}];

In[17]:= r = Reduce[Sequence @@ eqs, Backsubstitution -> True];

In[18]:= N[r]

Out[18]= a == 0.953819 && b == 0.00444532 && 
 c == 0.028011 && \[Mu][1.] == 0. && \[Mu][2.] == 0. && \[Mu][3.] == 
  0. && \[Mu][4.] == 190922.
POSTED BY: Frank Kampas
Posted 6 years ago

If you substitute 1-a-b-c for d, remove the unnecessary 0*a, use better starting values, and use FindMaximum, you'll get the following very quickly:

FindMaximum[{1000000/100*(25*b + 75*c + 100*(1 - a - b - c)),  a >= 0 && b >= 0 && c >= 0 &&
   a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 50*100*a^99*b^2 >= 1/20}, 
   {{a, 0.95}, {b, 0.004}, {c, 0.03}}]

(* {35843.9, {a -> 0.953819, b -> 0.00444532, c -> 0.028011}} *)
POSTED BY: Jim Baldwin
NMaximize gives a better result but there are a lot of warning messages.
In[7]:= NMaximize[{1000000/100*(0*a + 25*b + 75*c + 100*d), 
  r && b >= 0 && c >= 00 && d >= 0}, {a, b, c, d}]

During evaluation of In[7]:= NMaximize::incst: NMaximize was unable to generate any initial points satisfying the inequality constraints {0.0244168 -1. c-1. d<=0,-0.0461073+1. c+1. d<=0,-c+(1-100000 (0.0461073 +Times[<<1>>]+Times[<<2>>])^2 Root[<<1>>&,2,0]^99-2020 (<<21>> -1. c-1. d) Root[<<1>>&,2,0]^100-20 Root[<<1>>&,2,0]^101)/(2020 Root[1+<<1>>+Times[<<2>>]&,2,0]^99 (100 (0.0461073 +Times[<<2>>]+Times[<<2>>])+Root[Plus[<<3>>]&,2,0]))<=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.

During evaluation of In[7]:= NMaximize::incst: NMaximize was unable to generate any initial points satisfying the inequality constraints {-b+(Sqrt[(1000+Times[<<2>>])/Plus[<<4>>]^99]-101 (1. -1. b-1. c-1. d))/10000<=0,0. -1. c-1. d<=0,0.0292252 -1. b-1. c-1. d<=0,-0.0461073+1. b+1. c+1. d<=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.

During evaluation of In[7]:= NMaximize::incst: NMaximize was unable to generate any initial points satisfying the inequality constraints {b-Sqrt[(1000-9799 Power[<<2>>])/(1. +Times[<<2>>]+<<1>>+Times[<<2>>])^99]/10000+(101 (1. -1. b-1. c-1. d))/10000<=0,-c+(1-100000 b^2 (<<1>>)^99-2020 b (<<1>>)^100-20 (<<1>>)^101)/(2020 (1. -1. b-1. c-1. d)^99 (1. +99. b-1. c-1. d))<=0,0.0292252 -1. b-1. c-1. d<=0,-0.0461073+1. b+1. c+1. d<=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.

During evaluation of In[7]:= General::stop: Further output of NMaximize::incst will be suppressed during this calculation.

During evaluation of In[7]:= Less::nord: Invalid comparison with 0.260607 -4.68869*10^17 I attempted.

During evaluation of In[7]:= NMaximize::bcons: The following constraints are not valid: {d==1-a-b-c,b>=0,c>=0,d>=0,a<Root[1-2020 #1^100+2000 #1^101&,2,0],(5050-5050 a-Sqrt[5] Sqrt[(-51+<<1>><<1>><<1>>+Times[<<2>>])/a^99])/5100<b,Root[-1+100000 #1^99-197980 #1^100+98000 #1^101&,1,0]<a,b<=(-101 a+Sqrt[(1000-9799 Power[<<2>>])/a^99])/10000,(1-20 a^101-2020 a^100 b-100000 a^99 b^2)/(2020 a^99 (a+100 b))<=c,c<=1-a-b}. Constraints should be equalities, inequalities, or domain specifications involving the variables.

During evaluation of In[7]:= Less::nord: Invalid comparison with 0.358884 -5.76557*10^46 I attempted.

During evaluation of In[7]:= Less::nord: Invalid comparison with 0.686662 +2814.48 I attempted.

During evaluation of In[7]:= General::stop: Further output of Less::nord will be suppressed during this calculation.

During evaluation of In[7]:= NMaximize::bcons: The following constraints are not valid: {d==1-a-b-c,b>=0,c>=0,d>=0,(5050-5050 a-Sqrt[5] Sqrt[(-51+<<2>>+Times[<<2>>])/a^99])/5100<b,b<(5050-5050 a+Sqrt[5] Sqrt[(-51+<<1>><<1>><<1>>+<<1>>)/a^99])/5100,Root[-51+5100500 #1^99-10097980 #1^100+4998500 #1^101&,1,0]<a,a<=Root[-1+100000 #1^99-197980 #1^100+98000 #1^101&,1,0],(1-20 a^101-2020 a^100 b-100000 a^99 b^2)/(2020 a^99 (a+100 b))<=c,c<=1-a-b}. Constraints should be equalities, inequalities, or domain specifications involving the variables.

Out[7]= {35083.5, {a -> 0.950184, b -> 0.00455738, c -> 0.0452588, 
  d -> 0}}
POSTED BY: Frank Kampas
In[8]:= r = 
  Reduce[a + b + c + d == 1 && a >= 0 && b >= 0 && c >= 0 && d >= 0 &&
     a^101 + 101*a^100*b + 101*a^100*c + 101*100*a^99*b*c + 
      50*100*a^99*b^2 >= 1/20, {a, b, c, d}, Reals, 
   Backsubstitution -> True];

In[9]:= FindMaximum[{1000000/100*(0*a + 25*b + 75*c + 100*d), r}, {a, 
  b, c, d}]

Out[9]= {34580.5, {a -> 0.953893, b -> 0., c -> 0.0461073, d -> 0.}}

Note that I changed 0.05 to 1/20 to give Reduce an infinite precision input

POSTED BY: Frank Kampas
Posted 6 years ago

Thanks for your help! That was the case of dollar unit sampling with two errors in cents. But I wonder how people found four probabilities a=0.9538 b=0.0045 c=0.0280 d=0.0137 on old computer ibm370 in the late seventies. I'm pretty sure that I specify all conditions correctly. What do you think on it? In theory b and c should't be zero

POSTED BY: Alex Graham
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract