EDA365电子论坛网

标题: ode45和lsqcurvefit结合，对微分方程组的参数求最优解 [打印本页]
作者: IBB-EUT    时间: 2021-2-24 13:48
标题: ode45和lsqcurvefit结合，对微分方程组的参数求最优解
采用ode45和lsqcurvefit结合的方法，对微分方程组的参数求最优解，效果如图，第二张表不理想，
3 u9 @4 }& Q+ {: Y# y. S! t! B, `感觉只获得了局部最优解。如何优化，不知有没有大神能否提供有效解决方法。
) p" M+ J6 I4 `) ?1 ~( J( d+ D, ~/ Q. x8 |" G$ p

  G% P% b$ Y% V微分方程：; g5 g7 J1 C, N
N = 11000000/250;
5 n9 j- I2 v# mdydt(1)=-k(1)*(y(3)/N)*y(1)-k(2)*(y(4)/N)*y(1);
2 e- O# T: X4 y/ y+ `/ P- Bdydt(2)=k(1)*(y(3)/N)*y(1)+k(2)*(y(4)/N)*y(1)-k(3)*y(2);7 T. E- A* B" U  P, n4 j0 l
dydt(3)=k(3)*k(4)*y(2)-k(7)*y(3);. x! n5 D' s9 s* z( M
dydt(4)=k(3)*k(5)*y(2)-k(6)*y(4);, e  b4 @# s. @  m! D. ]! g- `, E
dydt(5)=k(3)*(1-k(4)-k(5))*y(2);0 T) J; G" g3 a1 x* v0 j$ R/ y
dydt(6)=k(6)*y(4)+k(7)*y(3)-k(8)*y(6)-k(9)*y(6);
; u' e) g) v- }dydt(7)=k(8)*y(6);. ~& u4 _; G) I5 Y) e+ @0 b
dydt(8)=k(9)*y(6);0 i  S/ v; C3 D8 o4 e& B! M$ o$ @
dydt(9)=dydt(3)+dydt(4)+dydt(6);
0 W+ E2 ^9 g/ |& @
) U/ V" _# @. k2 N0 ~3 Oy0=[43994;0;1;5;0;0;0;0;6]; % 初始状态
% `7 V2 C( h. \9 d. |k0=[7,5,0.5,0.58,0.2,0.94,0.7,0.2,0.03]; %猜测参数初值8 i9 f9 x7 v7 j9 K# {
! p1 ^  Z8 E- ], N) _2 Z/ W
lb=[2 1 0.01 0.01 0.001 0.001 0.001 0.001 0.001]; % 参数下限! Z7 o  ^5 L+ V
ub=[30 30 1 1 1 1 1 1 1];  % 参数上限
: {+ Y5 M1 ~3 w# Y- |' a1 V, C" m6 m4 j. A1 ~1 U, r0 S: H' l
代码如下：（如果不能运行可能是复制文本有错误，附件备用代码.m文件）; T9 E& ]+ d% f
clear;clc;close all;
* g% O0 s9 L9 C2 Nformat long) n# d- l2 W: X* j8 E

* ]3 D$ @2 m" u: j6 G4 O%方程9真实数据
8 R" `% Q- {: d% O. eydata1 = [6,  12,  19,  25,  31,  38,  44,  60,  80,  131,  131,  259,  467,  688,  776 ...,
, W# e  O9 m5 ^$ Q1776,  1460,  1739,  1984,  2101,  2590,  2827,  3233,  3892,  3697,  3151 ...,
+ R, ?) m1 e2 p8 ]3 f. o3387,  2653,  2984,  2473,  2022,  1820,  1998,  1506,  1278,  2051,  1772 ...,
& w/ E* n: w3 s- J: ^; l  z$ I1891,  399,  894,  397,  650,  415,  518,  412,  439,  441,  435,  579,  206 ...,
) q3 d, `$ n/ z2 D2 C( ^# S# G130,  120,  143,  146,  102,  46,  45,  20,  31,  26,  11,  18,  27,  29,  39,  39]';3 i$ M# n. y/ E# a
%方程8真实数据4 U8 X# @% |, W2 A  j, k/ B' N

5 X1 N  C( A3 o( Q& Uydata2= [0,  0,  0,  0,  0,  0,  0,  0,  4,  4,  4,  8,  15,  15,  25,  26,  26 ...,
; Z  l# W: \, q* j3 U* f9 r& |    38,  43,  46,  45,  57,  64,  66,  73,  73,  86,  89,  97,  108,  97,  254 ...,4 i" F. x' o9 \4 G0 {
    121,  121,  142,  106,  106,  98,  115,  118,  109,  97,  150,  71,  52,  29 ...,
, N  K- h. {! c% J, X7 [! [    44,  37,  35,  42,  31,  38,  31,  30,  28,  27,  23,  17,  22,  11,  7,  14 ...,6 N3 |! ?- P3 [. @
    10,  14,  13,  13]';4 }. K" m: U: F5 q4 s% I

: M" M4 `; w9 a7 J" q8 i5 G, k; sn=size(ydata1,1);%获取数据长度
8 b' o8 I9 X; p9 v9 ^/ s9 itspan=1:1:n;% size=n*1
( R9 k1 F4 ^  T! x' J3 `5 y, k% R7 i
N = 11000000/250;' o9 n$ J( ~3 N9 S
y0=[43994;0;1;5;0;0;0;0;6]; % 初始状态
, A" h* ?/ v2 q& N; ak0=[7,5,0.5,0.58,0.2,0.94,0.7,0.2,0.03]; %猜测参数初值
% R7 c+ p6 x; ^) Z
4 C) T1 o- T, L" o8 P& mlb=[2 1 0.01 0.01 0.001 0.02 0.02 0.001 0.001]; % 参数下限0 j* x: ^9 E$ ~( `- M
ub=[30 30 1 1 1 1 1 1 1];  % 参数上限! h1 O; v7 ]) h# Y0 y

