How to constrain the resulting equation from a polynomial surface fit to a positive range?
5 views (last 30 days)
Show older comments
I have a set of data , where . In application, .
clc; clear all; close all;
% minimal working problem
fid = fopen('xyz.txt');
formatspec = ['%f','%f','%f'];
data = textscan(fid,formatspec);
fid = fclose(fid);
x = data{:,1};
y = data{:,2};
z = data{:,3};
sfS = fit([x, y],z,'poly33');
figure;
plot(sfS,[x,y],z)
The resulting polynomial allows for . How do I modify the fit function call to include this constraint? I tried using the Lower option in fitoptions, but that does not seem to be permitted.
Thanks in advance.
2 Comments
dpb
on 4 Sep 2022
Not a viable constraint for a polynomial model -- at least as estimated by fit with OLS.
Let's see what the surface actually looks like and see if can figure out better model...
data=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1116665/xyz.txt');
plot3(data(:,1),data(:,2),data(:,3))
grid on
Well, it's reasonably well behaved excepting those additional peaks are going to create real problems.
What's the purpose in fitting, interpolation, I presume? In that case, I'd recomend forgetting about the surface fit and just use the <scatteredInterpolant>
Answers (3)
Torsten
on 4 Sep 2022
Use "lsqlin" and put the constraints
p33(x(i),y(j)) >= 0
for all (i,j)-combinations.
These constraints can be put as A*coeffs <= b in the matrix A and vector b that you can use in the call to "lsqlin".
1 Comment
dpb
on 4 Sep 2022
I'll have to try to remember the "how" -- about 50 year ago now, Dr. Clark showed me how to introduce a zero boundary condition on a quadratic surface fit we were using in the processing of the reactor incore detector readings to infer power distribution across the rest of the reactor core non-instrumented positions.
It's been so long ago I forget just how he did formulate it, but it worked like a charm there...
John D'Errico
on 4 Sep 2022
The lower bounds in fitoptions apply only to the coefficients of the polynomial. That has NOTHING to do with the value of the function itself.
As far as the use of lsqlin (as suggested by Torsten) this is a good idea, a necessary one, but not a sufficient one. That is, you will be constraining the value of the result at a list of specific points. But that does NOT constrain the polynomial from ever going beyond your goal over the entire domain. It may easily go below zero between the points where you have constrained it.
In fact, most polynomials will be unbounded. The example case, where you used poly33 as the fit type is one that has NO bounds on it. Such a polynomial will go to plus infinity in some direction, and minus infinity in another direction. At best you can hope it obeys your goal over some restricted domain.
1 Comment
Alex Sha
on 5 Sep 2022
Hi, try the function:
z = b0+b1*x+b2*x^2+b3*y+b4*y^2+b5*x*y+b6*x*y^2+b7*x^2*y+b8*x^2*y^2
Ignore the first 19 sets of data,the result will be like below, with all positive values of z
Sum Squared Error (SSE): 344417258173.384
Root of Mean Square Error (RMSE): 44879.1266928162
Correlation Coef. (R): 0.977862014045546
R-Square: 0.956214118513212
Parameter Best Estimate
--------- -------------
b0 29498.3033433469
b1 1534.47381229165
b2 -23.1645304350467
b3 6787.20106399848
b4 -34.2742291974444
b5 597.016772969477
b6 -2.99513376969146
b7 -7.70814745171561
b8 0.0386612650731015
Bruno Luong
on 5 Sep 2022
Edited: Bruno Luong
on 5 Sep 2022
This is implementation of Torsen's idea
data=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1116665/xyz.txt');
x = data(:,1);
y = data(:,2);
z = data(:,3);
% Stuff needed to normalize the data for better inversion
[xmin, xmax] = bounds(x);
[ymin, ymax] = bounds(y);
xynfun = @(x,y)deal((x(:)-xmin)/(xmax-xmin),(y(:)-ymin)/(ymax-ymin));
[xn,yn]=xynfun(x,y);
% cofficients of 2D polynomial 3d order
k = [0 0 1 0 1 2 0 1 2 3];
l = [0 1 0 2 1 0 3 2 1 0];
C = [xn.^k.*yn.^l]; % please no comment about my use of bracket here
d = z;
% Constraint positive of 3 x 3 points in the recatagular domain to be positive,
% it should be enough
[XNC,YNC] = meshgrid(linspace(0,1,3),linspace(0,1,3));
A = -[XNC(:).^k.*YNC(:).^l]; % please no comment ...
b = 0+zeros(size(A,1),1); % A*P<=0 means polynomial at (xnc,ync)>=0
P = lsqlin(C,d,A,b);
% Graphical check
% Create a grided model surface
xi=linspace(xmin,xmax,33);
yi=linspace(ymin,ymax,33);
[Xi,Yi]=meshgrid(xi,yi);
[Xin,Yin]=xynfun(Xi,Yi);
Zi=[Xin.^k.*Yin.^l]*P; % please no comment about my use of bracket here
Zi=reshape(Zi,size(Xi));
close all
surf(xi,yi,Zi);
hold on
plot3(x,y,z,'or')
xlabel('x')
ylabel('y')
zlabel('z')
2 Comments
Torsten
on 5 Sep 2022
Edited: Torsten
on 5 Sep 2022
Compare with Alex Sha's solution which looks promising.
data=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1116665/xyz.txt');
x = data(:,1);
y = data(:,2);
z = data(:,3);
% Stuff needed to normalize the data for better inversion
[xmin, xmax] = bounds(x);
[ymin, ymax] = bounds(y);
xynfun = @(x,y)deal((x(:)-xmin)/(xmax-xmin),(y(:)-ymin)/(ymax-ymin));
%[xn,yn]=xynfun(x,y);
xn = x;
yn = y;
% cofficients of 2D polynomial 3d order
%k = [0 0 1 0 1 2 0 1 2 3];
%l = [0 1 0 2 1 0 3 2 1 0];
k = [0 1 2 0 0 1 1 2 2];
l = [0 0 0 1 2 1 2 1 2];
C = [xn.^k.*yn.^l]; % please no comment about my use of bracket here
d = z;
% Constraint positive of 3 x 3 points in the recatagular domain to be positive,
% it should be enough
%[XNC,YNC] = meshgrid(linspace(0,1,3),linspace(0,1,3));
%[XNC,YNC] = meshgrid(linspace(xmin,xmax,3),linspace(ymin,ymax,3));
%A = -[XNC(:).^k.*YNC(:).^l]; % please no comment ...
A = -[x(:).^k.*y(:).^l];
b = 0+zeros(size(A,1),1); % A*P<=0 means polynomial at (xnc,ync)>=0
format long
P = lsqlin(C,d,A,b)
norm(C*P-d)
b0 = 29498.3033433469 ;
b1 = 1534.47381229165 ;
b2 = -23.1645304350467 ;
b3 = 6787.20106399848 ;
b4 = -34.2742291974444 ;
b5 = 597.016772969477 ;
b6 = -2.99513376969146 ;
b7 = -7.70814745171561 ;
b8 = 0.0386612650731015;
B= [b0;b1;b2;b3;b4;b5;b6;b7;b8];
norm(C*B-d)
% Graphical check
% Create a grided model surface
xi=linspace(xmin,xmax,33);
yi=linspace(ymin,ymax,33);
[Xi,Yi]=meshgrid(xi,yi);
%[Xin,Yin]=xynfun(Xi,Yi);
Xin = Xi(:);
Yin = Yi(:);
Zi=[Xin.^k.*Yin.^l]*P; % please no comment about my use of bracket here
Zi=reshape(Zi,size(Xi));
close all
surf(xi,yi,Zi);
hold on
plot3(x,y,z,'or')
xlabel('x')
ylabel('y')
zlabel('z')
Bruno Luong
on 5 Sep 2022
Edited: Bruno Luong
on 5 Sep 2022
+1
the tensorial of 1D polynomial looks better suit than isotropic 2D polynomial basis.
See Also
Categories
Find more on Polynomials in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!