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

IBB-EUT

采用ode45和lsqcurvefit结合的方法，对微分方程组的参数求最优解，效果如图，第二张表不理想，
- b8 a" I# S+ W! k' J感觉只获得了局部最优解。如何优化，不知有没有大神能否提供有效解决方法。
& e- {. ^& _  N7 A* P$ v* Q! _
* @8 k" o7 v" K
微分方程：
N = 11000000/250;
( Z) V+ j: l# `0 Z( B* ddydt(1)=-k(1)*(y(3)/N)*y(1)-k(2)*(y(4)/N)*y(1);
7 `7 g$ T2 q8 X3 O% pdydt(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);
: j0 R* ]8 ~. f. R( y9 \; e1 kdydt(8)=k(9)*y(6);
dydt(9)=dydt(3)+dydt(4)+dydt(6);
. D5 r- W- ]! l
" F: \- w5 }- A8 z4 O! M$ v2 X/ Ky0=[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]; %猜测参数初值
' }4 r9 ~5 K( t4 G4 s: q
lb=[2 1 0.01 0.01 0.001 0.001 0.001 0.001 0.001]; % 参数下限
  N; B3 S# I4 h" h) @$ ^0 ?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 ...,
* W* R  r5 a( ~0 s7 [! a) t$ `. i1776,  1460,  1739,  1984,  2101,  2590,  2827,  3233,  3892,  3697,  3151 ...,
1 }1 D8 w6 ~9 [3387,  2653,  2984,  2473,  2022,  1820,  1998,  1506,  1278,  2051,  1772 ...,
1891,  399,  894,  397,  650,  415,  518,  412,  439,  441,  435,  579,  206 ...,
* `+ \. Y% A( R/ a( F2 w130,  120,  143,  146,  102,  46,  45,  20,  31,  26,  11,  18,  27,  29,  39,  39]';
: ~8 E* X) W6 d; Z%方程8真实数据
" W' {2 R6 {6 E! l2 N5 m% j
) h6 {; Y$ }+ {ydata2= [0,  0,  0,  0,  0,  0,  0,  0,  4,  4,  4,  8,  15,  15,  25,  26,  26 ...,
# n- }& A, T/ a+ U5 b    38,  43,  46,  45,  57,  64,  66,  73,  73,  86,  89,  97,  108,  97,  254 ...,
4 S1 u1 e+ I/ m! m8 W3 u    121,  121,  142,  106,  106,  98,  115,  118,  109,  97,  150,  71,  52,  29 ...,
  Y$ V9 H5 l% q; ?    44,  37,  35,  42,  31,  38,  31,  30,  28,  27,  23,  17,  22,  11,  7,  14 ...,
2 D& `7 o5 r1 x1 [- T- H0 N) ]    10,  14,  13,  13]';
7 L- K0 W7 z! ^1 R
n=size(ydata1,1);%获取数据长度
tspan=1:1:n;% size=n*1
4 T$ J4 X! _6 Z/ g) q) V5 q- B
: @4 {* o# p2 H6 lN = 11000000/250;
y0=[43994;0;1;5;0;0;0;0;6]; % 初始状态
3 l+ z) c8 w% z8 Qk0=[7,5,0.5,0.58,0.2,0.94,0.7,0.2,0.03]; %猜测参数初值
: [9 q1 H5 ~& R" l, V7 R
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];  % 参数上限

%% ---------------------------------------------------------------------
; R0 |/ y# s. `7 b- S
6 k, y' E: l8 F/ ?% v' Z6 u% 使用函数lsqcurvefit()进行参数估计
options = optimset('MaxFunEvals',5000,'MaxIter',2000*12);
[k,resnorm]=lsqcurvefit(@(k,tspan)model(k,tspan,y0),k0,tspan,[ydata2 ydata1],lb,ub,options);
( f  C, n5 [% o5 }: N
7 ]1 z: L3 c5 K& L, Z%% ----------------绘制结果图-----------------------------------------------------
  d9 B# _+ @$ _! D! _8 `7 T& K[t,Y]=ode45(@(t,y)rigid(t,y,k),tspan,y0);
% 绘制y9人数随天数变化图以及模拟变化图
figure(1)
" `  d/ P  |/ {0 Q" w; @% G% }9 ysubplot(1,2,1);
! d$ I- u1 H' N1 S4 a6 lplot(t,Y(:,9),'k',tspan,ydata1,'-g'),legend('模拟值','真实数据','Location','best');
9 b; q8 [# I9 K: K. o7 K1 S& [xlabel('时间');ylabel('确诊人数');

4 I8 }) d% q4 T6 r: ]% 绘制y8人数随天数变化图以及模拟变化图
subplot(1,2,2);
* L9 s  T* Y4 z. [* Nplot(t,Y(:,8),'k',tspan,ydata2,'-r'),legend('模拟值','真实数据','Location','best');
% e3 {$ ^3 x+ u' b: T% N* L) _9 Rxlabel('时间');ylabel('人数');

) o. v+ P! Y7 i" p( K/ |8 sclearvars -except k %只保留k

%% ---------------------------------------------------------
3 J9 s3 R0 x9 }7 n
1 ^9 q, J8 g" p4 g8 C9 A' dfunction p=model(k,tspan,y0) % 目标函数
, b, h4 R' \( k  B  {% d[~,Y]=ode45(@(t,y)rigid(t,y,k),tspan,y0);% 调用语句
/ K6 u" W+ ^$ y, j4 Y3 i& ?) cp=Y(:,8:9);
# T" e/ Y" f) T/ A9 m/ C/ @end

%% ---------------------------------------------------------
+ ^7 m6 a( m4 p6 {
function dydt=rigid(~,y,k)    % 微分方
1 x: S! s7 @( T; _3 L$ yN = 11000000/250;
6 @6 D7 r" p/ b7 u' ldydt = zeros(9,1);
$ f8 O+ o/ k2 X/ A5 B) F! j7 z( B' l
dydt(1)=-k(1)*(y(3)/N)*y(1)-k(2)*(y(4)/N)*y(1);
  F% `  x4 w6 k* I+ Q$ H' Odydt(2)=k(1)*(y(3)/N)*y(1)+k(2)*(y(4)/N)*y(1)-k(3)*y(2);
$ K2 v( r# [& q" O! g- Rdydt(3)=k(3)*k(4)*y(2)-k(7)*y(3);
5 B) G! q/ Q" n+ Z  c/ S6 Bdydt(4)=k(3)*k(5)*y(2)-k(6)*y(4);
* T+ l4 p6 y' S2 k$ A8 zdydt(5)=k(3)*(1-k(4)-k(5))*y(2);
4 v# F9 M0 H( {0 Rdydt(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);
8 @3 e; J0 q  S# e5 d9 L. Bdydt(9)=dydt(3)+dydt(4)+dydt(6);
end
8 o, p" e4 B# o2 _- G( R% I

; H0 W' S0 Q" {) e$ P1 z6 L* e% O

zzz.dan

cichishia

一楼正解！

大小的小

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