Solving a non-linear least squares inverse problem

I have written the following forward problem. My ultimate goal is to solve the inverse problem for the parameter K. The equation is temperature variation in the half-space due to a time-periodic surface temperature.

I am going to evaluate at a single depth (1.3 m) so z (y) will be a constant. I have real data but am initially just confirming with synthetic data this forward problem produces. My current specific question is in evaluating the least squares solution -- can do this simply in the solver-based optimization (https://www.mathworks.com/help/optim/ug/fmincon.html) but am having trouble conceptualizing how to apply it to my forward problem as written. My analytical forward problem is:

%%ANALYTICAL MODEL
%PARAMETERS
conductivity=.0033; %W m-1 K-1
heat_capacity=671.8; %J kg-1K-1
density=1300; %kgm-1^3
diffusivity=conductivity/(heat_capacity*density);
synodic_period=2.55024e6; %seconds
simulation_period=3*synodic_period;
synodic_frequency=(2*pi)/synodic_period;
T_av=250; %K
T_amp=150; %K
skin_depth=sqrt(2*diffusivity/synodic_frequency); %m
t_list=linspace(0,66162100,1000); %time frame over (S_P)
%t_list = 0:5:synodic_period;
 
z_list=linspace(0,1.5,1000); %depth
T_an=nan(length(t_list),length(z_list)); %output vector
T_an(:, 1) = T_av;
 
%calculate temperature
for t_index=1:length(t_list)
    t=t_list(t_index); 
    for z_index=1:length(z_list)
        z=z_list(z_index);
     
        T_an(t_index,z_index)=T_av+T_amp*exp(-z*sqrt(synodic_frequency/(2*diffusivity)))*cos(synodic_frequency*t-z*sqrt(synodic_frequency/(2*diffusivity)));
        
    end
end

Please let me know if I can provide any additional clarrification. Thank you!

Answer 1

[Sulaymon Eshkabilov]

Hi,

I would suggest to employ curve fit models: (linear least squares or non-linear least squares method)

OPTs = fitoptions( 'Method', 'NonlinearLeastSquares');

MODEL = fittype( MODEL, 'coeff', {'K'});

[L]

Hi, my current update is I have the following for a nonlinear least squares method. (IM_data is my measured data). I am having difficulties with my function T_an -- I am not sure I have set up the functions correctly.

clear all

load IM_data.mat
diffusivity_list = logspace(-8, -3, 1000); %define boundaries
I_list = zeros(size(diffusivity_list));

for n=1:length(I_list)
    I_list(n) = calcI(diffusivity_list(n), IM_data);
end %for n

loglog(diffusivity_list, I_list)
options=optimset('Display','iter');
OPTIONS.TolX = 1e-12;
actual_diffusivity=fminbnd(@(kappa)calcI(kappa, IM_data), 2e-8, 5e-8, options)

function I = calcI(diffusivity, IM_data)

N_temps = size(IM_data,1);
z = 1.3; %meters

T_an = Analytical(IM_data(:,1)', z, diffusivity);
I  = 0.5*sum( (T_an-IM_data(:,2)).^2);

end %function calcI

%%ANALYTICAL MODEL
%PARAMETERS
function T_an = Analytical(t_list, z_list, diffusivity)
%conductivity=.0033; %W m-1 K-1
heat_capacity=671.8; %J kg-1K-1
density=1800; %kgm-1^3
%diffusivity=conductivity/(heat_capacity*density);
synodic_period=2.55024e6; %seconds
simulation_period=3*synodic_period;
synodic_frequency=(2*pi)/synodic_period;
T_av=256; %K
T_amp=150; %K
skin_depth=sqrt(2*diffusivity/synodic_frequency); %m
%t_list=linspace(0,86388,86388);
%t_list=linspace(0,66162100,1000); %time frame over (S_P)
%z_list=linspace(0,1.5,1000); %depth
T_an=nan(length(t_list),length(z_list)); %output vector
T_an(:, 1) = T_av;

%calculate temperature
T_an=T_av+T_amp*exp(-z_list/skin_depth).*cos(synodic_frequency*t_list'-z_list/skin_depth);
end