Caclulate Horizonzal Area Between Two Curves And Minimize This Area (Optimization Problem)
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Jan Sander
el 25 de Sept. de 2020
Comentada: Jan Sander
el 29 de Sept. de 2020
Hello everybody,
Orginally i have 2 curves, lets say curve A und curve B, that calculates curve C (A-B=C). Curve A and B can be shifted in x-direction. This curve C needs to be approximated to antoher curve D (reference curve).
i want Matlab to find the optimimum overlapping between Curve C and Curve D. Thus, Curve C needs to be shifted in x-direction. I think a simple quality factor is the area between Curve C and D. Since i am shifting the Curve C in x-direction, the horizonzal area between curve C and D is interesting and shall be minimized. Whith this optimization i want an output (x-value) on how curve C is shifted in x-direction.
1) How do i calculate the horzontal area between two curves?
(but i dont get it... the graph is just flipped, the value of my calculated area didnt change)
2) How do i get this area minimized?
For example i have simulated data in the following:
%Example Data, y1 = Curve D (reference), y2 = Curve C (C=A-B)
xscale = -11; %Scaling/Shifting parameter = optimization variable
x1=2:0.1:40;
x2=2:0.1:20;
y1_roh = -0.5*x1.^3+50000;
y2_roh = -0.08*(x2-xscale).^3+50000;
del1 = find(y1_roh>4.9973e+04);
y1_roh(del1) =[];
x1(del1)=[];
del2 = find(y1_roh<48750)
y1_roh(del2) =[];
x1(del2)=[];
del3 = find(y2_roh>4.9973e+04)
y2_roh(del3) =[];
x2(del3)=[];
del4 = find(y2_roh<48750)
y2_roh(del4) =[];
x2(del4)=[]
%Plot
figure(1)
plot(x1,y1_roh,'k')
hold on
plot(x2,y2_roh,'b')
hold off
%Example Fit (needed to be done in my work & for function building)
xfit1=x1(1):0.01:x1(end);
xfit2=x2(1):0.01:x2(end);
p1 = polyfit(x1,y1_roh,40);
y1 = polyval(p1,xfit1);
p2 = polyfit(x2,y2_roh,20);
y2 = polyval(p2,xfit2);
figure(2)
plot(xfit1,y1,'k')
hold on
plot(xfit2,y2,'b')
hold off
%horizontal Area caclulation?
%Optimization? Minimum Area?
(I already tried to do the primivitve integral of my polynomial fit via polyint and then calculate the area and tried to build in my optimazition function, but unfortunatly i faild. And this still gives me the area between the curve and x-axis)
Thanks in adavacne!
4 comentarios
J. Alex Lee
el 25 de Sept. de 2020
take the square, or absolute value of the differences b/w curves, then integrate that?
but if instead of integrate, you just add up the points, i think you are left with a typical least squares type statement, which you can usually make well defined solvers for.
it's also unclear if the curves you care about are inherently discrete in description, or if you always have analytical parametric exprssions for them. there may be better strategies for comparing analytical curves, rather than discretizing them on arbitrary domains and comparing the discrete values.
J. Alex Lee
el 25 de Sept. de 2020
another point to add is that your answer will depend on over what values of y you care about the curves. from your example code, it isn't clear that you're looking at the same range, though it looks roughly between 4.88e4 and 5e4.
Respuesta aceptada
J. Alex Lee
el 25 de Sept. de 2020
Take a look at this:
clc;
close all;
clear;
% it looks like these are ultimately the curves you want to compare
f1 = @(x)(-0.5*x.^3+50000);
f2 = @(x,s)(-0.08*(x-s).^3+50000);
% i don't think you need to look for cross overs, just decide on what
% values of y you want to match over by shifting in x
% y = a (x-s)^3 + b -> x = nthroot((y-b)/a,3)+s
g1 = @(y)(nthroot((y-50000)/-0.5,3));
g2 = @(y,s)(nthroot((y-50000)/-0.08,3)+s);
% decide, say (according to your original plots)
ylim = [4.88e4,4.98e4];
% which very roughly comes out to (according to your original plots)
xlim = [2,14];
ylist = linspace(ylim(1),ylim(2),100)';
% your first example / initial shift is
s0 = -11;
% here is how things look with the shift of s0=-11
figure(1); cla; hold on;
ax0 = gca;
% plot over x limits
fplot(f1,xlim,'-','Color',[1 1 1]*0.7)
fplot(@(x)f2(x,s0),xlim,'Color',[1 1 1]*0.7)
% plot over y limits
ph(1) = plot(g1(ylist),ylist,'--','LineWidth',1.5)
ph(2) = plot(g2(ylist,s0),ylist,'--','LineWidth',1.5)
yline(ylim(1))
yline(ylim(2))
% now define a bunch of different shift factors
N = 50;
sList = linspace(-18,7,N)';
% evaluate the different cost functions as function of shift s
for i = N:-1:1
res_hd(i) = horiz_cost_discrete(sList(i),ylist,g1,g2);
res_hi(i) = horiz_cost_int(sList(i),ylim,g1,g2);
res_wr(i) = horiz_cost_int_wrong(sList(i),ylim,g1,g2);
end
% determine optimum s for minimizing cost for different cost functions
% (except the one that won't work without the absolute value)
opt_hd = fminsearch(@(s)horiz_cost_discrete(s,ylist,g1,g2),s0)
opt_hi = fminsearch(@(s)horiz_cost_int(s,ylim,g1,g2),s0)
% these would change depending on what ylim you set
figure(2);
ax(1) = subplot(2,1,1); cla; hold on;
plot(sList,res_hd,'.')
xline(opt_hd)
title('Cost functions')
ylabel('based on discrete summing')
ax(2) = subplot(2,1,2); cla; hold on;
plot(sList,res_hi,'.')
plot(sList,res_wr,'.')
linkaxes(ax,'x');
xline(opt_hi)
ylabel('based on integration')
ph(3) = plot(ax0,g2(ylist,opt_hd),ylist,'-','LineWidth',1.5)
legend(ax0,ph,'String',{'Curve 1 (static)','Curve 2 (shifted by s0)','Curve 2 optimal'});
%% cost function based on discretized functions
function res = horiz_cost_discrete(s,ylist,g1,g2)
x1list = g1(ylist);
x2list = g2(ylist,s);
res = sum(abs(x1list-x2list));
end
%% cost function based on integration of absolute value
function res = horiz_cost_int(s,ylim,g1,g2)
df = @(y)( abs(g1(y)-g2(y,s)) );
res = integral(df,ylim(1),ylim(2));
end
%% cost function based on integration of just the difference between curves
function res = horiz_cost_int_wrong(s,ylim,g1,g2)
df = @(y)( (g1(y)-g2(y,s)) );
res = integral(df,ylim(1),ylim(2));
end
not sure how you'd implement Matt J's suggestion.
8 comentarios
J. Alex Lee
el 28 de Sept. de 2020
in retrospect, i guess this is very closely related to matching the centers of mass of the "cloud" of points defining the curves. In this sense, you can get pretty close to above by just doing
cntrC = sum([CD,CB])/numel(CD)
cntrD = sum([DD,DB])/numel(DD)
dCntrs = cntrC-cntrD
% return
% plot result of simple least squares
figure(1); cla; hold on;
title('simple least squares')
plot(CD, CB, '-', 'LineWidth', 4,'Color',[1 1 1]*0.8) % Curve C original
plot(DD, DB, '-', 'LineWidth', 4,'Color',[1 1 1]*0.8) % Curve D original
plot(CD-dCntrs(1), CB, '-b', 'LineWidth', 1) %Curve C shifted in x
plot(CD-dCntrs(1), CB-dCntrs(2), '--r', 'LineWidth', 1) %Curve C shifted in both x and y
Más respuestas (1)
Matt J
el 25 de Sept. de 2020
Editada: Matt J
el 26 de Sept. de 2020
EDIT: If you can accept the minimum squared area difference as a criterion, it has a known analytical solution, which is to take the shift x which maximizes the correlation of C and D.
[~,imax]=max(conv(flip(C),D,'same'));
x_optimal=x(imax);
5 comentarios
Matt J
el 26 de Sept. de 2020
I cannot say what the "best" error criterion is, since that depends on Jan's needs, but I can say that the L2 error criterion is the most standard, and if it is sufficient for Jan, it has the advantage of an analytical solution.
J. Alex Lee
el 26 de Sept. de 2020
Editada: J. Alex Lee
el 26 de Sept. de 2020
i concede I should have used L2 and taken the square, but this distinction is less interesting than the one that i mean, which is about taking differences in x direction vs y direction.
i guess you can make it equivalent to the horizontal area metric by interpolating on the y domain and convolving on the x coordinates. either way, despite the "analytical" solution, there's some nontrivial preprocessing necessary on the "signals"
of course i completely agree what's "better" depends on the need.
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!