Adding constraints to lsq fitting

8 vues (au cours des 30 derniers jours)
Benjamin
Benjamin le 20 Mar 2019
Modifié(e) : Benjamin le 25 Mar 2019
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

Réponses (3)

Adam Danz
Adam Danz le 20 Mar 2019
"...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 commentaires
Adam Danz
Adam Danz le 20 Mar 2019
Assuming I've found the correct polyVal2D() function you're using,
% Polynomial coefficients are in the following order.
%
% F(X,Y) = P_1 * X^N * Y^M + P_2 * X^{N-1} * Y^M + ... + P_{N+1} * Y^M + ...
% P_{N+2} * X^N * Y^{M-1} + P_{N+3} * X^{N-1} * Y^{M-1} + ... + P_{2*(N+1)} * Y^{M-1} + ...
% ...
% P_{M*(N+1)+1} * X^N + P_{M*(N+1)+2} * X^{N-1} + ... + P_{(N+1)*(M+1)}
Benjamin
Benjamin le 20 Mar 2019
So then it is the last coefficient that needs to be set to 1 right?

Connectez-vous pour commenter.


Matt J
Matt J le 21 Mar 2019
Modifié(e) : Matt J le 21 Mar 2019
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 commentaires
Matt J
Matt J le 21 Mar 2019
Can you remove the weighting function, and then I can just set my coeffs to what they need to be?
It's your model. I assume rdf_contact has a reason to be there. As you'll recall, that was the whole reason you couldn't simply do this with straight up polynomial fitting with the Curve Fitting Toolbox.
Do I add like r = [1 2] and eta=[0 1]?
You add whatever additional points you want included in the surface plot.
I have no idea where this came from:
By setting C(1,1) to some constant q, as you've been proposing, and C(2:Np+1,1)=0 and letting eta-->0 you will find that the surface approaches the limiting value
LimitingValue=(3+otherStuff.A1)/otherStuff.eta_c/4 * q;
Therefore, setting
q=4.*g0./((3+otherStuff.A1)/otherStuff.eta_c));
leads to
LimitingValue = g0
Benjamin
Benjamin le 21 Mar 2019
Modifié(e) : Benjamin le 21 Mar 2019
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?

Connectez-vous pour commenter.


Matt J
Matt J le 21 Mar 2019
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 commentaires
Matt J
Matt J le 22 Mar 2019
The fit that I get looks pretty good to me...
untitled.png
Benjamin
Benjamin le 25 Mar 2019
Modifié(e) : Benjamin le 25 Mar 2019
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.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Get Started with Curve Fitting Toolbox dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by