Computation of distance of a set of points from a fitting curve (of 2n grade like a parabola).

Question

Giuseppe Zumbo 2024-1-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2076541-computation-of-distance-of-a-set-of-points-from-a-fitting-curve-of-2n-grade-like-a-parabola

编辑： Torsten 2024-1-31

Hi everybody, It seems to be an easy problem but I m in trouble:

I have to compute the distance btw the black points (distance for each of these points)( which have vector coordinates (flowrate, pPHeadLoss) ) from the red curve that is the fitting curve which has as coordinates the row vectors ((min(flowrate):0.0001:max(flowrate), y)

Actually I am truing to use pdist function but it doesnt work for me, what am I wrong?

Thanks in advance

[y,delta] = polyval(p, (min(flowrate):0.0001:max(flowrate)),S,mu);

% FIT

plot(min(flowrate):0.0001:max(flowrate), y , 'r'); % red curve

The full code is the next one:

flowdata = readtable("flodata.txt");
% column 1: pressure at pipe start [bars]
% column 2: acoustic noise  at start [standard deviation]
% column 3: pressure at pipe end (10km from start) [bars]
% column 4: acoustic noise  at pipe end[standard deviation]
% column 5: flow rate [m3/hour]
load('flowdata.mat');
% summary(flowdata)
% Vectors of Pressures
pressureIN = flowdata.pressureAtPipeStart_bars_;
pressureOUT= flowdata.pressureAtPipeEnd_10kmFromStart__bars_;
% Vectors of Acoustic Noises at IN and at OUT
aNoiseIN=flowdata.acousticNoiseAtStart_standardDeviation_;
aNoiseOUT=flowdata.acousticNoiseAtPipeEnd_standardDeviation_;
% Vector of flowrate
flowrate = flowdata.flowRate_m3_hour_;
% % Table of flowRate
% flowRate = table(flowrate);
% rawData = [pressureIN aNoiseIN pressureOUT aNoiseOUT flowrate];
% t = (linspace(0, length(cleanData)*10, length(cleanData)))/86400;
sample_distance = 10;
%t = (linspace(0, 7948790, 794879))/86400;
[samples,channels] = size(flowdata);
t = linspace(0, samples*sample_distance, samples)/86400;
% all main curves overlaid 
figure()
plot(t, pressureIN, '-g',t, pressureOUT, '-r',t, flowrate, '-b');
axis tight;
title(' Input , Output Pressure , FR vs time');
legend('Input Pressure','Output Pressure','Flow Rate');
xlabel('Tempo [Days]');
ylabel('all');
PHeadLoss = pressureIN-pressureOUT;
densityplot(flowrate, PHeadLoss, [100,100])
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
rawData = [pressureIN pressureOUT flowrate];
% TOLGO I NaN
% Creo un array di indici logici che indica quali righe contengono almeno un NaN
idx_nan = any(isnan(rawData), 2);
cleanData = rawData;
% Elimino le righe selezionate
cleanData(idx_nan, :) = [];
[samples,channels] = size(cleanData);
t = linspace(0, samples*sample_distance, samples)/86400;
% Vectors of Pressures
pressureIN = cleanData(:,1);
pressureOUT= cleanData(:,2);
flowrate = cleanData(:,3);
PHeadLoss = pressureIN - pressureOUT;
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
hold on
%densityplot(flowrate, PHeadLoss, [50,50]);
[p, S, mu] = polyfit(flowrate, PHeadLoss, 2)
[y,delta] = polyval(p,flowrate,S,mu);
%LINEAR FIT
plot(flowrate, y, 'r-');
%hold on
% Linear Fit of Data with 95% Prediction Interval
%plot(flowrate,y+2*delta,'m--',flowrate,y-2*delta,'m--')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
densityplot(flowrate, PHeadLoss, [100,100])
% TOLGO I VALORI NEGATIVI
% Creo un array di indici logici che indica quali righe hanno almeno un valore negativo
idx_neg = any(cleanData < 0, 2);
cleanData = cleanData;
% Elimino le righe selezionate
cleanData(idx_neg, :) = [];
[samples,channels] = size(cleanData);
t = linspace(0, samples*sample_distance, samples)/86400;
pressureIN = cleanData(:,1);
pressureOUT= cleanData(:,2);
flowrate = cleanData(:,3);
PHeadLoss = pressureIN - pressureOUT;
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
hold on
%densityplot(flowrate, PHeadLoss, [50,50]);
[p, S, mu] = polyfit(flowrate, PHeadLoss, 2)
[y,delta] = polyval(p,flowrate,S,mu);
%LINEAR FIT
plot(flowrate, y, 'r-');
%hold on
% Linear Fit of Data with 95% Prediction Interval
%plot(flowrate,y+2*delta,'m--',flowrate,y-2*delta,'m--')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
densityplot(flowrate, PHeadLoss, [100,100])
% % Table of flowRate
% flowRate = table(flowrate);
% rawData = [pressureIN aNoiseIN pressureOUT aNoiseOUT flowrate];
% t = (linspace(0, length(cleanData)*10, length(cleanData)))/86400;
sample_distance = 10;
%t = (linspace(0, 7948790, 794879))/86400;
[samples,channels] = size(cleanData);
t = linspace(0, samples*sample_distance, samples)/86400;
% all main curves overlaid 
figure()
plot(t, pressureIN, '-g',t, pressureOUT, '-r',t, flowrate, '-b');
axis tight;
title(' Input , Output Pressure , FR vs time');
legend('Input Pressure','Output Pressure','Flow Rate');
xlabel('Tempo [Days]');
ylabel('all');
% Prima considerazione
% large variations on the output pressure can be seen at days 1 to 15 , 
% not related to input pression
% we should not use that portion of the data 
% It's a brutal cleaning but... 
% I am sceptical about what happens during the first 15 days : 
% lots of peaks are seen on the output pressure side that are 
% not reflected / correlated to the input side , 
% so I prefer to discard the first 20 days (yes it's brutal !)
ind = t>13;
t = t(ind);
pressureIN = pressureIN(ind);
pressureOUT = pressureOUT(ind);
flowrate = flowrate(ind);
PHeadLoss = pressureIN-pressureOUT;
figure()
plot(t, pressureIN, '-g',t, pressureOUT, '-r',t, flowrate, '-b');
axis tight;
title(' Input , Output Pressure , FR vs time');
legend('Input Pressure','Output Pressure','Flow Rate');
xlabel('Tempo [Days]');
ylabel('all');
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
hold on
%densityplot(flowrate, PHeadLoss, [50,50]);
[p, S, mu] = polyfit(flowrate, PHeadLoss, 2)
[y,delta] = polyval(p,flowrate,S,mu);
%LINEAR FIT
plot(flowrate, y, 'r-');
%hold on
% Linear Fit of Data with 95% Prediction Interval
%plot(flowrate,y+2*delta,'m--',flowrate,y-2*delta,'m--')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
densityplot(flowrate, y, [100,100]);
% seconda considerazione
% keep data (valid) for gradient below a certain threshold (and non zero
% FR too btw)
% i also don't like those sudden drops to zero , 
% because to me its like dividing zero by zero , 
% when you want to make a relationship between two variables ,
% I avoid using background noise data , it will blurr the result.
% This function gradient calculates the gradient of a signal, 
% which is the rate of change of the signal over time.
pressureIN_grad = abs(gradient(pressureIN));
pressureOUT_grad = abs(gradient(pressureOUT));
flowrate_grad = abs(gradient(flowrate));
threshold = 0.0073; 
ind1 = pressureIN_grad<threshold*max(pressureIN_grad);
ind2 = pressureOUT_grad<threshold*max(pressureOUT_grad);
ind3 = flowrate_grad<threshold*max(flowrate_grad);
ind4 =  flowrate> 0; 
ind = ind1 & ind2 & ind3  & ind4;
t = t(ind);
pressureIN = pressureIN(ind);
pressureOUT = pressureOUT(ind);
flowrate = flowrate(ind);
figure()
plot(t, pressureIN, '.g',t, pressureOUT, '.r',t, flowrate, '.b');
axis tight;
title(' Input , Output Pressure , FR vs time');
legend('Input Pressure','Output Pressure','Flow Rate');
xlabel('Tempo [Days]');
ylabel('all');
PHeadLoss = pressureIN - pressureOUT;
densityplot(flowrate, PHeadLoss, [50,50])
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
hold on
%densityplot(flowrate, PHeadLoss, [50,50]);
[p, S, mu] = polyfit(flowrate, PHeadLoss, 2)
[y,delta] = polyval(p,flowrate,S,mu);
%LINEAR FIT
plot(flowrate, y, 'r-');
%hold on
% Linear Fit of Data with 95% Prediction Interval
%plot(flowrate,y+2*delta,'m--',flowrate,y-2*delta,'m--')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
densityplot(flowrate, y, [50,50]);
% terza considerazione
% PRESSURE HEAD LOSS (APPROX _no fattore di attrito, flusso del fluido e
% densità del fluido)
PHeadLoss = pressureIN - pressureOUT;
% consider only positive values 
ind =  PHeadLoss> 0; 
t = t(ind);
pressureIN = pressureIN(ind);
pressureOUT = pressureOUT(ind);
flowrate = flowrate(ind);
PHeadLoss = pressureIN - pressureOUT;
densityplot(flowrate, PHeadLoss, [50,50])
ylabel('P Head Loss');
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
plot(flowrate, PHeadLoss, '*k');
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
hold on
%densityplot(flowrate, PHeadLoss, [50,50]);
[p, S, mu] = polyfit(flowrate, PHeadLoss, 2)
n=size(flowrate)
[y,delta] = polyval(p, (min(flowrate):0.0001:max(flowrate)),S,mu);
% FIT
plot(min(flowrate):0.0001:max(flowrate), y , 'r');
%hold on
% Linear Fit of Data with 95% Prediction Interval
%plot(flowrate,y+2*delta,'m--',flowrate,y-2*delta,'m--')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
ylabel('P Head Loss');
densityplot(flowrate, y, [50,50]);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-1-31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2076541-computation-of-distance-of-a-set-of-points-from-a-fitting-curve-of-2n-grade-like-a-parabola#answer_1400516

移动：Torsten 2024-1-31

Use "polyval" to compute the polynomial values p at the flow rates f of the "black points".

The distance is then abs(p-Head Loss of the black points).

2 个评论
显示无隐藏无

Giuseppe Zumbo 2024-1-31

Hi @Torsten, thanks for answeing but in this way this distance would not be the minimum btw the point of the dataset and the curve, but it's just the vertical distance, isn't it ?

To estimate the accuracy, DO I need this one or the minimum distance? THanks in advance

Torsten 2024-1-31

编辑：Torsten 2024-1-31

"Polyfit" minimizes the sum of vertical distances squared between your measurement data and the polynomial curve. So you also should use the vertical distance in your evaluation, shouldn't you ?

请先登录，再进行评论。