Solve system of simultaneous equations for only real numbers

Hello,

I am trying to solve the following system of equations and while I do get an answer, this is complex. For my application, only positive & real solutions are relevant. Furthermore, the solutions to "c1", "c2", "c4" should lie in the range 0 to 1. And c1 + c2 + c3 + c4 should sum to 1.

Is there a way to use vpasolve to find solutions for c1, c2, c4 that are positive, real, and between 0 to 1, or can anyone suggest another solver to use for my problem please?

Thank you.

Code:

clear all;

clc;

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];

c3 = 0.95;

syms n c1 c2 c4

init_param = [0; 1];

eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );

eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );

eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );

eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );

sols = vpasolve([eqn1 eqn2 eqn3 eqn4], [n; c1; c2; c4], init_param)

Gives:

Error using sym/vpasolve>checkMultiVarX

Incompatible starting points and variables.

Code:

clear all;

clc;

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];

c3 = 0.95;

syms n c1 c2 c4

init_param = [0; 1];

eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );

eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );

eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );

eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );

sols = vpasolve([eqn1 eqn2 eqn3 eqn4], [n; c1; c2; c4], init_param)

Gives:

n: 6.7426863581108857304747603024111 + 3.9194661375046172368785216288436i

c1: 0.030720057908933124309712773265263 + 0.027689129624897200337074934618641i

c2: 0.0022216968838103768366405641539907 + 0.0018149441042006610883472427462662i

c4: 0.017058245207256498853646662580747 - 0.029504073729097861425422177364907i

[Dyuman Joshi]

Do you know if a real solution exists for all variables?

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;
syms n c1 c2 c4
%4 variablees, need 4 set of initial parameters
init_param = repmat([0 1],4,1)

init_param = 4×2

0 1 0 1 0 1 0 1

eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );
eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );
eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );
eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );
sols = vpasolve([eqn1; eqn2; eqn3; eqn4], [n; c1; c2; c4], init_param)

sols = struct with fields:
     n: [0×1 sym]
    c1: [0×1 sym]
    c2: [0×1 sym]
    c4: [0×1 sym]

[Grass]

Hi @Dyuman Joshi, thank you for your answer. I believe that a real solution does exist for all variables... Please see my comment on Fabio's answer. Thank you for the help.

Answer 1

[Fabio Freschi]

I tried with a numerical solution

clear all;

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;

% variable mapping
% n -> x(1)
% c1 -> x(2)
% c2 -> x(3)
% c4 -> x(4)

fun = @(x)[log10((x(2)/a(1)) * a(4)/x(4)) / log10(1.57488) - x(1);
    log10((x(3)/a(2)) * a(4)/x(4)) / log10(1.55548) - x(1);
    log10((c3/a(3)) * a(4)/x(4)) / log10(1.30607) - x(1);
    x(2)+x(3)+x(4)+c3-1];

options = optimoptions('fsolve','MaxFunctionEvaluations',10000,'MaxIterations',10000,...
    'FunctionTolerance',1e-10);

x = fsolve(fun,[1 1 1 1],options);

The solution has only positive values

>> x

x =

    6.9618    0.0431    0.0030    0.0321

however the solution is not very accourate

>> fun(x)

ans =

   -0.0006
   -0.0000
    0.0006
    0.0282
    

Are you sure about the use of parenthesis in your original formulation?

[Dyuman Joshi]

You skipped over an important detail -

fsolve did not find a solution, it returned the last values before terminating the solver.

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;

% variable mapping
% n -> x(1)
% c1 -> x(2)
% c2 -> x(3)
% c4 -> x(4)

fun = @(x)[log10((x(2)/a(1)) * a(4)/x(4)) / log10(1.57488) - x(1);
    log10((x(3)/a(2)) * a(4)/x(4)) / log10(1.55548) - x(1);
    log10((c3/a(3)) * a(4)/x(4)) / log10(1.30607) - x(1);
    x(2)+x(3)+x(4)+c3-1];

options = optimoptions('fsolve','MaxFunctionEvaluations',10000,'MaxIterations',10000,...
    'FunctionTolerance',1e-10);

x = fsolve(fun,[1 1 1 1],options)

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

x = 1×4

6.9618 0.0431 0.0030 0.0321

[Fabio Freschi]

Yes but no. I showed fun(x) that reports how close the solution is to the desired result. This is why I asked if equations are correct

[Grass]

Hi @Fabio Freschi, @Dyuman Joshi, thank you for showing your approaches to the problem. The equation I am trying to solve is Fenske's equation, commonly used for distillation design. In my code, n = Nmin, c1 = xDi1, c2 = xDi2, c3 = xDi3, c4 = xDi4. The minimum number of stages should be identical no matter which component is used in the equation, as long as the xBr & xDr are the same for all equations & the molar composition (x-values) is constant. The last relationship is from continuity: all xDi - values should sum to 1.

