# Problem with a chemical optimization problem [closed]

First of all, I want to say that I'm quite new to Mathematica. I searched a lot on this site to learn how to solve my problem, which this is the problem of my life! I want to minimize a function. It is a chemical problem regarding a famous equation of state called the Peng-Robinson equation of state.

Perhaps I'm did something wrong when defining my function or there is an error in the numerical minimization, but I can't understand where where the problem is. So I'll appreciate any suggestions,

(*Preparing variables*)
T = {100};
tc = {126.1, 190.6};
pc = {33.94, 46.04};(*MPa*)
ω = {0.040, 0.011};
R = 83.14472;
P = 4.119;
x1 = 0.5;
x2 = 1 - x1;
y1 = 0.958;
y2 = 1 - 0.958;
k12 = 0;
lunghezza = Length[T];
m12 = -0.26992*ω^2 + 1.54226*ω + 0.37464;
a12vet = Table[
0.45724*(R*tc[[j]])^2/
pc[[j]]*(1 + m12[[j]] (1 - (T[[i]]/tc[[j]])^(1/2)))^2, {i,
lunghezza}, {j, {1, 2}}];
b12 = 0.0778*R*tc/pc;
b1 = b12[[1]];
b2 = b12[[2]];
a1 = a12vet[[All, 1]];
a2 = a12vet[[All, 2]];
a12 = (a1*a2)^(1/2) (1 - k12);
a = (x1)^2*a1 + 2*x1*x2*a12 + (x2)^2*a2;
b = x1*b1 + x2*b2;

(*Function to minimize*)

funcion[yy1_, PP_] :=
Module[{yy2, A11, A22, B1, B2, A12, BVap, AVap, ALiq, BLiq},
yy2 = 1 - yy1;
A11 = (a1*PP)/(R^2*T^2);
A22 = (a2*PP)/(R^2*T^2);
B1 = (b1*PP)/(R*T);
B2 = (b2*PP)/(R*T);
A12 = Sqrt[A11*A22];
BVap = yy1*B1 + yy2*B2;
AVap = yy1^2 A11 + 2 yy1 yy2 A12 + yy2^2 A22;
ALiq = x1^2 A11 + 2*(x1*x2)* A12 + x2^2 A22;
BLiq = x1*B1 + x2*B2;
solzv =
NSolve[Z - 1/(1 - BVap/Z) +
AVap/BVap*BVap/Z/(1 + 2 *BVap/Z - (BVap/Z)^2) == 0, Z];
solsv = Z /. {solzv[[1]], solzv[[3]]};
solv = solsv[[1]];
solzl =
NSolve[Z - 1/(1 - BLiq/Z) +
ALiq/BLiq*BLiq/Z/(1 + 2 *BLiq/Z - (BLiq/Z)^2) == 0, Z];
solsl = Z /. {solzl[[1]], solzl[[3]]};
soll = solsl[[2]];
Subscript[lnϕ, 1 v] =
B1/BVap (solv - 1) - Log[solv - BVap] -
AVap/(BVap Sqrt[8])
Log[(solv + (1 + Sqrt[2]) BVap)/(
solv + (1 - Sqrt[2]) BVap)] ((2 (yy1 A11 + yy2 A12))/AVap - B1/
BVap);
Subscript[lnϕ, 2 v] =
B2/BVap (solv - 1) - Log[solv - BVap] -
AVap/(BVap Sqrt[8])
Log[(solv + (1 + Sqrt[2]) BVap)/(
solv + (1 - Sqrt[2]) BVap)] ((2 (yy1 A12 + yy2 A22))/AVap - B2/
BVap);
Subscript[ϕ, 1 v] = Exp[Subscript[lnϕ, 1 v]];
Subscript[ϕ, 2 v] = Exp[Subscript[lnϕ, 2 v]];
Subscript[lnϕ, 1 l] =
B1/BLiq (soll - 1) - Log[soll - BLiq] -
ALiq/(BLiq Sqrt[8])
Log[(soll + (1 + Sqrt[2]) BLiq)/(
soll + (1 - Sqrt[2]) BLiq)] ((2 (x1 A11 + x2 A12))/ALiq - B1/
BLiq);
Subscript[lnϕ, 2 l] =
B2/BLiq (soll - 1) - Log[soll - BLiq] -
ALiq/(BLiq Sqrt[8])
Log[(soll + (1 + Sqrt[2]) BLiq)/(
soll + (1 - Sqrt[2]) BLiq)] ((2 (x1 A12 + x2 A22))/ALiq - B2/
BLiq);
Subscript[ϕ, 1 l] = Exp[Subscript[lnϕ, 1 l]];
Subscript[ϕ, 2 l] = Exp[Subscript[lnϕ, 2 l]];
Abs[(Subscript[ϕ, 1 v] - Subscript[ϕ,
1 l]) + (Subscript[ϕ, 2 v] - Subscript[ϕ, 2 l])]
]

