Error using fsolve (line 269) FSOLVE requires all values returned by functions to be of data type double.

clearvars
close all
clc

global nu F_in ai bi T_op P_op R Keq Ai Bi

syms T 

% Reference temperature [K]

Trif = 298.15; 

% Operating temperature [K] 

T_op = 250 + 273.15;

% Operating Pressure    [Pa]

P_op = 50*1e5;

% Gas constant [J/(mol*K)]

R = 8.31445;

% Total inlet molar flow rate [mol/s]

F_intot = 100;

% Inlet composition []

x_in = [0     0.5      0      0    0.5     0];
        %CO2  %H2   %CH3OH  %H2O   %CO %CH3OCH3 

% Inlet flow rate [mol/s]

F_in = F_intot*x_in;

% Pitzer acentric factor [] 2-138

w = [0.2236  -0.2160  0.5658  0.3449  0.0482  0.2002];

% Critical Temperatures [K] 2-138

Tc = [304.21   33.19   512.5   647.096   132.92   400.1];

% Reduced Temperatures

Tr = T_op./Tc;

% Critical Pressures [Pa] 2-138

Pc = 1e6*[7.383   1.313   8.084   22.064   3.499   5.37]; 

% Reduced Pressures

Pr = P_op./Tc;

% S, ai, bi, Ai, Bi factors

S = 0.48 + 1.54.*w - 0.176.*(w.^2);
bi = 0.08664*R.*Tc./Pc;
k = 1 + S.*(1 - sqrt(Tr)).^2;
ai = 0.42748*R.*(Tc.^2)./Pc;
Ai = (ai.*P_op)/(R*T_op)^2;
Bi = (bi.*P_op)/(R*T_op);

% Stoichiometric coefficients

nu = [-1   -3     1     1    0     0        %R1  CO2 + 3H2 -> CH3OH + H2O
       1    1     0    -1   -1     0        %R3  CO + H2O -> CO2 + H2 
       0    0    -2     1    0     1];      %R4  2CH3OH -> CH3OCH3 + H2O
      %CO2  %H2 %CH3OH %H2O  %CO %CH3OCH3

% Enthalpies of formation (Trif = 298.15 K) [J/mol]

DH_0f = 1e4*[-39.351  0  -20.094  -24.1814  -11.053  -18.41];
               %CO2  %H2  %CH3OH     %H2O     %CO   %CH3OCH3

DH_0f = nu*DH_0f';

% Gibbs free energies (Trif = 298.15 K) [J/mol]

DeltaG_0 = 1e4*[-39.437  0  -16.232  -22.859  -13.715  -11.28];
                  %CO2  %H2  %CH3OH     %H2O     %CO   %CH3OCH3

DeltaG_0 = nu*DeltaG_0';

% Equilibrium constant at reference conditions

Keq_0 = exp(DeltaG_0/(R*Trif));

% Specific Heats [J/(mol*K)]

cpi = [2.9370e+04   2.7617e+04   3.9252e+04   3.3363e+04   2.9108e+04   5.1480e+04     % C1
       3.4540e+04   9.5600e+03   8.7900e+04   2.6790e+04   8.7730e+03   1.4420e+05     % C2 
       1.4280e+03   2.4660e+03   1.9165e+03   2.6105e+03   3.0851e+03   1.6034e+03     % C3
       2.6400e+04   3.7600e+03   5.3654e+04   8.8960e+03   8.4553e+03   7.7470e+04     % C4
       5.8800e+02   5.6760e+02   8.9670e+02   1.1690e+03   1.5382e+03   7.2540e+02];   % C5
           %CO2         %H2         %CH3OH         %H2O         %CO       %CH3OCH3

for j = 1:6

cp(j) = ((cpi(1,j) + cpi(2,j)*((cpi(3,j)/T)/(sinh(cpi(3,j)/T)))^2 + cpi(4,j)*...
    ((cpi(5,j)/T)/(cosh(cpi(5,j)/T)))^2))*1e-3;

end

Delta_cp = nu*cp';

% Vant't Hoff equation (Temperature dependece)

DH_0 = DH_0f + int(Delta_cp,Trif,T);
f = (DH_0/(R*T^2));
Keq = Keq_0 + exp(int(f,Trif,T));

% Graphical Representation

figure(1)
ezplot(Keq(1), [250, 450, 0, 5])
hold on
ezplot(Keq(2), [250, 450, 0, 5])
hold on
ezplot(Keq(3), [250, 450, 0, 5])
xlabel('T [K]');
ylabel('Keq [-]');
legend('Keq1','Keq2','Keq3')
title('Equilibrium Reactions Constants');

