Finding parameters by fitting data to a system of ODEs with lsqnonlin

Dear Matlab Team,

I have a system of ODEs, with three parameters (s,d,gamma). I want to find the parameters by fitting data to this system:

First eq1: dT0/dt = s - d*T0 - gamma * T0

Second eq2: dT1/dt = 2*d*T0 + d*T1 - gamma * T1

I know that T0 converge to 2016 and T1 converge to 42 in steady state and I can find T0 and T1 values over time. Now, I want to find estimation for (s,d,gamma) by lsqnonlin.

Thanks for your help and time !!

Answer 1

[prabhat kumar sharma]

Heelo neda,

To estimate the parameters ( s ), ( d ), and ( \gamma ) for your system of ODEs using lsqnonlin in MATLAB, you can follow these steps. This approach involves defining your system of ODEs, simulating it over time, and using lsqnonlin to minimize the difference between the simulated and observed data.

You can use the below steps as reference.

1. Define the System of ODEs: Create a function that computes the derivatives based on your current estimates of the parameters.
2. Simulate the ODEs: Use ode45 or another ODE solver to simulate the system over time.
3. Objective Function: Define an objective function that computes the difference between the simulated results and your observed data.
4. Use lsqnonlin: Call lsqnonlin to minimize the objective function and estimate the parameters.

function estimate_parameters
    % Observed steady-state values
    T0_ss = 2016;
    T1_ss = 42;
    
    % Initial guess for parameters [s, d, gamma]
    initial_guess = [1000, 0.1, 0.1];
    
    % Time span for simulation
    tspan = [0, 100];
    
    % Initial conditions for T0 and T1
    T0_initial = 0;
    T1_initial = 0;
    initial_conditions = [T0_initial, T1_initial];
    
    % Define the objective function
    objective = @(params) model_error(params, tspan, initial_conditions, T0_ss, T1_ss);
    
    % Set options for lsqnonlin
    options = optimoptions('lsqnonlin', 'Display', 'iter', 'TolFun', 1e-6, 'TolX', 1e-6);
    
    % Estimate parameters using lsqnonlin
    [estimated_params, resnorm] = lsqnonlin(objective, initial_guess, [], [], options);
    
    % Display the estimated parameters
    fprintf('Estimated Parameters:\n');
    fprintf('s = %.4f\n', estimated_params(1));
    fprintf('d = %.4f\n', estimated_params(2));
    fprintf('gamma = %.4f\n', estimated_params(3));
end

function error = model_error(params, tspan, initial_conditions, T0_ss, T1_ss)
    % Unpack parameters
    s = params(1);
    d = params(2);
    gamma = params(3);
    
    % Solve the ODE system
    [~, T] = ode45(@(t, y) odesystem(t, y, s, d, gamma), tspan, initial_conditions);
    
    % Compute the error between the simulated steady-state values and the observed values
    T0_error = T(end, 1) - T0_ss;
    T1_error = T(end, 2) - T1_ss;
    
    % Return the error as a vector
    error = [T0_error; T1_error];
end

function dydt = odesystem(t, y, s, d, gamma)
    % Unpack the variables
    T0 = y(1);
    T1 = y(2);
    
    % Define the ODEs
    dT0dt = s - d * T0 - gamma * T0;
    dT1dt = 2 * d * T0 + d * T1 - gamma * T1;
    
    % Return the derivatives
    dydt = [dT0dt; dT1dt];
end

I hope it helps!

[Neda]

Thank you so much. It was so helpful. I have two questions:

Actually, I have values of T0 and T1 over time, for tspan=[0 500], with time step = 5, I have 100 T0 and 100 T1. How can I write error?

Second: I want to plot the observed and fitted solution?

Thanks

Answer 2

[Torsten]

syms T1(t) T0(t) s d g
eqn1 = diff(T0,t) == s - (d+g) *T0 ;
eqn2 = diff(T1,t) == 2*d*T0 + (d-g)*T1;
sol = dsolve([eqn1,eqn2])

sol = struct with fields:
    T1: exp(t*(- g + d))*(- (s*exp(- d*t + g*t))/(d - g) + C2) + exp(-t*(g + d))*(- (s*exp(d*t)*exp(g*t))/(d + g) + C1)
    T0: -exp(-t*(g + d))*(- (s*exp(d*t)*exp(g*t))/(d + g) + C1)

simplify(sol.T0)

ans = 

simplify(sol.T1)

ans = 

●

Now you have the symbolic solutions and you know that s/(g+d) equals 2016 and s/(g-d) equals 42+2016. So you are left with only one parameter to fit.

[Star Strider]

@Neda — Use the matlabFunction function to auttomatically create an anonymous function from your symbolic code.

You can then use that anonymous function with any optimisatiion function.