(*Minimization*)

(*Does my function pass all variables?*)

funcion[0.9186212148724066, 4]

(* {0.171251} *)

NMinimize[{funcion[yy1, PP], 0 < yy1 < 1 && 0 < PP < 10}, {yy1, PP}]


## closed as off-topic by Szabolcs, MarcoB, Yves Klett, m_goldberg, user9660 Jun 28 '16 at 4:25

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – MarcoB, Yves Klett
If this question can be reworded to fit the rules in the help center, please edit the question.

• Possible duplicate of What are the most common pitfalls awaiting new users? – Szabolcs Jun 27 '16 at 9:34
• Please see mathematica.stackexchange.com/questions/18393/… Next time please (1) Try to reduce your code to the smallest possible example that still shows the problem. The same problem can be demonstrated with a one-liner function that uses NSolve (2) state clearly what the problem is. – Szabolcs Jun 27 '16 at 9:36
• try this NMinimize[{funcion[yy1, PP][[1]], 0 < yy1 < 1 && 0 < PP < 10}, {yy1, PP}] – happy fish Jun 27 '16 at 9:40
• I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. – m_goldberg Jun 27 '16 at 23:49

Your problem is caused by your trying to use T both as vector and scalar at the same time. That is bound to cause problems. happy fish indicated a quick-and-dirty fix in a comment to your question. Here is an attempt to produce a solution the roots out the structural problem.

T = 100; (* modified *)
tc = {126.1, 190.6};
pc = {33.94, 46.04};
ω = {0.040, 0.011};
R = 83.14472;
P = 4.119;
x1 = 0.5;
x2 = 1 - x1;
y1 = 0.958;
y2 = 1 - 0.958;
k12 = 0;
(* lunghezza=Length[T]; *)
m12 = -0.26992*ω^2 + 1.54226*ω + 0.37464;
a12vet =
Table[
0.45724*(R*tc[[j]])^2/pc[[j]]*(1 + m12[[j]] (1 - (T/tc[[j]])^(1/2)))^2,
{j, {1, 2}}] (* modified *);
b12 = 0.0778*R*tc/pc;
b1 = b12[[1]];
b2 = b12[[2]];
a1 = a12vet[[1]] ;(* modified *)
a2 = a12vet[[2]]; (* modified *)
a12 = (a1*a2)^(1/2) (1 - k12);
a = (x1)^2*a1 + 2*x1*x2*a12 + (x2)^2*a2;
b = x1*b1 + x2*b2;

