MATLAB Answers

Adding constraints to lsq fitting

13 ビュー (過去 30 日間)
Benjamin
Benjamin 2019 年 3 月 20 日
編集済み: Benjamin 2019 年 3 月 25 日
I have a code that is of a function g(X,Y). The below code fits my data very well. However, I want to add some constraints at the boundaries where my data is not. For instance, how can I set the coefficient such that when Y=0, then the value of my function is 1?
Can anyone help with this?
clear
clc
close all
num=xlsread('C:data.xlsx',1,'A2:F18110');
eta = num(:,3);
r = num(:,4);
g = num(:,6);
%Do the surface fit
options=optimoptions(@lsqcurvefit,'SpecifyObjectiveGradient',false,'MaxFunctionEvaluations',50000,'MaxIterations',1000);
xdata={r,eta};
params0=linspace(0.1, 0.1, 50);
[params, resnorm]=lsqcurvefit(@modelfun,params0,xdata,g,[],[],options)
for i=1:50
C(i)=params(i);
end
xmin = 1; xmax = 2; dx = 0.01;
etamin=0; etamax=0.55; deta=0.01;
[Xgrid,etagrid]=ndgrid(xmin:dx:xmax,etamin:deta:etamax);
surf(Xgrid,etagrid,modelfun(params,{Xgrid,etagrid}))
zlim([0 8]);
hold on;
scatter3(r(:),eta(:),g(:),'MarkerEdgeColor','none',...
'MarkerFaceColor','b','MarkerFaceAlpha',.5); hold off
xlabel 'x', ylabel '\eta', zlabel 'g(x,\eta)'
function [out,Jacobian]=modelfun(params,xdata)
X=xdata{1};
Y=xdata{2};
for i=1:50
C(i)=params(i);
end
A1 = -0.4194;
A2 = 0.5812;
A3 = 0.6439;
A4 = 0.4730;
eta_c = 0.6452;
PV = 1 + 3*Y./(eta_c-Y)+ A1*(Y./eta_c) + 2*A2*(Y./eta_c).^2 + 3*A3*(Y./eta_c).^3 + 4*A4*(Y./eta_c).^4;
rdf_contact = (PV - 1)./(4*Y);
poly_guess = polyVal2D(C,X-1,Y/eta_c,4,4);
out = (poly_guess.*rdf_contact);
if nargout>1
%Jacobian=[X(:), Y(:), ones(size(X(:)))];
end
end

回答 (3 件)

Adam Danz
Adam Danz 2019 年 3 月 20 日
"...I want to add some constraints at the boundaries... "
In your call to lsqcurvefit(), you're not using the upper and lower bounds options (5th and 6th inputs are empty). I'd start out by imposing boundaries.
"For instance, how can I set the coefficient such that when Y=0, then the value of my function is 1"
I don't follow this part. The variable "Y" doesn't exist in your sample code. Also, what do you mean "the value of my function"? Do you mean the output of your nonlinear function or do you mean the coefficient estimates produced by the fit?
  13 件のコメント
Benjamin
Benjamin 2019 年 3 月 20 日
So then it is the last coefficient that needs to be set to 1 right?

サインインしてコメントする。


Matt J
Matt J 2019 年 3 月 21 日
編集済み: Matt J 2019 年 3 月 21 日
Here is an adaptation of the algebraic solution that I presented in your previous thread. I don't bother with the lsqcurvefit alternative, since as I showed you, it is both slower and less accurate. Below are also surface plots of the fitted surface, both in the neighborhood of the data points and near eta=0 where you can see that everything converges to the target value parameter, g0, as eta-->0. In this example, I have set g0=1.3.
clear
clc
close all
%% Data Set-up
num=xlsread('example.xlsx',1,'A2:F18110');
Np=4; %polynomial order
g0=1.3; %target value at eta=0
g = num(:,6);
r = num(:,4);
eta = num(:,3);
otherStuff.r = r;
otherStuff.eta = eta;
otherStuff.g = g;
otherStuff.eta_c = 0.6452;
otherStuff.Np=Np;
otherStuff.A1 = -0.4194;
otherStuff.A2 = 0.5812;
otherStuff.A3 = 0.6439;
otherStuff.A4 = 0.4730;
%% Fit using matrix algebra
A0=func2mat(@(p) modelfun(p,otherStuff), ones(Np+1,Np+1));
q=4*g0./((3+otherStuff.A1)/otherStuff.eta_c); %weight on limiting target value as eta-->0
b=g(:)-q*A0(:,1);
A=A0;
A(:,1:Np+1)=[];
coeffs =[zeros(Np+1,1); A\b]; %fitted coefficients
coeffs(1)=q;
figure(1); showFit(coeffs,otherStuff);
%% Also check surface near eta=0
[Rg,Eg]=ndgrid(unique(r),linspace(1e-8,0.01,100));
differentStuff=otherStuff;
differentStuff.r = Rg(:);
differentStuff.eta = Eg(:);
differentStuff.g=nan([numel(Rg),1]);
figure(2); showFit(coeffs,differentStuff);
function ghat=modelfun(C,otherStuff)
r = otherStuff.r;
eta = otherStuff.eta;
eta_c = otherStuff.eta_c;
Np = otherStuff.Np;
A1 = otherStuff.A1;
A2 = otherStuff.A2;
A3 = otherStuff.A3;
A4 = otherStuff.A4;
PV = 1 + 3*eta./(eta_c-eta)+ A1*(eta./eta_c) + 2*A2*(eta./eta_c).^2 +...
3*A3*(eta./eta_c).^3 + 4*A4*(eta./eta_c).^4;
rdf_contact = (PV - 1)./(4*eta);
Cflip=rot90(reshape(C,Np+1,Np+1),2);
poly_guess = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);
ghat = (poly_guess.*rdf_contact);
end
function showFit(coeffs,otherStuff)
r = otherStuff.r(:);
eta = otherStuff.eta(:);
g = otherStuff.g(:);
ghat=modelfun(coeffs,otherStuff);
tri = delaunay(r,eta);
%% Plot it with TRISURF
h=trisurf(tri, r, eta, ghat);
h.EdgeColor='b';
h.FaceColor='b';
axis vis3d
hold on;
scatter3(r,eta,g,'MarkerEdgeColor','none',...
'MarkerFaceColor','r','MarkerFaceAlpha',.05);
xlabel 'r', ylabel '\eta', zlabel 'g(r,\eta)'
hold off
end
  7 件のコメント
