Group Abstract Group Abstract

Message Boards Message Boards

0
|
5.8K Views
|
10 Replies
|
6 Total Likes
View groups...
Share
Share this post:

Create a list of integer variables to Minimize a particular quantity?

Elisa Boffelli
Posted 7 years ago

Hi everyone and thanks in advance for your help. I used a list of variables to Minimize a particular quantity

Flatten[{ Table[Subscript[P, i], {i, 1, m, 1}], Round[Q], val1}]

How can I define that 

 Table[Subscript[P, i], {i, 1, m, 1}]

should be integers?? Any quick idea?
POSTED BY: Elisa Boffelli
10 Replies
Sort By:

I believe there is a bug in the constraint checking code and I will file a bug report on that. A workaround involves pulling one of the variables outside a particular square root. I show the recoded variant below. I use p[1] and the like instead of subscripts, and made a few other alterations that I believe do not affect the outcome. Below is the altered loop. It terminates for me after two steps.

While[ctotparz <= ctotfin && m < 4, m++;
  pvars = Array[p, m];
  min = NMinimize[{(Av*Edom)/rQ + (m*Ar*Edom)/
       rQ + ((m*(m - 1))/Edom + 1/p[1] - 
         Sum[(Sum[2/p[j], {j, 2, i - 1}] + 1/p[i]), {i, 2, m}])*(Edom*
          hv*rQ)/(2*m^2) + 
      Edom/m*(Sum[(g/p[i] + b*(p[i])^beta + (w0/p[i] + k)*e), {i, 1, 
          m}]) + (w0 + Max[pvars]*k)*o + alfa*Edom*Cv + 
      Edom/rQ*(m*Ab + m*Z) + 
      Cb*Edom + (rQ/(2*m) + val1*dev*Sqrt[rQ/(m*p[1]) + delta + ton])*
       hb + (gamma*Edom*dev)/
        rQ*(Sqrt[
           rQ/(m*p[1]) + delta + 
            ton]*(PDF[NormalDistribution[0, 1], val1] - 
            val1*(1 - CDF[NormalDistribution[0, 1], val1])) + (m - 1)*
          Sqrt[delta]*(PDF[
             NormalDistribution[0, 
              1], (val1*
               Sqrt[rQ/(m*p[1]*delta) + ton/delta + 1])] - (val1*
               Sqrt[rQ/(m*p[1]*delta) + ton/delta + 1])*(1 - 
               CDF[NormalDistribution[0, 
                 1], (val1*
                  Sqrt[rQ/(m*p[1]*delta) + ton/delta + 1])]))), 
     Flatten[{rQ >= 1, Thread[pvars >= 101], Thread[pvars <= 300], 
       Element[pvars, Integers], 
       Element[rQ, 
        Integers], (1 - 

          CDF[NormalDistribution[0, 1], 
           val1] + (m - 1)*(1 - 
             CDF[NormalDistribution[0, 
               1], (val1*
                Sqrt[rQ/(m*delta) + p[1]*(ton/delta + 1)]/
                 Sqrt[p[1]])])) == (hb*rQ)/(gamma*Edom)}]}, 
    Flatten[{pvars, rQ, val1}]];
  ctotparz = min[[1]];
  y = min[[2]];
  Pparz = pvars /. y;
  Qparz = rQ /. y;
  val1parz = val1 /. y;
  If[ctotparz <= ctotfin, (ctotfin = ctotparz; Qfin = Qparz;
    Pfin = Pparz)]];

Suffice it to say that this sort of thing is impossible to diagnose in the absence of the unabridged code being used.
POSTED BY: Daniel Lichtblau
Elisa Boffelli
Elisa Boffelli
Posted 7 years ago

Thank you for your help, in this way that model works.

I tried to apply the same mechanism to an other model a little more complex but Mathematica reports the same error. Is it possible that for this model pulling one of the variables outside a particular square root doesn't work?

Are days that I'm trying to solve this problem, I'm going crazy so really thank you.

The model is:

