Adding constraints to lsq fitting

6 views (last 30 days)
Benjamin on 20 Mar 2019
Edited: Benjamin on 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?
close all
eta = num(:,3);
r = num(:,4);
g = num(:,6);
%Do the surface fit
params0=linspace(0.1, 0.1, 50);
[params, resnorm]=lsqcurvefit(@modelfun,params0,xdata,g,[],[],options)
for i=1:50
xmin = 1; xmax = 2; dx = 0.01;
etamin=0; etamax=0.55; deta=0.01;
zlim([0 8]);
hold on;
'MarkerFaceColor','b','MarkerFaceAlpha',.5); hold off
xlabel 'x', ylabel '\eta', zlabel 'g(x,\eta)'
function [out,Jacobian]=modelfun(params,xdata)
for i=1:50
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(:)))];

Answers (3)

Adam Danz
Adam Danz on 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?
Benjamin on 20 Mar 2019
So then it is the last coefficient that needs to be set to 1 right?

Sign in to comment.

Matt J
Matt J on 21 Mar 2019
Edited: Matt J on 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.
close all
%% Data Set-up
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.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
coeffs =[zeros(Np+1,1); A\b]; %fitted coefficients
figure(1); showFit(coeffs,otherStuff);
%% Also check surface near eta=0
differentStuff.r = Rg(:);
differentStuff.eta = Eg(:);
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);
poly_guess = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);
ghat = (poly_guess.*rdf_contact);
function showFit(coeffs,otherStuff)
r = otherStuff.r(:);
eta = otherStuff.eta(:);
g = otherStuff.g(:);
tri = delaunay(r,eta);
%% Plot it with TRISURF
h=trisurf(tri, r, eta, ghat);
axis vis3d
hold on;
xlabel 'r', ylabel '\eta', zlabel 'g(r,\eta)'
hold off
Benjamin on 21 Mar 2019
Edited: Benjamin on 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?

Sign in to comment.

Matt J
Matt J on 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
eta_c= 0.6452;
r = num(:,4);
eta = num(:,3);
H = num(:,5);
%% Set-up for fit
Terms= compose('C%u%u*r^%u*eta^%u', I(:),J(:),I(:),J(:)) ;
Terms=strjoin(Terms,' + ');
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 ]);
ft=fittype(Terms,'independent',independent, 'dependent', dependent,...
'coefficients', unknownCoeffs,'problem', knownCoeffs);
%% Fit and display
fobj = fit([r,eta],H,ft,'problem',knownVals) ;
xlabel 'r',
ylabel '\eta'
zlabel 'H(r,\eta)'
Benjamin on 25 Mar 2019
Edited: Benjamin on 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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!