- F$ l% T6 ]3 r; _%% ---------------------------------------------------------------------2 F- U) d7 m9 A0 x- s
6 h9 `7 H( K/ K: k' X7 ]
% 使用函数lsqcurvefit()进行参数估计
7 i3 }" l3 B/ |- P* i( Z1 ^options = optimset('MaxFunEvals',5000,'MaxIter',2000*12);
% }" Q3 Z0 g! V9 ]+ o[k,resnorm]=lsqcurvefit(@(k,tspan)model(k,tspan,y0),k0,tspan,[ydata2 ydata1],lb,ub,options);
; N1 X' u9 \) l# M- {5 P! }0 o) z& z
%% ----------------绘制结果图-----------------------------------------------------
! Y3 ^: e3 \, ?9 t[t,Y]=ode45(@(t,y)rigid(t,y,k),tspan,y0);& j; w6 X9 S7 S4 `6 t9 _; |7 V
% 绘制y9人数随天数变化图以及模拟变化图7 L/ Y) X6 a' H6 G$ B5 I: B
figure(1)
  s1 _! m6 L# k8 w8 t( Tsubplot(1,2,1);
' _/ t8 X! U6 G7 |3 y0 h% splot(t,Y(:,9),'k',tspan,ydata1,'-g'),legend('模拟值','真实数据','Location','best');
9 l. b3 ~3 i3 X6 E, ^* c( o6 x0 Lxlabel('时间');ylabel('确诊人数');
: Y2 B. X& z2 {+ T7 J
5 s' F& T0 r, Q% 绘制y8人数随天数变化图以及模拟变化图  }2 b5 T* E' ^& {
subplot(1,2,2);3 H" b0 \. N/ v- Q& i1 J- D. D
plot(t,Y(:,8),'k',tspan,ydata2,'-r'),legend('模拟值','真实数据','Location','best');
+ Y1 p8 H' G& p4 S4 P  q( M/ B3 ~xlabel('时间');ylabel('人数');0 A( E* Z% I, l
9 ?/ K# J0 Y( X$ e  Z% Z4 r7 p
clearvars -except k %只保留k
9 ^- U" V  B' `: d' f6 k
" t& q: d/ c3 ?6 d%% ---------------------------------------------------------' a* }2 G% J1 C& [# [) r# }- {

7 p" a) X1 d) V& [function p=model(k,tspan,y0) % 目标函数* ~* K) E* t% R. ^
[~,Y]=ode45(@(t,y)rigid(t,y,k),tspan,y0);% 调用语句0 z2 f& }1 I$ Z; D5 V1 M
p=Y(:,8:9);
2 ?+ a  }! q4 `" k" V- b; q$ b' @5 Eend0 J# t; f& S* l% u& l: |: C* J/ i

4 N1 |" t7 F1 ]3 |! A8 \3 c%% ---------------------------------------------------------; a: L& T. m7 t; L
" v4 ~6 R+ s, C
function dydt=rigid(~,y,k)    % 微分方0 F, ?: o3 \3 w' J, H/ ^1 h$ \9 R7 Q
N = 11000000/250;$ U" Z& k3 v. J. E
dydt = zeros(9,1);
1 R  _$ [9 Z/ Y" b* z: [' s  f! L) M' _9 N( E
dydt(1)=-k(1)*(y(3)/N)*y(1)-k(2)*(y(4)/N)*y(1);5 S0 \- Y4 ?7 B* f- m( b0 O( v* \/ x1 a* {
dydt(2)=k(1)*(y(3)/N)*y(1)+k(2)*(y(4)/N)*y(1)-k(3)*y(2);. @( a9 `+ V. r) ?% s, `  u+ ^5 _
dydt(3)=k(3)*k(4)*y(2)-k(7)*y(3);
5 _& B. \( I8 M/ U7 qdydt(4)=k(3)*k(5)*y(2)-k(6)*y(4);' P( d! E0 }) `( t1 H' }& q
dydt(5)=k(3)*(1-k(4)-k(5))*y(2);9 Y# n2 U* G6 f% Z9 K! p# `
dydt(6)=k(6)*y(4)+k(7)*y(3)-k(8)*y(6)-k(9)*y(6);
3 t! S' M4 L+ N! L( j# jdydt(7)=k(8)*y(6);
+ J) q" S% n( `9 I2 v) B( bdydt(8)=k(9)*y(6);
) d! J- y6 S& V* v9 Pdydt(9)=dydt(3)+dydt(4)+dydt(6);8 l: g5 U0 S. t4 g
end
: o& _& ]) u7 J; L& k4 V5 P& W% F$ r0 u7 n

3 u  b9 r/ o7 S# t% `
作者: 大小的小    时间: 2021-2-24 14:12
你的实验数据明显有几个毛刺或者说疑似异常值（outlier），现实中的实际数据的确不会完美地全部分布在模型对应的曲线上，但如果严重偏离，那么这种数据究竟是否仍旧有效，都是值得探讨的。没有对数据的有效性做判断就一股脑把原始数据全扔进去，寄望于模型能匹配出合适的参数，这种做法值得商榷。
作者: cichishia    时间: 2021-2-25 13:35
一楼正解！
作者: zzz.dan    时间: 2021-2-25 13:47




欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/) Powered by Discuz! X3.2