How to constrain the resulting equation from a polynomial surface fit to a positive range?

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.

[dpb]

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>

[pluie]

I am trying to find an equation that best fits the data.

Answer 1

[Torsten]

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".

[dpb]

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...

Answer 2

[John D'Errico]

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.

[Alex Sha]

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

Answer 3

[Bruno Luong]

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);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

% 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')

[Torsten]

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)

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

P = 9×1

1.0e+04 * 1.388899826775378 0.225628371614422 -0.003037370243114 0.825156361866147 -0.004128526736913 0.053668583711757 -0.000270948473846 -0.000718386473020 0.000003621703902

norm(C*P-d)

ans = 
     5.461779422561789e+05

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)

ans = 
     5.666626680992555e+05

% 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]

+1

the tensorial of 1D polynomial looks better suit than isotropic 2D polynomial basis.