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

1 visualizzazione (ultimi 30 giorni)
Giuseppe Zumbo il 31 Gen 2024
Modificato: Torsten il 31 Gen 2024
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?
[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:
% 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]
% 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');
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
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);
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
hold on
[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');
% 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);
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
hold on
[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');
% % 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);
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')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
hold on
[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');
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.
% which is the rate of change of the signal over time.
threshold = 0.0073;
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');
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
hold on
[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');
densityplot(flowrate, y, [50,50]);
% terza considerazione
% PRESSURE HEAD LOSS (APPROX _no fattore di attrito, flusso del fluido e
% densità del fluido)
% consider only positive values
t = t(ind);
pressureIN = pressureIN(ind);
pressureOUT = pressureOUT(ind);
flowrate = flowrate(ind);
%DATA
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
figure('Name',' PRESSURE HEAD LOSS vs FLOWRATE')
title('Pressure Head Loss vs FLOW RATE');
xlabel('flow rate');
hold on
[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');
densityplot(flowrate, y, [50,50]);
0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

Risposte (1)

Torsten il 31 Gen 2024
Spostato: Torsten il 31 Gen 2024
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 CommentiMostra NessunoNascondi Nessuno
Giuseppe Zumbo il 31 Gen 2024
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 il 31 Gen 2024
Modificato: Torsten il 31 Gen 2024
"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 ?

Accedi per commentare.

Categorie

Scopri di più su Fit Postprocessing in Help Center e File Exchange

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by