funcion[yy1_, PP_] :=
Module[{yy2, A11, A22, B1, B2, A12, BVap, AVap, ALiq, BLiq},
yy2 = 1 - yy1;
A11 = (a1 PP)/(R^2 T^2); A22 = (a2 PP)/(R^2 T^2);
B1 = (b1 PP)/(R T);
B2 = (b2 PP)/(R T);
A12 = Sqrt[A11 A22];
BVap = yy1 B1 + yy2 B2;
AVap = yy1^2 A11 + 2 yy1 yy2 A12 + yy2^2 A22;
ALiq = x1^2 A11 + 2 (x1 x2) A12 + x2^2 A22; BLiq = x1 B1 + x2 B2;
solzv =
NSolve[
Z - 1/(1 - BVap/Z) + (AVap BVap)/(BVap (Z (1 + (2 BVap)/Z -
(BVap/Z)^2))) == 0,
Z];
solsv = Z /. {solzv[[1]], solzv[[3]]}; solv = solsv[[1]];
solzl =
NSolve[
Z - 1/(1 - BLiq/Z) + (ALiq BLiq)/(BLiq (Z (1 + (2 BLiq)/Z -
(BLiq/Z)^2))) == 0,
Z];
solsl = Z /. {solzl[[1]], solzl[[3]]}; soll = solsl[[2]];
Subscript[lnϕ, 1 v] =
(B1 (solv - 1))/BVap - Log[solv - BVap] -
(AVap
Log[(solv + (1 + Sqrt[2]) BVap)/(solv + (1 - Sqrt[2]) BVap)]
((2 (yy1 A11 + yy2 A12))/AVap - B1/BVap)) /
(BVap Sqrt[8]);
Subscript[lnϕ, 2 v] =
(B2 (solv - 1))/BVap - Log[solv - BVap] -
(AVap
Log[(solv + (1 + Sqrt[2]) BVap)/(solv + (1 - Sqrt[2]) BVap)]
((2 (yy1 A12 + yy2 A22))/AVap - B2/BVap)) /
(BVap Sqrt[8]);
Subscript[ϕ, 1 v] = Exp[Subscript[lnϕ, 1 v]];
Subscript[ϕ, 2 v] = Exp[Subscript[lnϕ, 2 v]];
Subscript[lnϕ, 1 l] =
(B1 (soll - 1))/BLiq - Log[soll - BLiq] -
(ALiq
Log[(soll + (1 + Sqrt[2]) BLiq)/(soll + (1 - Sqrt[2]) BLiq)]
((2 (x1 A11 + x2 A12))/ALiq - B1/BLiq)) /
(BLiq Sqrt[8]);
Subscript[lnϕ, 2 l] =
(B2 (soll - 1))/BLiq - Log[soll - BLiq] -
(ALiq
Log[(soll + (1 + Sqrt[2]) BLiq)/(soll + (1 - Sqrt[2]) BLiq)]
((2 (x1 A12 + x2 A22))/ALiq - B2/BLiq)) /
(BLiq Sqrt[8]);
Subscript[ϕ, 1 l] = Exp[Subscript[lnϕ, 1 l]];
Subscript[ϕ, 2 l] = Exp[Subscript[lnϕ, 2 l]];
Abs[
(Subscript[ϕ, 1 v] - Subscript[ϕ, 1 l]) +
(Subscript[ϕ, 2 v] - Subscript[ϕ, 2 l])]]

funcion[0.9186212148724066, 4]
(* 0.171251 *)


Now NMinimize will find a minimum.

NMinimize[{funcion[yy1, PP], 0 < yy1 < 1 && 0 < PP < 10}, {yy1, PP}]


{3.25681*10^-6, {yy1 -> 0.833462, PP -> 9.44828}}

So there is no error in your definition of funcion, nor is there an error in your NMinimize expression. All your errors were committed in your preamble where you set up your constants and parameters.

I have taken the position here that it is better to treat T as a scalar. If you want to evaluate funcion for many values of T, I recommend introducing it as a argument to funcion.

funcion[yy1_, PP_, T_] := Module[...]


### Update

This result I show for NMinimize is simply wrong. NMinimize is not able handle funcion in the region defined by 0 < yy1 < 1 && 0 < PP < 10. Let's try to find out why. It is always helpful to visualize a computationally complex like funcion to get an understanding how it behaves over rhe domain of interest.

Plot3D[funcion[yy1, PP], {yy1, 0, 1}, {PP, 2.5, 10},
PlotRange -> {0, 1},
ImageSize -> Large]


Looking at the 3D plot from various angles, allows certain inferences.

• funcion depends rather weakly on yy1.

• funcion has infinite number of minima in the specified domain. The minima are all have value zero. They are all located along a curve where the derivative of fucion w.r.t. PP has a discontinuity caused by appearance of Abs in the last line of your definition of funcion. NMinimize was not designed to handle this kind of minimum. The result you are getting in whatever version of Mathematica is actully more informative than the wrong answer I got from V10.4.1.

• In addition to the discontinuity in the derivitive caused by

Abs[
(Subscript[ϕ, 1 v] - Subscript[ϕ, 1 l]) +
(Subscript[ϕ, 2 v] - Subscript[ϕ, 2 l])]


the subtractions cause loss of precision near the zeros of funcion, which results in numerical problems. It is unlikely, because funcion is defined in terms of machine arithmetic, that the zeros can be located very precisely.

Since funcion depends weakly on yy1, we can get some idea of its behavior near the zeros by looking at slices taken for some fixed values of yy1. These computations have the advantage of being quicker then computations based on the whole 3D surface.

Off[Power::infy]; Off[Infinity::indet];
With[{yy1 = {0., .5, 1.}},
GraphicsColumn[

I do not know if this behavior is what you expected, but I guess not. I guess you thought there would a nice local minimum of your function it the domain you defined. Well, there isn't. The problem may be that funcion doesn't model the chemistry correctly, but I don't know the chemistry involved and have nothing further to say on this issue.
• @Cecil. I am running Mathematica V10.4.1 on OS X. I have determined that the result NMinimize gave is wrong. I am looking further into the problem. Will post my findings when I feel they are complete.. – m_goldberg Jun 28 '16 at 22:29