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]