Clear["Global`*"];
Av = 10;
Ab = 4;
b = 0.0001;
beta = 1;
Cv = 0.01;
Cb = 0.1;
Edom = 100;
g = 8;
hv = 0.0025;
hvs = 0.000625;
hvf = 0.001875;
hb = 0.004;
hbs = 0.0015;
hbf = 0.0025;
Z = 5;
alfa = 2;
delta = 5;
gamma = 1;
dev = 20;
w0 = 30;
k = 0.01;
e = 0.1;
ton = 0.1;
o = 0.1;
Ar = 2;
imax = 2000;
ctotparz = 100000000000000;
ctotfin = 100000000000000;


n = 0;
While[ctotparz <= ctotfin,
  n++;
  pvar = Array[P, n];
  ctotparz = 100000000000000; 
  For[rit = 0, rit < n, rit++,
   x = NMinimize[{(Av*Edom)/Round[Q] + (n*Ar*Edom)/Round[Q] + 
       Edom/Round[
         Q]*((rit + 1)*Z) + ((rit*(rit + 1))/Edom + 
          Sum[1/P[i], {i, 1, (n - rit)}] - 
          Sum[(Sum[1/P[j], {j, (n - rit + 1), (i - 1)}] + 
             Sum[1/P[j], {j, (n - rit + 1), (i)}]), {i, (n - rit + 
              1), (n)}] )*(Edom*hv*Round[Q])/(2*n^2) + 
       Edom/n* Sum[(g/P[i] + b*(P[i])^beta + (w0/P[i] + k)*e), {i, 1, 
          n}] + (w0 + Max[pvar]*k)*o + alfa*Edom*Cv + 
       Edom/Round[Q]*(n*Ab) + 
       Cb*Edom + (Round[Q]*Edom)/(
        2*n^2)*(hbs + 
          hvf)*(Edom* Sum[1/(P[i])^2, {i, 2, (n - rit)}] + 
          2* Sum[(1/P[i]*((i - 1) - Edom*Sum[1/P[j], {j, 2, i}])), {i,
              2, (n - rit)}] + 
          rit/Edom*(2*n - 2*rit - 1 - 
             2*Edom* Sum[1/P[i], {i, 2, (n - rit)}]) + 
          1/Edom*(n - rit - Edom* Sum[1/P[i], {i, 2, (n - rit)}])^2) +
        val1*dev Sqrt[Round[Q]/(n*P[1]) + delta + ton]*hb + (
        gamma*Edom*dev)/
        Round[Q]*(Sqrt[
           Round[Q]/(n*P[1]) + delta + 
            ton]*(PDF[NormalDistribution[0, 1], val1] - 
             val1*(1 - CDF[NormalDistribution[0, 1], val1])) + 
          Sum[(Sqrt[Round[Q]/(
             n*P[i])]*(PDF[
                NormalDistribution[0, 1], ((
                 Round[Q]/n*((i - 1) - Edom* Sum[1/P[j], {j, 2, i}]) +
                   val1*dev*Sqrt[Round[Q]/(n*P[1]) + delta + ton])/(
                 dev*Sqrt[Round[Q]/(n*P[i])]))] - ((
                 Round[Q]/n*((i - 1) - Edom*Sum[1/P[j], {j, 2, i}]) + 
                  val1*dev*Sqrt[Round[Q]/(n*P[1]) + delta + ton])/(
                 dev*Sqrt[Round[Q]/(n*P[i])]))*(1 - 
                  CDF[NormalDistribution[0, 1], (1/(
                    dev*Sqrt[Round[Q]/(

                    n*P[i])])(Round[Q]/
                    n*((i - 1) - Edom*Sum[1/P[j], {j, 2, i}]) + 
                    val1*dev*Sqrt[
                    Round[Q]/(n*P[1]) + delta + ton]))]))), {i, 
            2, (n - rit)}] + 
          rit*Sqrt[
           delta]*(PDF[
              NormalDistribution[0, 1], (1/(
               dev*Sqrt[
                delta])(Round[Q]/
                  n*(n - rit - 1 - 
                    Edom*Sum[1/P[i], {i, 2, (n - rit)}]) + 
                 val1*dev*Sqrt[
                  Round[Q]/(n*P[1]) + delta + ton]))] - (1/(
               dev*Sqrt[
                delta])(Round[Q]/
                  n*(n - rit - 1 - 
                    Edom*Sum[1/P[i], {i, 2, (n - rit)}]) + 
                 val1*dev*Sqrt[Round[Q]/(n*P[1]) + delta + ton]))*(1 -
                 CDF[NormalDistribution[0, 1], (1/(
                  dev*Sqrt[
                   delta])(Round[Q]/
                    n*(n - rit - 1 - 
                    Edom*Sum[1/P[i], {i, 2, (n - rit)}]) + 
                    val1*dev*Sqrt[
                    Round[Q]/(n*P[1]) + delta + ton]))]))),  
      Flatten[{Round[Q] >= 1, Thread[pvar >= 101], 
        Thread[pvar <= 300], 
        Element[pvar, 
         Integers],  (val1*dev*Sqrt[Round[Q]/n + P[1]*(delta + ton)]/
            Sqrt[P[1]] + 
           Edom*Sqrt[Round[Q]/n + P[1]*(delta + ton)]/Sqrt[P[1]] + 
           Round[Q]/n*(n - rit) ) <= 
         imax, (1 - CDF[NormalDistribution[0, 1], val1] + 

           Sum[(1 - 
              CDF[NormalDistribution[0, 1], ((
                Round[Q]/n*((i - 1) - Edom*Sum[1/P[j], {j, 2, i}]) + 
                 val1*dev*Sqrt[Round[Q]/n + P[1]*(delta + ton)]/Sqrt[
                  P[1]])/(dev*Sqrt[Round[Q]/n]/Sqrt[P[i]]))]), {i, 
             2, (n - rit)}] + 
           rit*(1 - 
              CDF[NormalDistribution[0, 1], (1/(
                dev*Sqrt[
                 delta])(Round[Q]/
                   n*(n - rit - 1 - 
                    Edom*Sum[1/P[i], {i, 2, (n - rit)}]) + 
                  val1*dev*Sqrt[Round[Q]/n + P[1]*(delta + ton)]/Sqrt[
                   P[1]]))])) == (hb*Round[Q])/(gamma*Edom)}]}, 
     Flatten[{pvar, Round[Q], val1}]];
   ctotnew  = x[[1]];
   y = x[[2]];
   Pnew = pvar /. y;
   Qnew = Round[Round[Q] /. y];
   val1new = val1 /. y;
   If[ctotnew <= ctotparz, ( ctotparz = ctotnew;  Qparz = Qnew;  
     Pparz = Pnew; ritparz = rit), ctotnew = ctotnew]];
  If[ctotparz <= ctotfin, ( ctotfin = ctotparz;  Qfin = Qparz;  
    Pfin = Pparz; ritfin = ritparz), ctotparz = ctotparz]];
