matlab动力学参数,用matlab求动力学参数代码 求指教

程序我仿照写出来了　　用的是2010年数据　　用fmincon优化(文章中用fminsearch,但是我不知道用这个该怎么加lb和ub)　但是估计的ｋ跟文章中２０１０年相差很大

k0 = [0.5  0.5  0.5  0.5  0.5 0.5];         % 参数初值我是在lb和ub范围内随便选的,是不是会对结果有很大影响？

文章中E (0) in [1.8 × 104,1.4 × 107]  I(0) in [1.8 × 103,1.8 × 105], 但是我不知道怎么将ｙ的初值设成一个范围，于是用2.043*10^4和 1.8001*10^3代替了，有什么办法可以改一下吗

%6个参数 运行不报错  但是参数值和文章相差很大 结果拟合也不好

%function k1k2k3

format long

clear all

clc

tspan = [1:12];

k0 = [0.5  0.5  0.5  0.5  0.5 0.5];         % 参数初值

lb = [0.0001 0.0001 0.0001 0.1 0.01 0.1];                   % 参数下限

ub = [50 1 1 1 1 1];    % 参数上限

x0 = [1.4*10^8 2.043*10^4 1.8001*10^3 1212 0];    %y初值

ExpData=[

1 3756.7

2 2386.2

3 7775.6

4 24860.9

5 35434.7

6 34310

7 26126.3

8 11909.6

9 10165.4

10 8761.2

11 7959.1

12 6087.9

];

t=ExpData(:,1);

yexp = ExpData(:,2);  % yexp: CDC数据

% 使用函数fmincon()进行参数估计

[k,fval,flag] = fmincon(@MSS,k0,[],[],[],[],lb,ub,[],[],x0,yexp);

fprintf('\n使用函数fmincon()估计得到的参数值为:\n')

fprintf('\tk1 = %.4f\n',k(1))

fprintf('\tk2 = %.4f\n',k(2))

fprintf('\tk3 = %.4f\n',k(3))

fprintf('\tk4 = %.4f\n',k(4))

fprintf('\tk5 = %.4f\n',k(5))

fprintf('\tk6 = %.4f\n',k(6))

fprintf('  The sum of the squares is: %.1e\n\n',fval)

%% =======绘图显示结果=======

[tt xx] = ode45(@li,tspan,x0,[],k);

plot(t,yexp(:,1),'r.-')

hold on

plot(tt,xx(:,4),'ro-')

legend('exp y1', 'pre y1')

xlabel('t/月')

ylabel('人数')

function f = MSS(k,x0,yexp)           % fmincon()

tspan = [1:12];

[t x] = ode45(@li,tspan,x0,[],k);

y(:,1) = x(:,4);

f = sum((log2(1+ yexp(:,1))-log2(1+ y(:,1))).^2);

%微分方程组

function dydt=li(t,y,k)

u=3.9139*10^-5;

dydt=[(u*(y(1)+y(2)+y(3)+y(4)+y(5))-k(1)*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))+k(2)*y(5)-u*y(1))

(k(1)*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))-k(3)*y(2)-u*y(2))

(k(3)*(1-k(5))*y(2)-k(4)*y(3)-u*y(3))

(k(3)*k(5)*y(2)-k(6)*y(4)-u*y(4))

(k(4)*y(3)+k(6)*y(4)-k(2)*y(5)-u*y(5))];

使用函数fmincon()估计得到的参数值为:

k1 = 0.0001

k2 = 0.3354

k3 = 0.4065

k4 = 0.8597

k5 = 1.0000

k6 = 0.1000

The sum of the squares is: 1.2e+01