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

4 views (last 30 days)
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

Sign in to comment.

Answers (3)

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

Sign in to comment.


John D'Errico
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
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

Sign in to comment.


Bruno Luong
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);
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')
  2 Comments
Bruno Luong
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.

Sign in to comment.

Categories

Find more on Polynomials in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!