fsolve function give poor results for multiple equations

I want to get 7 parameters used in three different equations (let's say equation CP, S, and H), however, I am only happy for the fitted parameters for H equation. The fitted results for CP and S is really poor. How to improve it? What I wish is that the fitted R square is minimun for the three equations. The code is listed below. Thanks in advance.

function y = nonlinmodel(x,Temp,CP,S,H)
y = [CP - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
    S - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
    H - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];
end

Temp=[1200 1300 1500 1600 1700 1800 2000 2200 2500];
S=[520.042428000000 567.730080000000 582.090804000000 595.781640000000 608.844456000000 633.337236000000 655.862220000000 686.593332000000 731.182752000000];
CP=[206.367372000000 220.686228000000 224.245008000000 227.343240000000 230.064660000000 234.586404000000 238.145184000000 242.164512000000 246.644388000000];
H=[-121082.256000000 -56898.6120000000 -34624.8360000000 -12057.9840000000 10801.9440000000 57275.4240000000 104586.264000000 176641.092000000 298937.520000000];
f = @(x)nonlinmodel(x,Temp,CP,S,H);
optimal_x = fsolve(f,[1;0.2;0.03;0.004;0.00005;6;72]);

Answer 1

[Matt J]

Your equations are linear, so there is no reason to be using fsolve. I was able to obtain the coefficient matrix A for your linear system using func2mat from the File Exchange,

https://www.mathworks.com/matlabcentral/fileexchange/44669-func2mat-convert-linear-function-to-matrix

A=full(func2mat(  @(x)CoeffFun(x,Temp,CP,S,H) ,zeros(7,1))),

function y = CoeffFun(x,Temp,CP,S,H)

 Temp=Temp(:);

 y = -[ - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
     - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
     - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];
end

The condition number for you A matrix is very poor, so it's not too surprising that you don't get a good result.

>> cond(A)

ans =

   3.2818e+19

I wonder if Temp is expressed in the correct units? The current units result in an insane order of magnitude for the coefficients,

 >> max(abs(A(:)))

ans =

   1.6238e+17

[Peng Liu]

Dear Matt.

Thanks for the kind reply. The unit of temperatures is right. I try to use Excel to get the coefficients. The fitted results is not that bad, as the attached picture shows. The strategy is to calculate the square difference (like (Cp_exp-Cp_pred.)^2) for each data, then take the sum of all square difference (lets say it as SSE). The next is to use the solve function in excel to get the coefficients a1-a7 via setting the objective SSE to min. What trouble me is that you have to optimize the coefficient one by one, and you may need several loops to get the relative good results.

Therefore, I would like to use matlab to solve this problem. But I do not find the correct way.

[Matt J]

The unit of temperatures is right. I try to use Excel to get the coefficients. The fitted results is not that bad, as the attached picture shows.

The first Excel graph looks noticeably poorer than the others.

Therefore, I would like to use matlab to solve this problem. But I do not find the correct way.

Based on my analysis, it is clear that you will not find one. A condition number of 3.2818e+19 means you will get bad results no matter what method you use.

[Peng Liu]

Dear Matt,

Is it possible to add constrain conditions duiring the optimazation process? What I persue is a globle minimun error. For example, I can accept is the R2 to be 0.9, 0.9, and 0.9 for the fitted three curves, and I can not accept the R2 to be 0.999, 0.998, 0.5.

[Matt J]

You can add more weight to equations that you wish to give more priority, for example,

function yw = nonlinmodel(x,Temp,CP,S,H)

 y = [CP - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
     S - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
     H - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];

 yw=[1000*y(1), 100*y(2), 30*y(3)];

end

Explicit bounds -r(i)<=y(i)<=r(i) can be applied using lsqlin, but I think that will be a less practical approach. You do not know in advance for what combinations of bounds a feasible solution exists.

[Matt J]

You can add more weight to equations that you wish to give more priority

The implementation below follows this idea and gets pretty good results.

Temp=[1200 1300 1500 1600 1700 1800 2000 2200 2500];
S=[520.042428000000 567.730080000000 582.090804000000 595.781640000000 608.844456000000 633.337236000000 655.862220000000 686.593332000000 731.182752000000];
CP=[206.367372000000 220.686228000000 224.245008000000 227.343240000000 230.064660000000 234.586404000000 238.145184000000 242.164512000000 246.644388000000];
H=[-121082.256000000 -56898.6120000000 -34624.8360000000 -12057.9840000000 10801.9440000000 57275.4240000000 104586.264000000 176641.092000000 298937.520000000];

weight=[1e4,1,1];

A=full(func2mat(  @(x)CoeffFun(x,Temp,weight) ,zeros(7,1)));

a=vecnorm(A,2,1);

A=A./a;

b=[CP(:)*weight(1);S(:)*weight(2);H(:)*weight(3)];

x=(A\b)./a(:);

pred=reshape(A*(A\b),[],3)./weight;

%%plot
close all

subplot(3,1,1) 
plot(Temp,pred(:,1),'r',Temp,CP,'bo');
ylabel('CP'), xlabel('Temp')    

subplot(3,1,2)
plot(Temp,pred(:,2),'r',Temp,S,'bo');
ylabel('S'), xlabel('Temp') 

subplot(3,1,3)
plot(Temp,pred(:,3),'r',Temp,H,'bo');
ylabel('H'), xlabel('Temp') 

function y = CoeffFun(x,Temp,weight)

 Temp=Temp(:);

 y = -[ - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4).*weight(1);
     - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7)).*weight(2);
     - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6)).*weight(3)];

 
