% simulation of a work process for the model MAX stand

% Copyright 2023 Collaborative Research Centre Transregio 96 (CRC/TR 96)
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
% TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER
% OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%% load matrices
SA{1} = load('maxStandSA1_n25872m10q10.mat');
SA{2} = load('maxStandSA2_n39527m30q30.mat');
SA{3} = load('maxStandSA3_n13551m6q6.mat');
SA{4} = load('maxStandSA4_n4813m23q23.mat');

for j=1:4
    k=1;
    while k~=0
        fieldname = ['A' num2str(k)];
        if isfield(SA{j},fieldname)
            SA{j}.Ai{k} = SA{j}.(fieldname);
            k = k+1;
        else 
            k=0;
        end
    end
end

%% compute coductivity matrix A
% heat transfer coefficients may be varied for use as parametric system
alpha_conv=10;    % between stand and ambience
alpha_Boden=100;  % between stand and floor

kappa_ext=cell(4,1);
kappa_ext{1}=[alpha_Boden alpha_conv];
kappa_ext{2}=alpha_conv;
kappa_ext{3}=alpha_conv;
kappa_ext{4}=alpha_conv;

% compute (parameter dependent) system matrix A
A = cell(1,4);
for j=1:4
    A{j} = SA{j}.A;        
    for k=1:size(SA{j}.Ai,2)
        A{j} = A{j} - kappa_ext{j}(k)*SA{j}.Ai{k};
    end 
end
%%

% dt has to be in full seconds! (dt \in \mathbb{N})
t_start = 0*3600;           % start time
t_end   = 16*3600;          % end time
dt = 10;                    % time step size
N = (t_end-t_start)/dt;     % number of time steps

u = cell(1,4);
u_sim = cell(1,4);
y = cell(1,4);
T = cell(1,4);
T_orig = cell(1,4);

for j=1:4
    T{j} = zeros(size(SA{j}.A,1),1);
    T_orig{j}=zeros(size(SA{j}.C,1),N+1);
    y{j}=zeros(size(SA{j}.C,1),1);
    u_sim{j}=zeros(size(SA{j}.B,2),N);
end

t=t_start;

first = [1, 1, 1, 1];
t1 = tic;
for i = 1:N    
    
    % compute input u(t)
    
    % an average value of u(t) is computed for dt
    % (one segment of the guide rail and spindle has a length of ca. 0.04 m
    % -> with a feed speed of 0.04 m/s one scanning instance per second is
    % sufficient
    % for a feed speed of 0.833 m/s 20 scanning instances per second are 
    % necessary                 
    
    if t<=(4*3600) || t>(13*3600) 
       dt_u=1;           % time step size for new values of u(t)
    else
        if t<=(7*3600)
            dt_u=0.5;
        else
            if t<=(10*3600)
                dt_u=0.125;
            else
                if t<=(13*3600)
                    dt_u=0.05;
                end
            end
        end
    end
    
    step_anz=dt/dt_u;
    time_mesh=t:dt_u:(t+dt-dt_u);

    u=MAX_stand_input(time_mesh,y,alpha_conv,alpha_Boden); 
    
    if not(mod(i,360))
       fprintf(' Step %d / %d\n',i,N);
    end
    %disp(i);  
    for j = 1:4      
        u_sim{j}(:,i)=u{j};                  
        
        % compute temperature T(t)
        % for constant coefficient matrices and step sizes use precomputed
        % LU factorization
        if first(j)
           first(j) = 0;
           [L{j},U{j},P{j},Q{j}] = lu((SA{j}.E-dt*A{j})); %#ok
        end
        T{j} = Q{j}*(U{j}\(L{j}\(P{j}*(SA{j}.E*T{j}+dt*SA{j}.B*u{j}))));

        % compute output y(t)
        y{j} = SA{j}.C*T{j};
         
        % temperature in special areas/nodes of interest (for comparison, error analysis)
        T_orig{j}(:,i+1) = y{j};         
      
    end  
    
    t = t + dt;
    
end
tsim = toc(t1);