Solve nonlinear objective function with linear equality constrain (Chemical composition case)

I have a problem to solve this equation. I have 9 equations with 9 variables. Three equations above is a linear constrain, I get this equation from balance atom for chemical composition. Six equations below is the gibbs energy, as the objective function. So, I used MatLab R2013a, and for the first time I use fsolve. But, I get the wrong result when I validate the output. Am I do wrong? Should I change the function using fmincon? what is the better one? Please help me. Thankyou.

I have two m.files which I use. The first m.files is named "cobafun.m" where I state my inputs, and the second one is named "gasf_func.m" where I use fsolve.

.......................Here is my cobafun.m file....................

function [fx] = cobafun(x)

nH2 = x(1);

nC = x(2);

nCO2 = x(3);

nH2O = x(4);

nCO = x(5);

nCH4 = x(6);

lamda1 = x(7);

lamda2 = x(8);

lamda3 = x(9);

R = 0.008314; %kJ/(mol.K)

%that will be variated

a = 0.5; %molar ratio

T = 773; %Temperature

%total of mol

ntotal = nH2 + nC + nCO2 + nH2O + nCO + nCH4;

%balance atom equation

fx (7,1) = 2*nH2 + 4*nCH4 + 2*nH2O - 1.67 + 2*a; %hydrogen (H) balance

fx (8,1) = nCO2 + nCO + nC + nCH4 - 1; %carbon (C) balance

fx (9,1) = 2*nH2 + nCO + nH2O - 0.75 - a; %oksygen (O) balance

%Gibbs equation

fx (1,1) = log(nH2/ntotal)/R*T + 2*lamda1;

fx (2,1) = log(nC/ntotal)/R*T + lamda2;

fx (3,1) = -395.4928817/R*T + log(nCO2/ntotal)/R*T + lamda2 + 2*lamda3;

fx (4,1) = -205.0850656/R*T + log(nH2O/ntotal)/R*T + 2*lamda1 + lamda3;

fx (5,1) = -180.3948228/R*T + log(nCO/ntotal)/R*T + lamda2 + lamda3;

fx (6,1) = -4.922967331/R*T + log(nCH4/ntotal)/R*T + lamda2 + 4*lamda1;

.....................and there is my gasf_func.m file.......................

function gasf_func

clear, clc;

x0 = [0.01 0.01 0.01 0.01 0.01 0.01 1 1 1]; %initial guess (still confuse how to take the value)

disp('variable values at the initial estimate');

f0=cobafun(x0);

disp(' variabel value function value');

for i=1:size(x0,2);

    disp([' x' int2str(i) '                 ' num2str(x0(i)) '             ' num2str(f0(i))])

end

options=optimset('Display','off','MaxIter',10000,'MaxFunEvals',100000);

xsolv=fsolve(@cobafun,x0,options);

disp('variable at the solution');

fsolv=cobafun(real(xsolv));

disp(' variable value function value');

for i=1:size(x0,2);

    disp([' x' int2str(i) '                 ' num2str(real(xsolv(i))) '             ' num2str(fsolv(i))])

end

........................the output......................

variable values at the initial estimate

 variabel              value              function value

 x1                 0.01             -166588.0974

 x2                 0.01             -166589.0974

 x3                 0.01             -36937816.0551

 x4                 0.01             -19234515.3761

 x5                 0.01             -16938923.6789

 x6                 0.01             -624301.449

 x7                 1             -0.59

 x8                 1             -0.96

 x9                 1             -1.21

variable at the solution

 variable               value              function value

 x1                 -0.36624             2519150.0516

 x2                 -0.36624             2519149.044

 x3                 0.73908             -34186797.6308+292091.787494i

 x4                 0.20735             -16601669.0831+292091.787494i

 x5                 0.13854             -14343569.6736+292091.787494i

 x6                 -0.35249             2057879.9846

 x7                 1             -2.3978

 x8                 1             -0.84111

 x9                 1             -1.6366

>>

____________________________________________________________________________________

[Alireza Kakoee]

did you find any answer for your gibbs function?

[Alireza Kakoee]

did you find any answer for your gibbs function?

Answer 1

[Alireza Kakoee]

theren is a problem in your gibbs function, I am working on gibbs minimazation, did you have any gibbs fuction ?