Optimization problem: Convolution of two functions to dublicate a third one

Hello,

I am looking for a solution for a convolution problem. I have two time dependent functions (F1 and F2). One is a target function (F2). I want to convolute the other function (F1) with a mathematical model with two variables (Axial Dispersion model) so that I end up with a exact dublicate of the F2 function.

F2=F1 conv. AD

Within the axial dispersion model one position (t_mean) and one scattering parameter (Bo) shall be varriied so that the residual sum of squares should be minimized (optimization problem)

This is my code:

tic;

% Parameters of functions Y_Value_Fkt_XX
%BK
c0_tc1 = 0.06582;
k1 = 0.07932;
s1 = 0.09395;
tdead1 = 6.3016;
%KB80
c0_tc2 = 0.06466;
k2 = 0.08148;
s2 = 0.88691;
tdead2 = 6.8116;

% Time and function values
Time_vector = linspace(1, 100, 1200);
Time_matrix = reshape(Time_vector,[],1);
[rows, ~] = size(Time_matrix);
Y_Value_Fkt_BK = zeros(rows, 1);
Y_Value_Fkt_KB80 = zeros(rows, 1);

% Calculate function values
for i = 1:rows
    Y_Value_Fkt_BK(i,1) = 0.5 * c0_tc1 * exp(k1 * (0.5 * k1 * s1 * s1 - (Time_matrix(i,1) - tdead1))) * erfc((k1 * s1 * s1 - (Time_matrix(i,1) - tdead1)) / (sqrt(2) * s1));
    Y_Value_Fkt_KB80(i,1) = 0.5 * c0_tc2 * exp(k2 * (0.5 * k2 * s2 * s2 - (Time_matrix(i,1) - tdead2))) * erfc((k2 * s2 * s2 - (Time_matrix(i,1) - tdead2)) / (sqrt(2) * s2));
end

Y_Value_Fkt_KB80_Vektor=Y_Value_Fkt_KB80(:);
Y_Value_Fkt_BK_Vektor=Y_Value_Fkt_BK(:);

%Initialization of variables of Axial Dispersion Model
Bo_Values_Vector = linspace(0.0001, 0.2, 300);
Bo_Values = reshape(Bo_Values_Vector,[],1);
t_mean_Values_Vector = linspace(0.0001, 5, 300);
t_mean_Values = reshape(t_mean_Values_Vector,[],1);

% Resuts should be saved in
results = table();

% Calculation of Convolution
for j = 1:size(Bo_Values,1)
    for k = 1:size(t_mean_Values,1)
        Bo = Bo_Values(j);
        t_mean = t_mean_Values(k);
        
        % Calculation Y_Value_AD
        Y_Value_AD=zeros(rows,1);
        if Bo ~= 0
            for l=1:size(Time_matrix,1)
                Y_Value_AD(l) = 0.5 * sqrt(Bo / (pi * (Time_matrix(l) / t_mean))) .* exp(-(1 - (Time_matrix(l)/t_mean))^2 * Bo / (4 * (Time_matrix(l) / t_mean)));
            end
        else
            Y_Value_AD(l)=NaN;
        end
        Y_Value_AD_Vektor=Y_Value_AD(:);

        %Convolution
        Conv_AD_BK=conv(Y_Value_AD_Vektor,Y_Value_Fkt_BK_Vektor,'same');

        %normalizing convolution
        Integral_Conv_AD_BK=0;
        for m=2:length(Time_matrix)
            Area=0.5*(Time_matrix(m)-Time_matrix(m-1))*(Conv_AD_BK(m)+Conv_AD_BK(m-1));
            Integral_Conv_AD_BK=Integral_Conv_AD_BK+Area;
        end
        Norm_Conv_AD_BK=Conv_AD_BK./Integral_Conv_AD_BK;
                    
        
        % RSS Calculation (Residual sum of squares)
        RRsum = sum((Norm_Conv_AD_BK - Y_Value_Fkt_KB80).^2);
        
        % Save 
        results = [results; table(Bo, t_mean, RRsum)];
    end
end


%Find row where RSS=min
[minRRsum, minIndex] = min(results.RRsum);
minRow = results(minIndex, :);

%Find Bo and t_mean for RSS=min
minBo = minRow.Bo;
mint_mean = minRow.t_mean;

%Calculation of Y_Values and convolution for RSS=min values
Conv_min=zeros(rows,1);
Y_Value_AD_min=zeros(rows,1);
for i = 1:size(Time_matrix,1)
    if mint_mean ~= 0
        Y_Value_AD_min(i,1)= 0.5 * sqrt(minBo / (pi * (Time_matrix(i, 1) / mint_mean))) * exp(-(1 - Time_matrix(i, 1))^2 * minBo / (4 * (Time_matrix(i, 1) / mint_mean)));
    else
        Y_Value_AD_min(i, 1) = NaN;
    end
end
    Conv_min=conv(Y_Value_AD_min,Y_Value_Fkt_BK,'same');
    Integral_Conv_min=0;
    for n=2:length(Time_matrix)
        Area_min=0.5*(Time_matrix(n)-Time_matrix(n-1))*(Conv_min(n)+Conv_min(n-1));
    Integral_Conv_min=Integral_Conv_min+Area_min;
    end
    Norm_Conv_min=Conv_min./Integral_Conv_min;

% control Plot
plot(Time_matrix, Y_Value_Fkt_BK, 'b', Time_matrix, Y_Value_Fkt_KB80, 'r', Time_matrix, Norm_Conv_min, 'g');

time=toc;
disp(time);

The following picture shows the result. Red is the target function, blue is the function which needs to be convoluted. green should be a duplicate of red if everything goes well.

Do you have any suggestions where a fault may be hidden or how to get a better fit?

Thanks in advance for your help :)

[Harald]

Hi,

it might be that within the range of values you have tried, this just is the best solution.

If you have Optimization Toolbox, you could let MATLAB take care of the optimization portion. To get started, consider the free Optimization Onramp. There is also a one-day training covering this more in depth:

https://de.mathworks.com/learn/training/optimization-techniques-in-matlab.html

Using the conv function of MATLAB and avoiding nested for-loops may speed up your code.

Best wishes,

Harald