ode45和lsqcurvefit结合，对微分方程组的参数求最优解

IBB-EUT · 发表于 2021-2-24 13:48

EDA365欢迎您登录！

您需要登录才可以下载或查看，没有帐号？注册

x

采用ode45和lsqcurvefit结合的方法，对微分方程组的参数求最优解，效果如图，第二张表不理想，
感觉只获得了局部最优解。如何优化，不知有没有大神能否提供有效解决方法。

微分方程：
N = 11000000/250;
dydt(1)=-k(1)*(y(3)/N)*y(1)-k(2)*(y(4)/N)*y(1);
dydt(2)=k(1)*(y(3)/N)*y(1)+k(2)*(y(4)/N)*y(1)-k(3)*y(2);
dydt(3)=k(3)*k(4)*y(2)-k(7)*y(3);
dydt(4)=k(3)*k(5)*y(2)-k(6)*y(4);
dydt(5)=k(3)*(1-k(4)-k(5))*y(2);
dydt(6)=k(6)*y(4)+k(7)*y(3)-k(8)*y(6)-k(9)*y(6);
dydt(7)=k(8)*y(6);
dydt(8)=k(9)*y(6);
dydt(9)=dydt(3)+dydt(4)+dydt(6);

y0=[43994;0;1;5;0;0;0;0;6]; % 初始状态
k0=[7,5,0.5,0.58,0.2,0.94,0.7,0.2,0.03]; %猜测参数初值

lb=[2 1 0.01 0.01 0.001 0.001 0.001 0.001 0.001]; % 参数下限
ub=[30 30 1 1 1 1 1 1 1];  % 参数上限

代码如下：（如果不能运行可能是复制文本有错误，附件备用代码.m文件）
clear;clc;close all;
format long

%方程9真实数据
ydata1 = [6,  12,  19,  25,  31,  38,  44,  60,  80,  131,  131,  259,  467,  688,  776 ...,
1776,  1460,  1739,  1984,  2101,  2590,  2827,  3233,  3892,  3697,  3151 ...,
3387,  2653,  2984,  2473,  2022,  1820,  1998,  1506,  1278,  2051,  1772 ...,
1891,  399,  894,  397,  650,  415,  518,  412,  439,  441,  435,  579,  206 ...,
130,  120,  143,  146,  102,  46,  45,  20,  31,  26,  11,  18,  27,  29,  39,  39]';
%方程8真实数据

ydata2= [0,  0,  0,  0,  0,  0,  0,  0,  4,  4,  4,  8,  15,  15,  25,  26,  26 ...,
38,  43,  46,  45,  57,  64,  66,  73,  73,  86,  89,  97,  108,  97,  254 ...,
121,  121,  142,  106,  106,  98,  115,  118,  109,  97,  150,  71,  52,  29 ...,
44,  37,  35,  42,  31,  38,  31,  30,  28,  27,  23,  17,  22,  11,  7,  14 ...,
10,  14,  13,  13]';

n=size(ydata1,1);%获取数据长度
tspan=1:1:n;% size=n*1

N = 11000000/250;
y0=[43994;0;1;5;0;0;0;0;6]; % 初始状态
k0=[7,5,0.5,0.58,0.2,0.94,0.7,0.2,0.03]; %猜测参数初值

lb=[2 1 0.01 0.01 0.001 0.02 0.02 0.001 0.001]; % 参数下限
ub=[30 30 1 1 1 1 1 1 1];  % 参数上限

%% ---------------------------------------------------------------------

% 使用函数lsqcurvefit()进行参数估计
options = optimset('MaxFunEvals',5000,'MaxIter',2000*12);
[k,resnorm]=lsqcurvefit(@(k,tspan)model(k,tspan,y0),k0,tspan,[ydata2 ydata1],lb,ub,options);

%% ----------------绘制结果图-----------------------------------------------------
[t,Y]=ode45(@(t,y)rigid(t,y,k),tspan,y0);
% 绘制y9人数随天数变化图以及模拟变化图
figure(1)
subplot(1,2,1);
plot(t,Y(:,9),'k',tspan,ydata1,'-g'),legend('模拟值','真实数据','Location','best');
xlabel('时间');ylabel('确诊人数');

% 绘制y8人数随天数变化图以及模拟变化图
subplot(1,2,2);
plot(t,Y(:,8),'k',tspan,ydata2,'-r'),legend('模拟值','真实数据','Location','best');
xlabel('时间');ylabel('人数');

clearvars -except k %只保留k

%% ---------------------------------------------------------

function p=model(k,tspan,y0) % 目标函数
[~,Y]=ode45(@(t,y)rigid(t,y,k),tspan,y0);% 调用语句
p=Y(:,8:9);
end

%% ---------------------------------------------------------

function dydt=rigid(~,y,k) % 微分方
N = 11000000/250;
dydt = zeros(9,1);

dydt(1)=-k(1)*(y(3)/N)*y(1)-k(2)*(y(4)/N)*y(1);
dydt(2)=k(1)*(y(3)/N)*y(1)+k(2)*(y(4)/N)*y(1)-k(3)*y(2);
dydt(3)=k(3)*k(4)*y(2)-k(7)*y(3);
dydt(4)=k(3)*k(5)*y(2)-k(6)*y(4);
dydt(5)=k(3)*(1-k(4)-k(5))*y(2);
dydt(6)=k(6)*y(4)+k(7)*y(3)-k(8)*y(6)-k(9)*y(6);
dydt(7)=k(8)*y(6);
dydt(8)=k(9)*y(6);
dydt(9)=dydt(3)+dydt(4)+dydt(6);
end

大小的小 · 发表于 2021-2-24 14:12

你的实验数据明显有几个毛刺或者说疑似异常值（outlier），现实中的实际数据的确不会完美地全部分布在模型对应的曲线上，但如果严重偏离，那么这种数据究竟是否仍旧有效，都是值得探讨的。没有对数据的有效性做判断就一股脑把原始数据全扔进去，寄望于模型能匹配出合适的参数，这种做法值得商榷。

cichishia · 发表于 2021-2-25 13:35

一楼正解！

zzz.dan · 发表于 2021-2-25 13:47

帐号		自动登录	找回密码
密码			注册

ode45和lsqcurvefit结合，对微分方程组的参数求最优解

EDA365欢迎您登录！

浏览过的版块

推荐内容 /1