% Open model
modelName = 'ee_supercapacitor_identification';
open_system(modelName);
set_param(find_system('ee_supercapacitor_identification','FindAll', 'on','type','annotation','Tag','ModelFeatures'),'Interpreter','off')

% Configure the supercapacitor using known parameter values.
%
% Supercapacitor block Cell Characteristics parameters
Kv = 190;                   % Voltage-dependent capacitor gain
R = [2.5e-3 0.9 5.2];       % Fixed resistances, [R1 R2 R3]
C = [270 100 220];          % Fixed capacitances, [C1 C2 C3]

% Store a copy of the known parameter values for an eventual comparison to
% parameter values identified using the methodology described in [1].
Kv_known = Kv;
R_known = R;
C_known = C;

% Specify input current by configuring the Step block.
%
% Step block parameters
stepTime = [40 1900 1917];  % Step time
initialValue = [28 0 0];    % Initial value
finalValue = [0 -25 25];    % Final value

% Specify the physical characteristics of the supercapacitor.
%
% Supercapacitor block Cell Characteristics and Configuration parameters
R_discharge = 9e3;          % Self-discharge resistance
N_series = 1;               % Number of series cells
N_parallel = 1;             % Number of parallel cells

% Specify the simulation duration
stopTime = 2100;

% Simulate and store voltage and current waveform data
sim(modelName);
t = simlog_ee_supercapacitor_identification.Sensing_current.It.i.series.time;
i = simlog_ee_supercapacitor_identification.Sensing_current.It.i.series.values;
v = simlog_ee_supercapacitor_identification.Sensing_voltage.V.V.series.values;

% Event n=1
t1 = 20e-3;
i1 = interp1(t,i,t1);
v1 = interp1(t,v,t1);

Ri = v1/i1;

% Extract charging data that interp1 can use to find time values rather
% than voltage values
[v3, v_max_idx] = max(v);
[v_charge, v_charge_idx] = unique(v(1:v_max_idx));
t_charge = t(v_charge_idx);

% Event n=2
delta_v = 50e-3;
v2 = v1+delta_v;
t2 = interp1(v_charge,t_charge,v2);
i2 = interp1(t,i,t2);

Ci0 = i2*(t2-t1)/delta_v;

% Event n=3
t3 = t(v_max_idx);

% Event n=4
delta_t = 20e-3;
t4 = t3+delta_t;
v4 = interp1(t,v,t4);
Qtotal = i1*(t4-t1);
Cq = Qtotal/v4;

Ci1 = (2/v4)*(Cq-Ci0);

% Event n=5
v5 = v4-delta_v;
discharge_idx = find(i<-20,1,'first');
t_discharge = t(v_max_idx:discharge_idx);
v_discharge = v(v_max_idx:discharge_idx);
t5 = interp1(v_discharge,t_discharge,v5);
delta_t = t5-t4;
Vci = v4-(delta_v/2);
Cdiff = Ci0+Ci1*Vci;
Rd = (v4-(delta_v/2))*delta_t/(Cdiff*delta_v);

% Event n=6
TypicallyRdTimesCd = 100;
t6 = t5 + 3*TypicallyRdTimesCd;
v6 = interp1(t,v,t6);
Cd = (Qtotal/v6)-(Ci0+((Ci1/2)*v6));

% Event n=7
v7 = v6-delta_v;
t7 = interp1(v_discharge,t_discharge,v7);
delta_t = t7-t6;
Vci = v6-(delta_v/2);
Cdiff = Ci0+Ci1*Vci;
Rl = (v6-(delta_v/2))*delta_t/(Cdiff*delta_v); % Rl value large

% Event n=8
t8 = 30*60;
v8 = interp1(t,v,t8);
Cl = (Qtotal/v8)-(Ci0+(v8*Ci1/2))-Cd; % Cl value too large

DataTable = table((1:8)',...
    [t1 t2 t3 t4 t5 t6 t7 t8]',...
    [v1 v2 v3 v4 v5 v6 v7 v8]',...
    'VariableNames',{'Event','Time','Voltage'}) %#ok<NOPTS>

% Create parameters for supercapacitor block using identified values and
% create copies for easy reference.
Kv = Ci1;
R = [Ri Rd Rl];
C = [Ci0 Cd Cl];
Kv_identified = Kv;
R_identified = R;
C_identified = C;

% Simulate the model and store voltage and current waveforms.
sim(modelName);
t_ = simlog_ee_supercapacitor_identification.Sensing_current.It.i.series.time;
v_ = simlog_ee_supercapacitor_identification.Sensing_voltage.V.V.series.values;

% Plot the data generated from simulating with known parameters on the
% same axis as the data generated from simulating with identified
% parameters.
figure;
plot(t,v,DataTable.Time,DataTable.Voltage,'o',t_,v_);
grid('on');
legend('Known parameters simulation',...
    'Measurements for parameter identification',...
    'Identified parameters simulation');

openfig('ee_supercapacitor_speresults');