nfin = n - 1; 
Print[ctotfin];
Print[Qfin] ;
Print[Pfin];
Print[nfin];
Print[ritfin]
POSTED BY: Elisa Boffelli
POSTED BY: Daniel Lichtblau
Elisa Boffelli
Elisa Boffelli
Posted 7 years ago

Thank you, as you suggested I used Round on the variable P[i] and Q, I removed also a constraint and I replaced it with an other one simpler, in this way I obtained a better result.

I have still problems with the other model, the one more complex, I adopted the same procedure and it works for n=1, but when I try to implement n, the program is not able to evaluate a result.

Thank ou very much for your help.

This is the modified code:

Clear["Global`*"];
Av = 10;
Ab = 4;
b = 0.0001;
beta = 1;
Cv = 0.01;
Cb = 0.1;
Edom = 100;
g = 8;
hv = 0.0025;
hvs = 0.000625;
hvf = 0.001875;
hb = 0.004;
hbs = 0.0015;
hbf = 0.0025;
Z = 5;
alfa = 2;
delta = 5;
gamma = 1;
dev = 20;
w0 = 30;
k = 0.01;
e = 0.1;
ton = 0.1;
o = 0.1;
Ar = 2;
imax = 3000;
ctotparz = 100000000000000;
ctotfin = 100000000000000;


