Fitting Numerical Solution of PDE to Given Data with Multiple Parameters

4 views (last 30 days)
I'm trying to fit my solution to be as close as possible to given data for my PDE system. I've been using the following as a reference (https://www.mathworks.com/matlabcentral/answers/406938-fitting-the-numerical-solution-of-a-pde-with-one-parameter-to-some-data) which has been super helpful. I was able to reproduce the results from the example (shown below), and I've adapted the code to my PDE system (which I'm solving for 10 parameters instead of one).
But as I run my code, the parameters are not adjusting at all, and it produce the following result.
As the message said, I adjusted the step tolerance to be lower (10^(-10)) so that it may continue find the optimal parameters. But as its being solved, the parameters don't adjust from the initialized ones. Is there a possible explanation and a remedy for this issue?
(Note: I'll provide my code as a comment to this question post. But of course, I won't be able to provide the data I'm using to compare to train the optimization algorithim )
global exper
exper=[-0.4,-112,-121,-128,-131,-131.1,-133.2,-133.1,-134.4,-134.7,-135.12,-135.7,-135.8,-135.1,-135.4,-135.8,-135.9,-134,-135.7,-135.1,-135.3,-135.4,-136,-139,-136];
% p=[0.218,0.52,0.0115,2.03,31.6];
% qr f a n ks
p0=[0.218,0.52,0.0115,2.03,31.6];
lb=[0.01,0.4,0,1,15]; %lower bound
ub=[Inf,Inf,10,10,100]; %upper bound
options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter');
[p,resnorm,residual,exitflag,output] = lsqnonlin(@richards1,p0,lb,ub,options)
uu=richards1(p);
global t
plot(t,uu+exper')
hold on
plot(t,exper,'ko')
function uu=richards1(p)
% Solution of the Richards equation
% using MATLAB pdepe
%
% $Ekkehard Holzbecher $Date: 2006/07/13 $
%
% Soil data for Guelph Loam (Hornberger and Wiberg, 2005)
% m-file based partially on a first version by Janek Greskowiak
%--------------------------------------------------------------------------
L = 200; % length [L]
s1 = 0.5; % infiltration velocity [L/T]
s2 = 0; % bottom suction head [L]
T = 4; % maximum time [T]
% qr = 0.218; % residual water content
% f = 0.52; % porosity
% a = 0.0115; % van Genuchten parameter [1/L]
% n = 2.03; % van Genuchten parameter
% ks = 31.6; % saturated conductivity [L/T]
x = linspace(0,L,100);
global t
t = linspace(0,T,25);
options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7)
u = pdepe(0,@unsatpde,@unsatic,@unsatbc,x,t,options,s1,s2,p);
global exper
uu=u(:,1,1)-exper';
% figure;
% title('Richards Equation Numerical Solution, computed with 100 mesh points');
% subplot (1,3,1);
% plot (x,u(1:length(t),:));
% xlabel('Depth [L]');
% ylabel('Pressure Head [L]');
% subplot (1,3,2);
% plot (x,u(1:length(t),:)-(x'*ones(1,length(t)))');
% xlabel('Depth [L]');
% ylabel('Hydraulic Head [L]');
% for j=1:length(t)
% for i=1:length(x)
% [q(j,i),k(j,i),c(j,i)]=sedprop(u(j,i),p);
% end
% end
%
% subplot (1,3,3);
% plot (x,q(1:length(t),:)*100)
% xlabel('Depth [L]');
% ylabel('Water Content [%]');
% -------------------------------------------------------------------------
function [c,f,s] = unsatpde(x,t,u,DuDx,s1,s2,p)
[q,k,c] = sedprop(u,p);
f = k.*DuDx-k;
s = 0;
end
% -------------------------------------------------------------------------
function u0 = unsatic(x,s1,s2,p)
u0 = -200+x;
if x < 10 u0 = -0.5; end
end
% -------------------------------------------------------------------------
function [pl,ql,pr,qr] = unsatbc(xl,ul,xr,ur,t,s1,s2,p)
pl = s1;
ql = 1;
pr = ur(1)-s2;
qr = 0;
end
%------------------- soil hydraulic properties ----------------------------
function [q,k,c] = sedprop(u,p)
qr = p(1); % residual water content
f = p(2); % porosity
a = p(3); % van Genuchten parameter [1/L]
n = p(4); % van Genuchten parameter
ks = p(5); % saturated conductivity [L/T]
m = 1-1/n;
if u >= 0
c=1e-20;
k=ks;
q=f;
else
q=qr+(f-qr)*(1+(-a*u)^n)^-m;
c=((f-qr)*n*m*a*(-a*u)^(n-1))/((1+(-a*u)^n)^(m+1))+1.e-20;
k=ks*((q-qr)/(f-qr))^0.5*(1-(1-((q-qr)/(f-qr))^(1/m))^m)^2;
end
end
end
  1 Comment
Thomas Rodriguez
Thomas Rodriguez on 31 Dec 2022
Initial Parameters are given as which are the same values called upon from (PDE_Coefficients(9)).
Code:
clc; clear; close all;
% Defining Actual Data:clc; clear; close all;
% Defining Actual Data:
global actual
Actual_1D = load('Actual_1D.mat', 'Actual_1D');
Actual_1D = Actual_1D.Actual_1D;
actual = Actual_1D;
% Defining Initial Parameters:
[D1,D2,B_e,B_s,B_a,B_v,theta,lambda,gamma,Omega,rho_e,rho_s,rho_a] = PDE_Coefficients(9);
% Storing Initial Parameters that are being varied:
p0 = [theta, lambda, gamma, B_e, B_v, D1, D2, rho_e, rho_s, rho_a];
% lower bound
lb=[10^(-4), 10^(-4), 10^(-4), 10^(-15), 10^(-20), 10^(-8), 10^(-8), 1000, 10^(6), 10^(4)];
% upper bound
ub=[1.0, 1.0, 1.0, 10^(-4), 10^(-4), 10^(-1), 10^(-1), 10^(4), 9*10^(8), 10^(6)];
options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter'); % ,'StepTolerance',10^(-9),'MaxIterations',30 ,'FunctionTolerance',10^(-8)
[p,resnorm,residual,exitflag,output] = lsqnonlin(@PDE_Model,p0,lb,ub,options)
% Plotting:
global t x
global loc
loc = 12; % Specific Location being examined
UU = PDE_Model(p);
figure; hold on; grid on; box on;
plot(t,UU + actual(:,loc))
hold on;
plot(t,actual(:,loc),'ro')
xlabel('t')
ylabel('U')
title(['Population Density at ' num2str(x(loc)) '$$^\circ$$'], 'Interpreter','latex')
legend('Model','Actual Data','Location','best')
figure; hold on;
global Approx
imagesc(x,t,Approx); % Contour Plot
grid on; box on;
xlabel('Latitude$^\circ$',"Interpreter","latex","FontSize",12)
ylabel('Time (Weeks)',"Interpreter","latex","FontSize",12)
title('Approximated: $$\theta \lambda \hat{E}$$',"Interpreter","latex","FontSize",12)
axis tight
colormap jet
colorbar
hold off;
figure; hold on;
imagesc(x,t,actual); % Contour Plot
grid on; box on;
xlabel('Latitude$^\circ$',"Interpreter","latex","FontSize",12)
ylabel('Time (Weeks)',"Interpreter","latex","FontSize",12)
title('Actual: $$\theta \lambda E$$',"Interpreter","latex",'FontSize',12)
axis tight
colormap jet
colorbar
hold off;
% Printing Table:
% p = [theta, lambda, gamma, B_e, B_v, D1, D2, rho_e, rho_s, rho_a];
fprintf('\n')
fprintf('Optimal Parameter Values \n')
fprintf('-----------------------\n')
fprintf('D1 %3.4e \n', p(6))
fprintf('D2 %3.4e \n', p(7))
fprintf('theta %3.4e \n', p(1))
fprintf('lambda %3.4e \n', p(2))
fprintf('gamma %3.4e \n', p(3))
fprintf('B_e %3.4e \n', p(4))
fprintf('B_v %3.4e \n', p(5))
fprintf('rho_e %3.4e \n', p(8))
fprintf('rho_s %3.4e \n', p(9))
fprintf('rho_a %3.4e \n', p(10))
fprintf('-----------------------\n')
function UU = PDE_Model(p)
% Defining Space and Time:
global A B
[A,B] = Loading_Data(); % Loading Dataframes
global x t t0 tf
% Space (x):
x = sort(A(:,1));
t0 = 9; tf = 37;
% Time (t):
t = [t0:1:tf]';
% Solving Equation:
fprintf('Solving PDE System: \n')
options = odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t,options,p);
% Solutions for the 6 Populations:
E = sol(:,:,2);
global actual loc
fprintf('Integrating Newly Infected Symptomatic Data, theta*lambda*E: \n')
N = length(t); theta = p(1); lambda = p(2);
global Approx
Approx = RK4(E,t0,tf,N,theta,lambda);
fprintf('Computing Error Matrix, UU: \n\n')
UU = Approx - actual;
% disp(UU)
UU = UU(:,loc);
% pdepe form:
% c(x,t,u,u_x)u_t = x^(-m)*(x^(m)*f(x,t,u,u_x))_x + s(x,t,u,u_x)
function [c,f,s] = pdefun(x,t,u,dudx,p)
% Parameters:
theta = p(1);
lambda = p(2);
gamma = p(3);
B_e = p(4);
B_v = p(5);
D1 = p(6);
D2 = p(7);
rho_e = p(8);
rho_s = p(9);
rho_a = p(10);
[B_s, B_a] = transmission_constant(t);
Omega = Omega_Surface_Trans(7*t);
% pdepe form:
% c(x,t,u,u_x)*u_t = x^(-m)*(x^(m)*f(x,t,u,u_x))_x + s(x,t,u,u_x)
% c(x,t,u,u_x):
% 6 populations for 6 system of eq'ns
c = ones([6 1]);
% f(x,t,u,u_x,u_y):
f = [D1; D1; D1; D1; D1; D2].*dudx;
% s(x,t,u,u_x):
F1 = -(B_e*u(2) + B_s*u(3) + B_a*u(4) + B_v*u(6)).*u(1);
F2 = (B_e*u(2) + B_s*u(3) + B_a*u(4) + B_v*u(6)).*u(1) - lambda*u(2);
F3 = theta*lambda*u(2) - gamma*u(3);
F4 = (1 - theta)*lambda*u(2) - gamma*u(4);
F5 = gamma*(u(3) + u(4));
F6 = rho_e*u(2) + rho_s*u(3) + rho_a*u(4) - Omega*u(6);
F = [F1;F2;F3;F4;F5;F6];
s = F;
end
% Initial Conditions:
function [u0] = pdeic(x,p)
% Input: scalar of space array (x) i.e. 1x1
% Initial Population:
[S0,E0,I_S0,I_A0,R0,V0] = Population(x);
% Storing Initial Pop. into u0:
u0 = [S0; E0; I_S0; I_A0; R0; V0];
end
% Establishing Boundary Conditions:
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t,p)
% Standard Boundary Conditions Solver:
% p(x,t,u) + q(x,t)*f(x,t,u,u_x) = 0
% No Flux Boundary Conditions:
% Domain: (ymin = 32.55, ymax = 33.03)
% Boundary Conditions on the Left:
% S_x(32.55,t) = E_x(32.55,t) = I_s_x(32.55,t) =
% I_a_x(32.55,t) = R_x(32.55,t) = V_x(32.55,t) = 0, etc.
pl = zeros([6 1]);
ql = ones([6 1]);
% Boundary Conditions on the Right:
% S_x(33.03,t) = E_x(33.03,t) = I_s_x(33.03,t) =
% I_a_x(33.03,t) = R_x(33.03,t) = V_x(33.03,t) = 0, etc.
pr = pl;
qr = ql;
end
% Determing B_s and B_a:
function [B_s, B_a] = transmission_constant(t)
% B_s and B_a are defined as such in
% Angleica Bloomquist, Naveen Vaidya, et. al. ”Modeling Risk of Environmental Reservior of SARS-CoV-2”
% Input: t - time (scalar)
% t0 = 9; % initial time
% tf = 38; % final time
tol = 10^(-2); % tolerance condition
n = 1000; % number of points in stage arrays
% First Stage: 03/07/2020 - May 10, 2020 (~ 64 days ~ 9 weeks)
t1 = linspace(t0,18,n);
% Second Stage: May 11, 2020 - July 14,2020 (~ 64 days ~ 9 weeks)
t2 = linspace(18.01,28,n);
% Third Stage: July 15, 2020 - September 16, 2020 (~ 63 days ~ 9 weeks)
t3 = linspace(28.01,tf,n);
for i = 1:n
if (t1(i) == t) | (abs(t - t1(i)) <= tol)
% First Stage: 03/07/2020 - May 10, 2020 (~ 64 days ~ 9 weeks)
B_s = 8.92*10^(-8);
B_a = 2.67*10^(-8);
break
elseif (t2(i) == t) | (abs(t - t2(i)) <= tol)
% Second Stage: May 11, 2020 - July 14,2020 (~ 64 days ~ 9 weeks)
B_s = 3.32*10^(-8);
B_a = 2.49*10^(-8);
break
elseif (t3(i) == t) | (abs(t - t3(i)) <= tol)
% Third Stage: July 15, 2020 - September 16, 2020 (~ 63 days ~ 9 weeks)
B_s = 1.64*10^(-8);
B_a = 1.54*10^(-8);
break
else
continue
end
end
end
% Omega function:
function [omega] = Omega_Surface_Trans(t)
Final_tmp = Final_Temperature(t);
% Defining a and b values:
a = 0.094519;
b = -6.4829;
omega = a*Final_tmp + b;
% Note: a and b values are the coefficient averages over all surfaces
end
% Runge-Kutta (4th Order):
function Approx = RK4(E,t0,tf,N,theta,lambda)
h = (tf - t0)/N;
tn = [t0:h:tf]';
% Calculating Cumulative Data:
% f(t) = \theta*\lambda*E(t,x)
% E(t,x) assumed to be 29x30 size
% theta = 0.6;
% lambda = 0.21;
f = @(t,E) theta*lambda*E;
% Initial Condition: y(1) = y0
Y(1,:) = theta*lambda*E(1,:);
% Runge-Kutta Method:
for i = 1:N-1
% k1 = f( t_n, y_n )
k1 = f(tn(i), E(i,:));
% k2 = f( t_n + 0.5*h, y_n + h*0.5*k1 )
k2 = f(tn(i) + 0.5*h, E(i,:) + h*0.5*k1);
% k3 = f( t_n + 0.5*h, y_n + h*0.5*k2 )
k3 = f(tn(i) + 0.5*h, E(i,:) + h*0.5*k2);
% k4 = f( t_n + h, y_n + h*k3 )
k4 = f(tn(i) + 0.5*h, E(i,:) + h*k3);
% f(t_n+1, y_n+1)
Y(i+1,:) = Y(i,:) + (1/6)*h*(k1 + 2*k2 + 2*k3 + k4);
end
Approx = Y;
end
function [S0,E0,I_S0,I_A0,R0,V0] = Population(x)
% fprintf('x input: %2.4f \n', x)
[S,E,I_S,I_A,R,V] = Initial_Pop();
tol = 10^(-8); % Defined Tolerance:
% [A,~] = Loading_Data();
Lat = sort(A(:,1)); % Latitude Coordinates:
for i = 1:length(Lat)
if (abs(Lat(i) - x) <= tol)
S0 = S(i);
E0 = E(i);
I_S0 = I_S(i);
I_A0 = I_A(i);
R0 = R(i);
V0 = V(i);
break
else
continue
end
end
end
end
% Calculating Initial Population Function:
function [S0,E0,I_S0,I_A0,R0,V0] = Initial_Pop()
global A B
Lat = A(:,1); Zip = string(A(:,4));
Pop = A(:,3); Case_Data = B(2,:); % Week 1: 03/1/2020 - 03/07/2020
% % Calculating Initial Population Distribution:
% Initial Susceptible Population Density:
[S0] = Initial_Population_Density(Lat,Pop,Zip,"Initial Susceptible $$S(0,X)$$ Population Density",0);
n = length(S0);
% Initial Exposed Population Density:
E0 = zeros(1,n);
% Initial Symptomatic Infected Population Density:
I_S0 = Initial_Population_Density(Lat,Case_Data,Zip,"Initial Susceptible $$S(0,X)$$ Population Density",0);
% Initial Asymptomatic Infected Population Density:
I_A0 = zeros(1,n);
% Initial Recovered Population Density:
R0 = zeros(1,n);
% Initial Environment Population Density:
V0 = zeros(1,n);
end
Final Temperature Function:
% Temperature Functions:
function [F] = Final_Temperature(x)
% Note: x = [0:(1/24):365]' gives the Net Temperature fluctation for One Year
f_s = Seasonal_Temp(x);
f_d = Dirunal_Temp(x);
F = f_s + f_d;
end
% Seasonal Temperature Fluctuation Function:
function [f_s] = Seasonal_Temp(x)
t_0 = 64.712;
epsilon_m = -6.9178;
tau_m = 365;
phi_m = 0.99605;
f_s = t_0 + epsilon_m*sin((2*pi/tau_m)*x + phi_m);
end
% Dirunal Temperature Fluctuation Function:
function [f_d] = Dirunal_Temp(x)
epsilon_d = -4.3372;
tau_d = 1;
phi_d = -5.0872;
m = 0;
f_d = m + epsilon_d*sin((2*pi/tau_d)*x + phi_d);
end

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!