parameter estimation and minimization

6 次查看(过去 30 天)
I have to estimate parameters(k1 ,k2) of the following equation
dx/dt=k1(3-x)-k2*ca0*x^2
I have data available of x vs t at different inital ca0.I have minimization function as follows.
where N=total no. of runs at different ca0(=6) and N'= total no. of runs at same ca0(=12). kindly help me to build a matlab code for the same bcoz i have done only with one summation error minimization.
  4 个评论
Devyani
Devyani 2014-6-28
i am sorry actually i was just trying with four random data points.. i posted tht code only,change the loop to 4 thn.and the function paraesti looks like this
function dxdt=paraesti(x,t,theta)
global theta;
k1=theta(1);
k2=theta(2);
dxdt=k1*(3-x)-k2*ca0*x^2;
Devyani
Devyani 2014-6-28
and i have actually 6 different values of ca0 and each ca0 has 12 data points.thts y the loop 6 and 12

请先登录,再进行评论。

采纳的回答

Star Strider
Star Strider 2014-6-28
Your ODE is not well-behaved. I had problems fitting it with synthesised data with even a small noise component, particularly with respect to k1 (or p(1) in my code). It is very sensitive to the initial parameter estimates. I have no idea what your parameters actually are, so you will have to experiment with the starting estimates to get a good fit to your data. I also used ‘0’ as an initial condition for x in the data synthesis and the DevyaniODEFIT function. Change this if it is not correct. I used a time vector of [0:12] for the independent variable. Change this to reflect the values of your actual independent variable.
I attached the objective function that integrates your ODE and is used by the main script as an objective function that the parameter estimation routine uses to fit your data. I used fminsearch here because I do not know if you have access to the Statistics Toolbox function nlinfit or the Optimization Toolbox function lsqcurvefit. They are more robust and will give a better fit than fminsearch, so you can easily change my code here to use them instead. Note that your paraesti function does not pass Ca0 as a parameter, so either include it in yours, or use my function instead.
The main code:
Ca0=[6 9 12 18];
DE = @(p,t,x,Ca0) p(1).*(3-x)-p(2).*Ca0.*x.^2; % ODE function to create data
% CREATE DATA:
for k1 = 1:length(Ca0)
p(:,k1) = rand(2, 1)*0.1; % Random parameters k1, k2
[t x(:,k1)] = ode45(@(t,x) DE(p(:,k1),t,x,Ca0(k1)), [0:12], 0);
x(:,k1) = x(:,k1) + 0.005*randn(size(x(:,k1))); % Add noise for realism
end
% ESTIMATE PARAMETERS:
for k1 = 1:length(Ca0)
xest = @(B) DevyaniODEFIT(B,t,Ca0(k1)); % Calls ODE function
OLS = @(B) sum( (x(:,1) - xest(B) ).^2); % Least squares objective function
B0 = rand(2,1)*0.1; % Initial parameter estimates
B(:,k1) = fminsearch(OLS, B0); % Optimise parameters to fit data
xgrf(:,k1) = DevyaniODEFIT(B,t,Ca0(k1)); % Generated fitted data to plot
end
figure(1)
plot(t, x, '-x', 'MarkerSize',2, 'LineWidth',1) % Plot original data
hold on % Plot estimated data
plot(t, xgrf, '-p', 'LineWidth',1.5, 'MarkerSize',3);
hold off
grid
The function DevyaniODEFIT.m is attached.
The code definitely works. I documented it with comments as well as I can, so you can easily understand it. You will probably have to experiment with different starting values it to fit your data.
  15 个评论
Devyani
Devyani 2014-7-8
thanks star!! i got to know the prblem and have estimated the parameters with my data :) again thanku !

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Nonlinear Regression 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by