-
Notifications
You must be signed in to change notification settings - Fork 0
/
para.m
64 lines (60 loc) · 1.84 KB
/
para.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
function error = para(value)
warning('off');
actual = readmatrix('./data/currentconfirmed.csv');
actual_used = actual(1:135)';
value = cell2mat(value);
i_I = value(1);
i_Q =value(2);
beta_iq =value(3);
beta_ir =value(4);
beta_qr =value(5);
gamma_2 =value(6);
beta_bd =value(7);
phi = 0.014 ;
omega = 0.00020602739;
mu = 0.0001967123;
K = 1500000000;
r = omega - mu;
ddefun = @(t,x,Z)[
(omega-phi*r*x(6)/K)*x(6)-(beta_bd*(Z(3,1)+Z(4,1))/Z(6,1))*x(1)-(mu+(1-phi)*r*x(6)/K)*x(1);
(beta_bd*(Z(3,1)+Z(4,1))/Z(6,1))*x(1) - (mu+(1-phi)*r*x(6)/K)*x(2) - (Z(3,1)*i_I+Z(4,1)*i_Q)*x(2)/(Z(3,1)+Z(4,1));
(Z(3,1)*i_I+Z(4,1)*i_Q)*x(2)/(Z(3,1)+Z(4,1)) - beta_ir*x(3) - beta_iq*x(3) - (mu+(1-phi)*r*x(6)/K)*x(3);
beta_iq*x(3) - beta_qr*x(4) - (mu+(1-phi)*r*x(6)/K)*x(4) - gamma_2*x(4);
beta_qr*x(4) + beta_ir*x(3) - (mu+(1-phi)*r*x(6)/K)*x(5);
(omega-phi*r*x(6)/K)*x(6)-(mu+(1-phi)*r*x(6)/K)*x(6);
];
options = ddeset('RelTol',1,'AbsTol',10);
sol = dde23(ddefun,[1 1 1 1 1 1],[1154633644,288658411,1218,776,34,1443497378],[1,135]);
tint = linspace(1,135,135);
l = length(sol.y(4,:));
x = linspace(1,l,l);
% if l < 135
% solv = sol.y(4,:);
% sols = pchip(x,solv,tint);
% end
%
% if l == 135
% sols = sol.y(4,:);
% end
%
% if l > 135
% sols = 1:135;
% n = idivide(int32(l),int32(135),'round');
% i = 1;
% while i < 135
% sols(i) = sol.y(4,i);
% i = i+n;
% end
% if length(sols) < 135
% sols(135) = sol.y(4,l);
% end
% end
if l == 135
sols = sol.y(4,:);
else
solv = sol.y(4,:);
sols = pchip(x,solv,tint);
end
sols = round(sols);
error = sqrt(mse(sols,actual_used));
end