end

[Peng Liu]

Wow! The curves look very beautiful. You made it, this is what I want.

Thank you again.

When I run the code you paste, I have a error "Undefined function or variable 'func2mat'". Is there anything I missed?

[Peng Liu]

I got the func2mat function from your shared link. However, new error comes "Undefined function or variable 'vecnorm'".

[Matt J]

If you can, upgrade to R2017b

https://www.mathworks.com/help/matlab/ref/vecnorm.html#responsive_offcanvas

Otherwise,

a=sqrt(sum(A.^2,1)) %vecnorm for older Matlab versions

[Peng Liu]

Yeah, It works now. Thank you!

Having a nice day!

Answer 2

[Alan Weiss]

I wasn't able to find a very good answer either. I'm not sure that one exists. I think that lsqnonlin is a more appropriate solver. Here is what I did using MultiStart to search:

Temp=[1200 1300 1500 1600 1700 1800 2000 2200 2500];
S=[520.042428000000 567.730080000000 582.090804000000 595.781640000000 608.844456000000 633.337236000000 655.862220000000 686.593332000000 731.182752000000];
CP=[206.367372000000 220.686228000000 224.245008000000 227.343240000000 230.064660000000 234.586404000000 238.145184000000 242.164512000000 246.644388000000];
H=[-121082.256000000 -56898.6120000000 -34624.8360000000 -12057.9840000000 10801.9440000000 57275.4240000000 104586.264000000 176641.092000000 298937.520000000];
f = @(x)nonlinmodel(x,Temp,CP,S,H);
[optimal_x, resnorm] = lsqnonlin(f,[1;0.2;0.03;0.004;0.00005;6;6]);

lb = [-5000;-2;-.01;-.01;-.01;-10;0];
ub = [2;2;.01;.01;.01;10;20000];
ms = MultiStart;
x0 = lb + rand(size(lb)).*(ub - lb);
prob = createOptimProblem('lsqnonlin','x0',x0,'objective',f);

prob.lb = lb;
prob.ub = ub;
[x,fval,exitflag,output,solutions] = run(ms,prob,1000);

disp(x(2:6))

function y = nonlinmodel(x,Temp,CP,S,H)
y = [CP - 8.314*(x(1)+x(2)*Temp+x(3)*Temp.^2+x(4)*Temp.^3+x(5)*Temp.^4);
    S - 8.314*(x(1)*log(Temp)+x(2)*Temp+x(3)/2*Temp.^2+x(4)/3*Temp.^3+x(5)/4*Temp.^4+x(7));
    H - 8.314*(x(1)*Temp+x(2)*Temp.^2/2+x(3)/3*Temp.^3+x(4)/4*Temp.^4+x(5)/5*Temp.^5+x(6))];
end

The smallest residual I found was near 1e9. Not so good.

Alan Weiss

MATLAB mathematical toolbox documentation

[Peng Liu]

Dear Alan,

I am very appreciate for your feedback. Is it possible to constrain the optimazation by minimizing the sum of squared difference between the estimated values and the actual values? I find that this constrain condition is very useful in Excel for my case.