n = 2;
rit = 0;
pvar = Array[P, n];

NMinimize[{(Av*Edom)/Round[Q] + (n*Ar*Edom)/Round[Q] + 
   Edom/Round[
     Q]*((rit + 1)*Z) + ((rit*(rit + 1))/Edom + 
      Sum[1/Round[P[i]], {i, 1, (n - rit)}] - 
      Sum[(Sum[1/Round[P[j]], {j, (n - rit + 1), (i - 1)}] + 
         Sum[1/Round[P[j]], {j, (n - rit + 1), (i)}]), {i, (n - rit + 
          1), (n)}] )*(Edom*hv*Round[Q])/(2*n^2) + 
   Edom/n* Sum[(g/Round[P[i]] + 
       b*(Round[P[i]])^beta + (w0/Round[P[i]] + k)*e), {i, 1, 
      n}] + (w0 + Max[Thread[Round[pvar]]]*k)*o + alfa*Edom*Cv + 
   Edom/Round[Q]*(n*Ab) + 
   Cb*Edom + (Round[Q]*Edom)/(
    2*n^2)*(hbs + 
      hvf)*(Edom* Sum[1/(Round[P[i]])^2, {i, 2, (n - rit)}] + 
      2* Sum[(1/
          Round[P[i]]*((i - 1) - 
            Edom*Sum[1/Round[P[j]], {j, 2, i}])), {i, 2, (n - rit)}] +
       rit/Edom*(2*n - 2*rit - 1 - 
         2*Edom* Sum[1/Round[P[i]], {i, 2, (n - rit)}]) + 
      1/Edom*(n - rit - 
         Edom* Sum[1/Round[P[i]], {i, 2, (n - rit)}])^2) + 
   val1*dev Sqrt[Round[Q]/(n*Round[P[1]]) + delta + ton]*hb + (
    gamma*Edom*dev)/
    Round[Q]*(Sqrt[
       Round[Q]/(n*Round[P[1]]) + delta + 
        ton]*(PDF[NormalDistribution[0, 1], val1] - 
         val1*(1 - CDF[NormalDistribution[0, 1], val1])) + 
      Sum[(Sqrt[Round[Q]/(
         n*Round[P[i]])]*(PDF[
            NormalDistribution[0, 1], ((
             Round[Q]/
               n*((i - 1) - Edom* Sum[1/Round[P[j]], {j, 2, i}]) + 
              val1*dev*Sqrt[Round[Q]/(n*Round[P[1]]) + delta + ton])/(
             dev*Sqrt[Round[Q]/(n*Round[P[i]])]))] - ((
             Round[Q]/
               n*((i - 1) - Edom*Sum[1/Round[P[j]], {j, 2, i}]) + 
              val1*dev*Sqrt[Round[Q]/(n*Round[P[1]]) + delta + ton])/(
             dev*Sqrt[Round[Q]/(n*Round[P[i]])]))*(1 - 

              CDF[NormalDistribution[0, 
                1], (1/(
                 dev*Sqrt[Round[Q]/(
                  n*Round[P[i]])]) (Round[Q]/
                    n*((i - 1) - Edom*Sum[1/Round[P[j]], {j, 2, i}]) +
                    val1*dev*Sqrt[
                    Round[Q]/(n*Round[P[1]]) + delta + ton]))]))), {i,
         2, (n - rit)}] + 
      rit*Sqrt[
       delta]*(PDF[
          NormalDistribution[0, 
           1], (1/(
            dev*Sqrt[
             delta]) (Round[Q]/
               n*(n - rit - 1 - 
                 Edom*Sum[1/Round[P[i]], {i, 2, (n - rit)}]) + 
              val1*dev*Sqrt[
               Round[Q]/(n*Round[P[1]]) + delta + ton]))] - (1/(
            dev*Sqrt[
             delta]) (Round[Q]/
               n*(n - rit - 1 - 
                 Edom*Sum[1/Round[P[i]], {i, 2, (n - rit)}]) + 
              val1*dev*Sqrt[
               Round[Q]/(n*Round[P[1]]) + delta + ton]))*(1 - 
            CDF[NormalDistribution[0, 
              1], (1/(
               dev*Sqrt[
                delta]) (Round[Q]/
                  n*(n - rit - 1 - 
                    Edom*Sum[1/Round[P[i]], {i, 2, (n - rit)}]) + 
                 val1*dev*Sqrt[
                  Round[Q]/(n*Round[P[1]]) + delta + ton]))]))),  
  Flatten[{Q >= 0.5000001, Thread[pvar >= 100.5000001], 
    Thread[pvar <= 
      300.5], (val1*dev*Sqrt[Round[Q]/n + Round[P[1]]*(delta + ton)]/
        Sqrt[Round[P[1]]] + 
       Edom*Sqrt[Round[Q]/n + Round[P[1]]*(delta + ton)]/Sqrt[
        Round[P[1]]] + Round[Q]/n*(n - rit) ) <= imax, val1 >= -3, 
    val1 <= 3}]}, Flatten[{pvar, Q, val1}]]
