-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain_unified_SRSA.m
140 lines (129 loc) · 4.77 KB
/
main_unified_SRSA.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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
%% Author : TAO ZHANG * [email protected] *
% Created Time : 2022-10-01 08:58
% Last Revised : TAO ZHANG ,2023-03-19
% Remark : main program: fractional order unified system
% Model ref: Parameter Estimation of Fractional Chaotic Systems
% Based on Stepwise Integration and Response Sensitivity Analysis
clear;clc;close all;
addpath('./SRSA');
global parameter_a x_cal_data Tdata x_cal_iden
%% Model setting and data
parameter_a=[0.3;0.7];
ini_parameter_a=parameter_a;
parameter_a_judge=@(parameter_a)(min(parameter_a(1))>0);
load simple_fre_data.mat; % load observed data --- Tdata and Xdata
gammaT=1.414;rhob=0.5; % parameter for trust-region algorithm
parameter_a_record=parameter_a;
TR_record=[]; % recording the parameter_a during trust region
x_cal_start=x_cal_data(1:end-1,:);
x_cal_end=x_cal_data(2:end,:);
Tstart=Tdata(1:end-1);
Tend =Tdata(2:end);
Psize=size(x_cal_data);
Psize=[Psize(1),9];
%% Response sensitivity iteration
Nmax=1000; % maximum number for response sensitivity iteration
Ntr=20; % maximum number for trust region iteration
%% response sensitivity Solution by ode45
NT=length(x_cal_data(:,1))-1;
for iii=1:Nmax
tic;
% compute response and response sensitivity for each incremental
Etol=1e-10; % Relative error tolerance for convergence of the RS algorithm
%%
x_cal_iden=zeros(Psize(1),Psize(2));
x_cal_iden(1,1:3)=x_cal_data(1,1:3);
for jjj=1:NT
x_cal_iden_j=cal_unified_step(parameter_a,[Tstart(jjj),Tend(jjj)],[x_cal_start(jjj,1:3),zeros(1,3)]',jjj);
x_cal_iden(jjj+1,:)=x_cal_iden_j;
end
%% SSSÏìÓ¦ÁéÃô¶È¾ØÕó
SSS1=[x_cal_iden(2:end,4);x_cal_iden(2:end,5);x_cal_iden(2:end,6)];
SSS2=[x_cal_iden(2:end,7);x_cal_iden(2:end,8);x_cal_iden(2:end,9)];
SSS=[SSS1,SSS2];
% determine initial lambda by L-curve method
dR1=x_cal_end(:,1)-x_cal_iden(2:end,1);
dR2=x_cal_end(:,2)-x_cal_iden(2:end,2);
dR3=x_cal_end(:,3)-x_cal_iden(2:end,3);
dR=[dR1;dR2;dR3];
[U,s,V]=csvd(SSS);
lambda_inverse=l_curve(U,s,dR);
atemp=parameter_a;
% trust-region algorithm
for trust=1:Ntr
da=tikhonov(U,s,V,dR,lambda_inverse);
if ~parameter_a_judge(atemp+da)
lambda_inverse=lambda_inverse*gammaT; % update of lambda
continue;
end
da=tikhonov(U,s,V,dR,lambda_inverse);
parameter_a=atemp+da;
%% new parameter_a=atemp+da
x_cal_da=zeros(Psize(1),Psize(2));
x_cal_da(1,1:3)=x_cal_data(1,1:3);
for jjj=1:NT
x_cal_da_j=cal_unified_step(parameter_a,[Tstart(jjj),Tend(jjj)],[x_cal_start(jjj,1:3),zeros(1,6)]',jjj);
x_cal_da(jjj+1,:)=x_cal_da_j;
end
dR1temp=x_cal_end(:,1)-x_cal_da(2:end,1);
dR2temp=x_cal_end(:,2)-x_cal_da(2:end,2);
dR3temp=x_cal_end(:,3)-x_cal_da(2:end,3);
dRtemp=[dR1temp;dR2temp;dR3temp];
LdR=SSS*da-dR;
rhos=(dR'*dR-dRtemp'*dRtemp)/(dR'*dR-LdR'*LdR); % agreement indicator
if rhos>=rhob
break;
end
lambda_inverse=lambda_inverse*gammaT;
end
da=tikhonov(U,s,V,dR,lambda_inverse);
parameter_a=atemp+da;
tolt=norm(da)/norm(parameter_a);
parameter_a_record=[parameter_a_record,parameter_a];
TR_record=[TR_record;lambda_inverse];
parameter_a
if tolt<=Etol
break;
end
every_a(iii).parameter_a=parameter_a;
iii
toc;timejishi(iii)=toc;
end
% system real parameter
real_parameter_a=[0.5;0.95];
relative_error=(parameter_a-real_parameter_a)./real_parameter_a*100
a1=zeros(1,iii);
a2=a1;
relative_error_1=a1;
relative_error_2=a1;
relative_error_1(1,1)=(ini_parameter_a(1)-real_parameter_a(1))/real_parameter_a(1)*100;
relative_error_2(1,1)=(ini_parameter_a(2)-real_parameter_a(2))/real_parameter_a(2)*100;
a1(1)=ini_parameter_a(1);
a2(1)=ini_parameter_a(2);
for i=2:iii
a=every_a(i-1).parameter_a;
a1(i)=a(1);
a2(i)=a(2);
relative_error_1(i)=(a(1)-real_parameter_a(1))/real_parameter_a(1)*100;
relative_error_2(i)=(a(2)-real_parameter_a(2))/real_parameter_a(2)*100;
end
figure;
plot(1:1:iii,relative_error_1,'r-','LineWidth',1.5);
hold on;
plot(1:1:iii,relative_error_2,'b--','LineWidth',1.5);
set(gca,'LineWidth',1.5);set(gca,'FontSize',15);
figure;
plot(1:1:iii,a1,'r-','LineWidth',1.5);
hold on;
plot(1:1:iii,a2,'b--','LineWidth',1.5);
set(gca,'LineWidth',1.5);set(gca,'FontSize',15);
h1=legend('$$a$$','$$q$$');
set(h1,'Interpreter','latex','FontSize',15);
set(gca,'FontName','Times New Roman','FontSize',15,'LineWidth',1.5);
ini_parameter_a
timejishi_record(1)=timejishi(1);
for i=1:iii-2
timejishi_record(i+1)=timejishi_record(i)+timejishi(i);
end
figure;
plot(1:iii-1,timejishi_record,'r*-');