I must admit however that my MATLAB skills are not very well refined as of yet.

The problem with the solution: Nmin = 6.9618, { xDc1 = 0.0431, xDc2 = 0.030, xDc4 = 0.0321 } is that 0.0431 + 0.030 + 0.0321 ≠ 0.05. @Dyuman Joshi @Fabio Freschi, Do you think there's a way to satisfy this condition?

Thanks again.

[Dyuman Joshi]

The use of parenthesis was incorrect, as Fabio mentioned.

Fabio used the correct formula and the result obtained did not satisfy the constraints. The result obtained does not qualify as a solution as the message from fsolve indicates.

I have corrected it the formula in the original code and ran the code, still no solution was obtained. I have changed the base of log to natural base, as in the formula.

Are you sure all the numerical values are correct? I would suggest you to cross-check once again.

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;
syms n c1 c2 c4
%Starting points instead of range
init_param = zeros(4,1);
eqn1 = n == log((c1/a(1,1)) * (a(4,1)/c4)) / log(1.57488);
eqn2 = n == log((c2/a(2,1)) * (a(4,1)/c4)) / log(1.55548);
eqn3 = n == log((c3/a(3,1)) * (a(4,1)/c4)) / log(1.30607);
eqn4 = 1 - (c1 + c2 + c3 + c4) == 0;
sols = vpasolve([eqn1; eqn2; eqn3; eqn4], [n; c1; c2; c4], init_param)

sols = struct with fields:
     n: [0×1 sym]
    c1: [0×1 sym]
    c2: [0×1 sym]
    c4: [0×1 sym]

[Grass]

Hi @Dyuman Joshi, thank you for continuing to assist me with this. With regards to the log base 10 vs log base e, I don't think this is the issue e.g. log10(3) / log10(4) == log30(3) / log30(4) == loge(3) / loge(4) etc.

The values in "a" & {1.57488, 1.55548, 1.30607} are the mole fractions in the bottoms & relative volatilities, respectively - these I think are correct. The "0.95" is the desired mole fraction of the product (which I "chose").

I'm really stuck with this problem... :( - are there any other solvers that could be used other than vpasolve or fsolve? I know that quite roughly: n = 7, c1 = 0.02, c2 = 0.02, c4 = 0.01.

[Fabio Freschi]

Maybe your choice for c3 is too demanding. In fact, if you use my code and set c3 = 0.9, fsolve exit with a good solution

clear all;
clc;

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.9;

fun = @(x)[log10((x(2)/a(1)) * a(4)/x(4)) / log10(1.57488) - x(1);
    log10((x(3)/a(2)) * a(4)/x(4)) / log10(1.55548) - x(1);
    log10((c3/a(3)) * a(4)/x(4)) / log10(1.30607) - x(1);
    x(2)+x(3)+x(4)+c3-1];

options = optimoptions('fsolve','MaxFunctionEvaluations',10000,'MaxIterations',10000,...
    'FunctionTolerance',1e-10);

[x,fval] = fsolve(fun,[7 .02 .02 .01],options)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×4

10.7873 0.0835 0.0055 0.0110

fval = 4×1

1.0e-12 * -0.2700 -0.3446 -0.4459 0

Note that I also used your initial guess, that usually helps in such problems.

With c3 = 0.9, also vpasolve is successful

a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.9;

syms n c1 c2 c4

eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );
eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );
eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );
eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );

sols = vpasolve([eqn1 eqn2 eqn3 eqn4], [n; c1; c2; c4],[6.9618    0.0431    0.0030    0.0321])

sols = struct with fields:
     n: 10.787264649094921719080297249444
    c1: 0.083524988284633557898145278973197
    c2: 0.0055105259988317881219826499549446
    c4: 0.010964485716534653979872071071858

[Alex Sha]

if taking c3=0.9, there will be one more solution:

n	3.47821500649581
c1	0.021268137459732
c2	0.00153621135446466
c4	0.0771956511857335

if taking c3=0.92, there will be also two solution:

No.	1	                2
n	5.49455717047147	8.48765845854466
c1	0.031707727103106	0.055519739381755
c2	0.00223373979164014	0.00376879846933087
c4	0.0460585331052497	0.0207114621489085

[Grass]

Hi @Alex Sha, thanks for this, may I ask how you arrived at these two solutions please? Did you use the same code? Thanks

Answer 2

[Grass]

@Dyuman Joshi @Fabio Freschi . I believe you were correct in saying one of the numeric values was incorrect - it seems that if c3 = 0.92 instead of 0.95 - a real, positive soution is found as required. I think I was being too ambitious with what molar composition I could achieve with the distillation. Thanks for the help!