Benjamin
Benjamin 2019 年 3 月 21 日
So, I've changed the excel file. I now have already divided out rdf_contact. Here is the new excel file. The new g is column H. It already had divided out rdf_contact so I am just fitting the polynomial.
Any chance you could modify the code without the weighting function and with the new data? Note that this means rdf_contact and all the coeffs (A1, A2, ..) with it can be removed.
I apologize if I should have done this from the beginning.
When I use cftool on (r,eta,g), I feel like the fits don't look as good...It actually looks a lot worse than what you were able to do in your code. Unless I am missing something. I also don't feel like I have as much control over changing things.
Could I also add data points at eta =0 , that g=1 for all x? That might help the fit too...I keep trying wtih cftool, and the fits dont look as good as in your matrix algebra code. Any chance you could change the code with the new things? I still can use the old code for when I have rdf_contact. But I want to try it without.
So to sum up, new excel file, column H has new g data with rdf_contact already accounted for. Rdf_contact can be removed from matlab code. And I want to set:
So to be clear. I need (where x = i, eta = j)
Cij format
C00 = 1;
Ci0's (other than the one listed above) = 0
C01 = 8/eta_c
C11 = -6/eta_c
C31 = 1/2/eta_c
Ci1's (that are not listed above) = 0
So it would be good to remove weighting function I think and just set these coeffs. Then I should be good to go!
Edit: I keep trying to get good fits in Cftool... I think your matrix algebra is better. Any chance you could update the code to ignore the rdf_contact, just read in column h for g(x,eta), and set the above coefficients?
Edit: When I change your code, my surface fit is completely messed up. I'm not sure if I'm not putting the coefficients in correctly?
Does coeff(9) correspond to C31?
Also with these coeffs, r-1 needs to be just r.
Sorry for all the comments...
Also, you pass p to model fun. What is p? It's not defined anywhere that I can tell?

サインインしてコメントする。


Matt J
Matt J 2019 年 3 月 21 日
Any chance you could modify the code without the weighting function and with the new data? Note that this means rdf_contact and all the coeffs (A1, A2, ..) with it can be removed.
It's virtually trivial now, with the Curve Fitting Toolbox:
%% Load Data
num=xlsread('example2.xlsx',1,'A2:F18110');
eta_c= 0.6452;
r = num(:,4);
eta = num(:,3);
H = num(:,5);
%% Set-up for fit
[I,J]=ndgrid(0:4);
Terms= compose('C%u%u*r^%u*eta^%u', I(:),J(:),I(:),J(:)) ;
Terms=strjoin(Terms,' + ');
independent={'r','eta'};
dependent='H';
knownCoeffs= {'C00','C10','C20','C30','C40', 'C01','C11','C21','C31','C41'};
knownVals=num2cell([ [1 , 0 , 0 , 0 , 0 ], [ 8 , -6 , 0 , 0.5 , 0 ]/eta_c ]);
allCoeffs=compose('C%u%u',I(:),J(:));
[unknownCoeffs,include]=setdiff(allCoeffs,knownCoeffs);
ft=fittype(Terms,'independent',independent, 'dependent', dependent,...
'coefficients', unknownCoeffs,'problem', knownCoeffs);
%% Fit and display
fobj = fit([r,eta],H,ft,'problem',knownVals) ;
hp=plot(fobj,[r,eta],H);
set(hp(1),'FaceAlpha',0.5);
set(hp(2),'MarkerFaceColor','r');
xlabel 'r',
ylabel '\eta'
zlabel 'H(r,\eta)'
untitled.png
  4 件のコメント
Benjamin
Benjamin 2019 年 3 月 25 日
If you notice, the surface still does not extend to eta=0. How do I get it to do that? Currently, the surface is only plotted over the data. I would like the surface to go all the way to eta/eta_c = 0.
Edit: I also created a new question that is related to this in another thread.

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by