POSTED BY: Elisa Boffelli

Since P[1] is guaranteed to be positive, just clear the denominator in the inequality constraint:

(val1*dev*Sqrt[Round[Q]/n + Round[P[1]]*(delta + ton)] + 
Edom*Sqrt[Round[Q]/n + Round[P[1]]*(delta + ton)] + 
Sqrt[Round[P[1]]]*Round[Q]/n*(n - rit)) <= imax*Sqrt[Round[P[1]]]

With that change it seems to run fine.
POSTED BY: Daniel Lichtblau
Elisa Boffelli
Elisa Boffelli
Posted 7 years ago

Perfect, so when I have some complex costraints I have to try to semplify them as much as possibile and avoid that the variables are in the denominator position.

Thank you for your help, now I understand better how Minimize works!
POSTED BY: Elisa Boffelli
Rohit Namjoshi
Rohit Namjoshi
Posted 7 years ago

Is this what you are looking for?

And @@ Table[Subscript[P, i] \[Element] Integers, {i, 1, 10, 1}]

enter image description here
POSTED BY: Rohit Namjoshi
Elisa Boffelli
Elisa Boffelli
Posted 7 years ago

yes, thankyou for the suggestion, but if I try to insert it in my full code the programm reports an error: 

Remove["Global`*"];
Av = 10;
Ab = 4;
b = 0.0001;
beta = 1;
Cv = 0.01;
Cb = 0.1;
Edom = 100;
g = 8;
hv = 0.0025;
hb = 0.004;
Z = 5;
alfa = 2;
delta = 5;
gamma = 1;
dev = 20;
w0 = 30;
k = 0.01;
e = 0.1;
ton = 0.1;
o = 0.1;
Ar = 2;
ctotfin = 100000000000000;
ctotparz = 100000000000000;

