% Algorithm to estimate optimal model paramteres of a differential
% equation. PJH, 2017
% The parameters should be estimated in the way that all measurement data
% is (loading at different points in time) taken into account.

clear all; clc;

% Time-Arrays, corresponding to measurement Data
t_xls(:,1) = xlsread('cc.xlsx','curves_1xls','A4:A675');
t_xls(:,2) = xlsread('cc.xlsx','curves_1xls','A4:A675');
t_xls(:,3) = xlsread('cc.xlsx','curves_1xls','A4:A675');
t_xls(:,4) = xlsread('cc.xlsx','curves_1xls','A4:A675');
t_xls(:,5) = xlsread('cc.xlsx','curves_1xls','A4:A675');
t_xls(:,6) = xlsread('cc.xlsx','curves_1xls','A4:A675');

% strain (epsilon) data from measurement
eps_xls(:,1) = xlsread('cc.xlsx','curves_1xls','D4:D675');
eps_xls(:,2) = xlsread('cc.xlsx','curves_1xls','E4:E675');
eps_xls(:,3) = xlsread('cc.xlsx','curves_1xls','F4:F675');
eps_xls(:,4) = xlsread('cc.xlsx','curves_1xls','G4:G675');
eps_xls(:,5) = xlsread('cc.xlsx','curves_1xls','H4:H675');
eps_xls(:,6) = xlsread('cc.xlsx','curves_1xls','I4:I675');

% applied unit-load (compression)
sigma0 = -1.0;

% adapt time- and strain-vectors, so that they only contain points in time,
% when load is applied (= when strain-data epsilon != 0)
for ii=1:size(eps_xls,2)
    ctr = 1;
    for nn=1:size(eps_xls,1)        
       if eps_xls(nn,ii)~=0
           epsvec_adap(ctr,ii) = eps_xls(nn,ii);
           timevec_adap(ctr,ii) = t_xls(nn,ii);
           ctr = ctr+1;
       end
    end
end

[maxval, index] = max(timevec_adap(:,size(timevec_adap,2)));
[row, col] = ind2sub(size(timevec_adap(:,size(timevec_adap,2))),index);
timevec = timevec_adap(1:row,:);
epsvec = epsvec_adap(1:row,:);

% to obtain a better fit, generate more data-points by linear interpolation
dbl_dt_early = 0.001;
dbl_dt_late = 0.001;

for ii=1:size(timevec,2)
  xi1 = timevec(1,ii):dbl_dt_early:(min(timevec(:,ii))+4);
  xi2 = (min(timevec(:,ii))+4):dbl_dt_late:max(timevec(:,ii));
  xi2(1) = [];
  xi = [xi1, xi2];  
  yi = interp1(timevec(:,ii), epsvec(:,ii), xi);
  
  tref(:,ii) = xi';
  epsref(:,ii) = yi';
end

timevec = tref;
epsvec = epsref;

% Initial conditions: values taken from literature
pf    = [28000;-1.95;-1.25;];
pE1   =  2900;
pd11  = [247;0.643];
pd12  = [3239; 0.643];
pd13  = [858; 0.643];
pd14  = [49394016; 0.643];

pf4 =  910;  % --> 0 in calibration
pd31 = 100;  % --> 0 in calibration
pd32 = 690;  % --> 0 in calibration
pd33 = 1400; % --> 0 in calibration
pd34 = 0;    % --> 0 in calibration

% optimization: only calibrate model to measurement data where t >= 2d
% -------------------------------------------------------------------------

% a) first try: do not set any boundaries
% ######

% create vector with initial guess for parameters
InitialGuess1 = [pf(1); pf(2); pf(3); pE1; pd11(1); pd11(2); pd12(1); ...
    pd13(1); pd14(1)]; 
fit_lb1 = [];
fit_ub1 = [];

% constant parameters that have to remain zero (they'd only affect
% solution, if curves t < 2d are used. but this is not the case here.
fixvec1 = [0;0;0;0;0];

calcParam1 = lsqnonlin(@(PAR) objectiveFunction_lsq_p1(PAR, fixvec1, ...
    timevec(:,3:end), epsvec(:,3:end),sigma0), InitialGuess1, fit_lb1, fit_ub1);

% b) second try: set some boundaries (acc. to dimensions of initial conds)
% ######
InitialGuess2 = InitialGuess1;
fit_lb2 = [1000;  -3;  -3;  1000; 100;  0;   1000;  500;  1e5]; 
fit_ub2 = [34100;  0;   0; 34100; 500;  2;  40000; 1000;  1e8];
fixvec2 = [0;0;0;0;0];

