fsolve function give poor results for multiple equations

6 views (last 30 days)
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]);

Accepted Answer

Matt J
Matt J on 20 Oct 2020
Edited: Matt J on 20 Oct 2020
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,
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
  9 Comments
Matt J
Matt J on 21 Oct 2020
Edited: Matt J on 21 Oct 2020
If you can, upgrade to R2017b
Otherwise,
a=sqrt(sum(A.^2,1)) %vecnorm for older Matlab versions

Sign in to comment.

More Answers (1)

Alan Weiss
Alan Weiss on 20 Oct 2020
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
  1 Comment
Peng Liu
Peng Liu on 21 Oct 2020
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!