m = 0;
While[ctotparz <= ctotfin,
  m++;
  x = NMinimize[{(Av*Edom)/Round[Q] + (m*Ar*Edom)/
      Round[Q] + ((m*(m - 1))/Edom + 1/Subscript[P, 1] - 
         Sum[(Sum[1/Subscript[P, j], {j, 2, i - 1}] + 
            Sum[1/Subscript[P, j], {j, 2, i}]), {i, 2, m}])*(
       Edom*hv*Round[Q])/(2*m^2) + 
      Edom/m*(Sum[(g/Subscript[P, i] + 
           b*(Subscript[P, i])^beta + (w0/Subscript[P, i] + k)*
            e), {i, 1, m}]) + (w0 + 
         Max[Table[Subscript[P, i], {i, 1, m, 1}]]*k)*o + 
      alfa*Edom*Cv + Edom/Round[Q]*(m*Ab + m*Z) + 
      Cb*Edom + (Round[Q]/(2*m) + 
         val1*dev*Sqrt[Round[Q]/(m*Subscript[P, 1]) + delta + ton])*
       hb + (gamma*Edom*dev)/
       Round[Q]*(Sqrt[

          Round[Q]/(m*Subscript[P, 1]) + delta + 
           ton]*(PDF[NormalDistribution[0, 1], val1] - 
            val1*(1 - CDF[NormalDistribution[0, 1], val1])) + (m - 1)*
          Sqrt[delta]*(PDF[
             NormalDistribution[0, 
              1], (val1*Sqrt[
               Round[Q]/(m*Subscript[P, 1]*delta) + ton/delta + 
                1])] - (val1*Sqrt[
               Round[Q]/(m*Subscript[P, 1]*delta) + ton/delta + 
                1])*(1 - 
               CDF[NormalDistribution[0, 
                 1], (val1*Sqrt[
                  Round[Q]/(m*Subscript[P, 1]*delta) + ton/delta + 
                   1])]))),  
     Round[Q] > 0   && Table[Subscript[P, i], {i, 1, m, 1}] >= 101 && 
      Table[Subscript[P, i], {i, 1, m, 1}] <= 300 && 
      And @@ Table[
        Subscript[P, i] \[Element] Integers, {i, 1, m, 1}] && (1 - 
         CDF[NormalDistribution[0, 1], 
          val1] + (m - 1)*(1 - 
            CDF[NormalDistribution[0, 
              1], (val1*Sqrt[
               Round[Q]/(m*Subscript[P, 1]*delta) + ton/delta + 
                1])])) == (hb*Round[Q])/(gamma*Edom) }, 
    Flatten[{ Table[Subscript[P, i], {i, 1, m, 1}], Round[Q], val1}]];
  ctotparz = x[[1]];
  y = x[[2]];
  Pparz = Table[Subscript[P, i], {i, 1, m, 1}] /. y[[1 ;; m]];
  Qparz = Round[Round[Q] /. y[[m + 1]]];
  val1parz = val1 /. y[[m + 2]];
  If[ctotparz <= ctotfin, ( ctotfin = ctotparz;  Qfin = Qparz;  
    Pfin = Pparz), ctotparz = ctotparz]];
mfin = m - 1; 
Print[ctotfin];
Print[Qfin] ;
Print[Pfin];
Print[mfin]

Power::infy: Infinite expression 1/0 encountered.

NMinimize:: The constraints are not valid. Constraints should be equalities, inequalities, or domain specifications involving the variables.

Any idea of the possible reason?
POSTED BY: Elisa Boffelli
Elisa Boffelli
Elisa Boffelli
Posted 7 years ago

Ok, so I want to add the constraint that the elements of the vector

Table[Subscript[P, i], {i, 1, m, 1}]

are Integers. This is a part of my model:

Clear["Global`*"];
Edom = 100;
hv = 0.0025;
Av = 10;
Ar = 2;
ctotfin = 100000000000000;
ctotparz = 100000000000000;

m = 0;
While[ctotparz <= ctotfin,
  m++;
  x = NMinimize[{ ((Av*Edom)/Round[Q] + (m*Ar*Edom)/Round[Q] + (
        m*(m - 1))/Edom + 1/Subscript[P, 1] - 
        Sum[(Sum[1/Subscript[P, j], {j, 2, i - 1}] + 
           Sum[1/Subscript[P, j], {j, 2, i}]), {i, 2, m}])*(
      Edom*hv*Round[Q])/(2*m^2),  
     Round[Q] > 0   && Table[Subscript[P, i], {i, 1, m, 1}] >= 101 && 
      Table[Subscript[P, i], {i, 1, m, 1}] <= 300 }, 
    Flatten[{ Table[Subscript[P, i], {i, 1, m, 1}], Round[Q], val1}]];
  ctotparz = x[[1]];
  y = x[[2]];
  Pparz = Table[Subscript[P, i], {i, 1, m, 1}] /. y[[1 ;; m]];
  Qparz = Round[Round[Q] /. y[[m + 1]]];
  val1parz = val1 /. y[[m + 2]];
  If[ctotparz <= ctotfin, ( ctotfin = ctotparz;  Qfin = Qparz;  
    Pfin = Pparz), ctotparz = ctotparz]];
POSTED BY: Elisa Boffelli

This is missing some details. Please provide a full minimal example of what you are trying to do.
POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard