Problem with a chemical optimization problem

Question

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}]

Answer 1

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[
    MapThread[
      Plot[funcion[#1, PP], {PP, 4.25, 4.75}, 
        ImageSize -> Medium, 
        PlotLabel -> #2] &,
      {yy1, Style[Row[{"yy1 = ", #}], "SR"] & /@ yy1}]]

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.