calcParam2 = lsqnonlin(@(PAR) objectiveFunction_lsq_p1(PAR, fixvec2, ...
    timevec(:,3:end), epsvec(:,3:end),sigma0), InitialGuess2, fit_lb2, fit_ub2);

% Alternative: use Multistart for global minimization: no real improvement
% of result noticeable
% problem = createOptimProblem('lsqnonlin','x0',InitialGuess2, ...
%         'lb', fit_lb2, 'ub', fit_ub2, ...
%         'objective',@(PAR) objectiveFunction_lsq_p1(PAR, fixvec2, ...
%         timevec(:,3:end), epsvec(:,3:end), sigma0));
% mstart = MultiStart('PlotFcns', @gsplotbestf);
% [x, fval, eflag, output, manymins] = run(mstart, problem, 50);
% calcParam2 = manymins(1, 1).X;


% c) third thing: set boundaries, use arbitrary initial conditions (some-
%    where between lb3 and ub3
% ######
InitialGuess3 = [20000;  -1.5;  -1.5;  15000; 250;   1;   20000; 800; 1e5]; 
fit_lb3 = fit_lb2;
fit_ub3 = fit_ub2;
fixvec3 = [0;0;0;0;0];

calcParam3 = lsqnonlin(@(PAR) objectiveFunction_lsq_p1(PAR, fixvec3, ...
    timevec(:,3:end), epsvec(:,3:end),sigma0), ... 
    InitialGuess3, fit_lb3, fit_ub3);

% Test of estimated Parameters
% -------------------------------------------------------------------------
for ii=3:size(timevec,2)

 % no boundary conditions used
 solu_Fit1(:,ii) = calc_epsilon(timevec(:,ii), sigma0, [0 0 0 0], ...
     [calcParam1(1); calcParam1(2); calcParam1(3)], calcParam1(4), ...
     [calcParam1(5); calcParam1(6)], [calcParam1(7); calcParam1(6)], ...
     [calcParam1(8); calcParam1(6)], [calcParam1(9); calcParam1(6)], ...
      0, 0, 0, 0, 0, 1 );
  
 % with boundaries
 solu_Fit2(:,ii) = calc_epsilon(timevec(:,ii), sigma0, [0 0 0 0], ...
     [calcParam2(1); calcParam2(2); calcParam2(3)], calcParam2(4), ...
     [calcParam2(5); calcParam2(6)], [calcParam2(7); calcParam2(6)], ...
     [calcParam2(8); calcParam2(6)], [calcParam2(9); calcParam2(6)], ...
      0, 0, 0, 0, 0, 1 );
  
 % with boundaries and arbitrary initial conds
  solu_Fit3(:,ii) = calc_epsilon(timevec(:,ii), sigma0, [0 0 0 0], ...
     [calcParam3(1); calcParam3(2); calcParam3(3)], calcParam3(4), ...
     [calcParam3(5); calcParam3(6)], [calcParam3(7); calcParam3(6)], ...
     [calcParam3(8); calcParam3(6)], [calcParam3(9); calcParam3(6)], ...
      0, 0, 0, 0, 0, 1 );
 
 % result with initial conditions
 solu_Orig(:,ii) = calc_epsilon(timevec(:,ii), sigma0, [0 0 0 0], ...
      pf, pE1, pd11, pd12, pd13, pd14, pf4, pd31, pd32, pd33, pd34, 1 ); 
  
end

epsvec = [ zeros(1,size(epsvec,2)); epsvec];
solu_Fit1 = [ zeros(1,size(solu_Fit1,2)); solu_Fit1];
solu_Fit2 = [ zeros(1,size(solu_Fit2,2)); solu_Fit2];
solu_Fit3 = [ zeros(1,size(solu_Fit3,2)); solu_Fit3];
solu_Orig = [ zeros(1,size(solu_Orig,2)); solu_Orig];
timevec = [ timevec(1,1:size(timevec,2)); timevec];

% plot results
% -------------------------------------------------------------------------
figure;
plot(timevec(:,3:6), epsvec(:,3:end), 'k', 'Linewidth', 2);
hold on;
plot(timevec(:,3:6), solu_Orig(:,3:end), 'r', 'Linewidth', 1.0, 'LineStyle', '--');
plot(timevec(:,3:6), solu_Fit1(:,3:end), 'b', 'Linewidth', 3.0, 'LineStyle', ':');
plot(timevec(:,3:6), solu_Fit2(:,3:end), 'm', 'Linewidth', 1.5, 'LineStyle', '-.');
plot(timevec(:,3:6), solu_Fit3(:,3:end), 'g', 'Linewidth', 1.5, 'LineStyle', '-.');
hold off;
