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

Posted 5 months ago
691 Views
|
10 Replies
|
6 Total Likes
|
 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?
10 Replies
Sort By:
Posted 5 months ago
 This is missing some details. Please provide a full minimal example of what you are trying to do.
Posted 5 months 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 5 months ago
 Is this what you are looking for? And @@ Table[Subscript[P, i] \[Element] Integers, {i, 1, 10, 1}] 
Posted 5 months 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 5 months ago
 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 5 months 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 5 months ago
 I'd start by getting rid of Round[Q] and making a new symbol, say rQ, to denote it. Then specify that it must be integer valued. The reason is that Round[Q] will function in one of two (bad) ways. Either it will round something called Q, which is not in your variable list and hence will throw off NMinimize, or else it will make a "variable" out of Round[Q], ignoring that the expression involves an integer-valued function, and so solutions will not use integer values. (As best I can tell, it is this second behavior that is taking place.)After that it would be good to remove all the integer valued variables from inner denominators, to the extent possible. Even then, I'm not sure this will satisfy the constraint checking code. In the worst case: remove the specification that each p[j] s an integer, and replace all uses in objective and constraints with Round[p[j]] (but leave as p[j] in the variable list). Adjust range accordingly, adding 1/2 at each end, so that the endpoint integers are as likely to appear in rounding as the interior ones.Here is an example of what I have in mind, modifying the example I showed previously. m = 2; pvars = Array[p, m]; NMinimize[{12. + 3200/rQ + 0.1*(30 + 0.01*Max[Round[p[1]], Round[p[2]]]) + 0.004*(rQ/4 + 20*val1*Sqrt[5.1 + rQ/(2*Round[p[1]])]) + (2000*(Sqrt[ 5]*(1/(E^((1/2)*val1^2*(1.02 + rQ/(10*Round[p[1]])))* Sqrt[2*Pi]) - val1*(1 - (1/2)* Erfc[-((val1*Sqrt[1.02 + rQ/(10*Round[p[1]])])/ Sqrt[2])])* Sqrt[1.02 + rQ/(10*Round[p[1]])]) + (1/(E^(val1^2/2)* Sqrt[2*Pi]) - val1*(1 - (1/2)*Erfc[-(val1/Sqrt[2])]))* Sqrt[5.1 + rQ/(2*Round[p[1]])]))/rQ + 0.03125*rQ*(1/50 + 1/Round[p[1]] - 1/Round[p[2]]) + 50*(0.1*(0.01 + 30/Round[p[1]]) + 8/Round[p[1]] + 0.0001*Round[p[1]] + 0.1*(0.01 + 30/Round[p[2]]) + 8/Round[p[2]] + 0.0001*Round[p[2]]), rQ >= 1, {p[1] >= 100.5, p[2] >= 100.5}, {p[1] <= 300.5, p[2] <= 300.5}, Element[rQ, Integers], 2 - (1/2)*Erfc[-(val1/Sqrt[2])] - (1/2)* Erfc[-((val1*Sqrt[1.02 + rQ/(10*Round[p[1]])])/Sqrt[2])] == 0.00004*rQ}, {p[1], p[2], rQ, val1}] (* Out[13]= {44.0651204574, {p[1] -> 299.944762632, p[2] -> 294.935213071, rQ -> 151, val1 -> 2.70188836409}} *) It turns out this is a better result than what I could obtain by specifying integrality constraints for the p[j] variables. SO I decided to also treat rQ this way, and got further improvement. m = 2; pvars = Array[p, m]; NMinimize[{12. + 3200/Round[rQ] + 0.1*(30 + 0.01*Max[Round[p[1]], Round[p[2]]]) + 0.004*(Round[rQ]/4 + 20*val1*Sqrt[5.1 + Round[rQ]/(2*Round[p[1]])]) + (2000*(Sqrt[ 5]*(1/(E^((1/2)*val1^2*(1.02 + rQ/(10*Round[p[1]])))* Sqrt[2*Pi]) - val1*(1 - (1/2)* Erfc[-((val1*Sqrt[1.02 + Round[rQ]/(10*Round[p[1]])])/ Sqrt[2])])* Sqrt[1.02 + Round[rQ]/(10*Round[p[1]])]) + (1/(E^(val1^2/2)* Sqrt[2*Pi]) - val1*(1 - (1/2)*Erfc[-(val1/Sqrt[2])]))* Sqrt[5.1 + Round[rQ]/(2*Round[p[1]])]))/Round[rQ] + 0.03125*rQ*(1/50 + 1/Round[p[1]] - 1/Round[p[2]]) + 50*(0.1*(0.01 + 30/Round[p[1]]) + 8/Round[p[1]] + 0.0001*Round[p[1]] + 0.1*(0.01 + 30/Round[p[2]]) + 8/Round[p[2]] + 0.0001*Round[p[2]]), rQ >= .50001, {p[1] >= 100.5, p[2] >= 100.5}, {p[1] <= 300.5, p[2] <= 300.5}, 2 - (1/2)*Erfc[-(val1/Sqrt[2])] - (1/2)* Erfc[-((val1*Sqrt[1.02 + Round[rQ]/(10*Round[p[1]])])/Sqrt[2])] == 0.00004*Round[rQ]}, {p[1], p[2], rQ, val1}] (* Out[62]= {32.7615577316, {p[1] -> 114.053105392, p[2] -> 291.4999941, rQ -> 438.497460088, val1 -> 2.21897691355}} *) Don't forget to round the integer-valued variables to get the desired values.
Posted 5 months 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}]] 
 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.