% Vant't Hoff equation 

Keq = Keq_0 + exp(int(f,Trif,T_op));

% Equation solution

inval = [1 2 3];
lambda = fsolve(@TRONC_2,inval)

function F = TRONC_2(x)

global nu F_in ai bi T_op P_op R Keq Ai Bi

% Unknowns (progress grades of reactions)

l1 = x(1);
l2 = x(2);
l3 = x(3);

% Outlet streams [mol/s]

l = [l1;l2;l3];
l_out = nu'*l;
F_out = F_in' + l_out;
F_outtot = sum(F_out);

% Outlet compositions

y_out = F_out/F_outtot;

% Mixing Rule ("a")

for i = 1:6
    for j = 1:6
aa(i,j) = y_out(i)*y_out(j)*sqrt(ai(i)*ai(j));
    end
end
a = sum(sum(aa));

% Mixing Rule ("b")

for i = 1:6
bb(i) = y_out(i)*bi(i); 
end

b = sum(bb);

% Parameters A and B

A = (a.*P_op)/(R*T_op)^2;
B = (b.*P_op)/(R*T_op);

% Compressibility factor

alpha = -1;
beta = A - B - B^2;
gamma = -A*B;
p = [1 alpha beta gamma];
Z = roots(p);  
Z = Z(3); % A cazzo (deve essere il maggiore dei tre, come fare?)

% Fugacity coefficient

phi = exp((Z-1)*(Bi./B) + (A/B)*(Bi./B - 2*sqrt(Ai/A))*log((Z + B)/Z) - log(Z - B));

% Equations

F(1) = Keq(1) - 1/P_op*(y_out(3)*y_out(4))/(y_out(1)*(y_out(2))^3);
F(2) = Keq(2) - (y_out(1)*y_out(2))/(y_out(3)*y_out(4));
F(3) = Keq(3) - (y_out(6)*y_out(4))/(y_out(3))^2;

F = F';

end

Answer 1

[Walter Roberson]

In the line

Keq = Keq_0 + exp(int(f,Trif,T_op));

you are doing a symbolic integration (f involves the symbol T) and getting a closed-form symbolic result . This hints that you are looking for closed form solutions for the equations.

Keq is a global variable.

In the line

lambda = fsolve(@TRONC_2,inval)

TRONC_2 is a function that involves the global variable Keq. Calculations that have symbolic values in the mix give symbolic outputs. Therefore TRONC_2 returns symbolic outputs. fsolve() refuses to work with symbolic values.

If you are looking for a numeric output, then you need to consider whether you need the extended precision of symbolics through the body of TRONC_2 but want to double() as the last step of that function. If not, then you should consider double() of what you assign to Keq. Or you should consider doing a numeric integral of f. As you are constructing f symbolically, if you do want to switch to numeric integral for it, consider using matlabFunction() to generate the function to be numerically integrated.

You should also consider not using fsolve. You have a function that returns symbolic values at present, so you can use

x = sym('x',[1 3]);
temp = TRONC_2(x);
lamba = vpasolve(temp);

One thing I do not understand is the Z = Z(3) "deve essere il maggiore dei tre, come fare" line. I see you take roots, but then you count on the result you want being the third root. Why that one? Are there particular properties of the root that are needed? This makes a difference especially if you do try the symbolic route, since the 3rd root of symbolic might not correspond to the third root of whatever you were expecting numerically.

If you do take the symbolic route then I suggest you change the last line of your function from

F = F';

to

F = F(:);

... Though now that I do that myself, I see it switches to finding multiple roots on the vpasolve(), implying that it is fundamentally a polynomial system, with multiple solutions. And that leads to the question of which of the solutions you want. Your inval of [1 2 3] looks a little too contrived for me to expect that you just happen to want the exact solution that is the least distance from [1 2 3].

[Walter Roberson]

There are solutions or near-solutions near

24.62181480998063, 24.246422699, -5.058521441
24.87106907564950, 24.990443679, 2.769055659
25.12105104502035, 24.989669429, 2.698575099
25.39440916065875, 25.813847449, -4.658031121
47.05280359343892, 91.114167519, -1.051418981
49.93676423883585, 99.809343539, 50.910571139
50.06405943862820, 100.191216809, 51.175884399
53.30387805348748, 109.862036629, -1.045739371

The second of those is the one that is closest to [1 2 3], by a small bit.