Nano-Hexapod - Test Bench

%% Load Identification Data
meas_data_lf = {};

for i = 1:6
    meas_data_lf(i) = {load(sprintf('mat/frf_data_exc_strut_%i_noise_lf.mat', i), 't', 'Va', 'Vs', 'de')};
    meas_data_hf(i) = {load(sprintf('mat/frf_data_exc_strut_%i_noise_hf.mat', i), 't', 'Va', 'Vs', 'de')};
end

%% Setup useful variables
% Sampling Time [s]
Ts = (meas_data_lf{1}.t(end) - (meas_data_lf{1}.t(1)))/(length(meas_data_lf{1}.t)-1);

% Sampling Frequency [Hz]
Fs = 1/Ts;

% Hannning Windows
win = hanning(ceil(1*Fs));

% And we get the frequency vector
[~, f] = tfestimate(meas_data_lf{1}.Va, meas_data_lf{1}.de, win, [], [], 1/Ts);

i_lf = f < 250; % Points for low frequency excitation
i_hf = f > 250; % Points for high frequency excitation

%% Coherence
coh_dvf = zeros(length(f), 6, 6);

for i = 1:6
    coh_dvf_lf = mscohere(meas_data_lf{i}.Va, meas_data_lf{i}.de, win, [], [], 1/Ts);
    coh_dvf_hf = mscohere(meas_data_hf{i}.Va, meas_data_hf{i}.de, win, [], [], 1/Ts);
    coh_dvf(:,:,i) = [coh_dvf_lf(i_lf, :); coh_dvf_hf(i_hf, :)];
end

%% DVF Plant (transfer function from u to dLm)
G_dvf = zeros(length(f), 6, 6);

for i = 1:6
    G_dvf_lf = tfestimate(meas_data_lf{i}.Va, meas_data_lf{i}.de, win, [], [], 1/Ts);
    G_dvf_hf = tfestimate(meas_data_hf{i}.Va, meas_data_hf{i}.de, win, [], [], 1/Ts);
    G_dvf(:,:,i) = [G_dvf_lf(i_lf, :); G_dvf_hf(i_hf, :)];
end

%% Coherence for the IFF plant
coh_iff = zeros(length(f), 6, 6);

for i = 1:6
    coh_iff_lf = mscohere(meas_data_lf{i}.Va, meas_data_lf{i}.Vs, win, [], [], 1/Ts);
    coh_iff_hf = mscohere(meas_data_hf{i}.Va, meas_data_hf{i}.Vs, win, [], [], 1/Ts);
    coh_iff(:,:,i) = [coh_iff_lf(i_lf, :); coh_iff_hf(i_hf, :)];
end

%% IFF Plant
G_iff = zeros(length(f), 6, 6);

for i = 1:6
    G_iff_lf = tfestimate(meas_data_lf{i}.Va, meas_data_lf{i}.Vs, win, [], [], 1/Ts);
    G_iff_hf = tfestimate(meas_data_hf{i}.Va, meas_data_hf{i}.Vs, win, [], [], 1/Ts);
    G_iff(:,:,i) = [G_iff_lf(i_lf, :); G_iff_hf(i_hf, :)];
end

save('matlab/mat/identified_plants_enc_struts.mat', 'f', 'Ts', 'G_iff', 'G_dvf')

%% Load data
load('identified_plants_enc_struts.mat', 'f', 'Ts', 'G_iff', 'G_dvf')

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', '2dof');

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');   io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dum'],  1, 'openoutput'); io_i = io_i + 1; % Force Sensors

Giff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Initialization of the Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', 'flexible');

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Encoders

Gdvf = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Tested bending dampings [Nm/(rad/s)]
cRs = [1e-3, 5e-3, 1e-2, 5e-2, 1e-1];

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Encoders

%% Idenfity the transfer function from actuator to encoder for all bending dampins
Gs = {zeros(length(cRs), 1)};

for i = 1:length(cRs)
    n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                           'flex_top_type', '4dof', ...
                                           'motion_sensor_type', 'struts', ...
                                           'actuator_type', 'flexible', ...
                                           'flex_bot_cRx', cRs(i), ...
                                           'flex_bot_cRy', cRs(i), ...
                                           'flex_top_cRx', cRs(i), ...
                                           'flex_top_cRy', cRs(i));

    G = exp(-s*Ts)*linearize(mdl, io, 0.0, options);
    G.InputName  = {'Va1', 'Va2', 'Va3', 'Va4', 'Va5', 'Va6'};
    G.OutputName = {'dL1', 'dL2', 'dL3', 'dL4', 'dL5', 'dL6'};

    Gs(i) = {G};
end

%% Tested bending dampings [Nm/(rad/s)]
xis = [1e-3, 5e-3, 1e-2, 5e-2, 1e-1];

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Encoders

%% Idenfity the transfer function from actuator to encoder for all bending dampins
Gs = {zeros(length(xis), 1)};

for i = 1:length(xis)
    n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                           'flex_top_type', '4dof', ...
                                           'motion_sensor_type', 'struts', ...
                                           'actuator_type', 'flexible', ...
                                           'actuator_xi', xis(i));

    G = exp(-s*Ts)*linearize(mdl, io, 0.0, options);
    G.InputName  = {'Va1', 'Va2', 'Va3', 'Va4', 'Va5', 'Va6'};
    G.OutputName = {'dL1', 'dL2', 'dL3', 'dL4', 'dL5', 'dL6'};

    Gs(i) = {G};
end

m_coef = 1e1;

%% Tested bending dampings [Nm/(rad/s)]
k_coefs = [1e-6, 5e-6, 1e-5, 5e-5, 1e-4];

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Encoders

%% Idenfity the transfer function from actuator to encoder for all bending dampins
Gs = {zeros(length(k_coefs), 1)};
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', 'flexible');

for i = 1:length(k_coefs)
    k_coef = k_coefs(i);

    G = exp(-s*Ts)*linearize(mdl, io, 0.0, options);
    G.InputName  = {'Va1', 'Va2', 'Va3', 'Va4', 'Va5', 'Va6'};
    G.OutputName = {'dL1', 'dL2', 'dL3', 'dL4', 'dL5', 'dL6'};

    Gs(i) = {G};
end

k_coef = 1e-6;

%% Tested bending dampings [Nm/(rad/s)]
m_coefs = [1e1, 5e1, 1e2, 5e2, 1e3];

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Encoders

%% Idenfity the transfer function from actuator to encoder for all bending dampins
Gs = {zeros(length(m_coefs), 1)};
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', 'flexible');

for i = 1:length(m_coefs)
    m_coef = m_coefs(i);

    G = exp(-s*Ts)*linearize(mdl, io, 0.0, options);
    G.InputName  = {'Va1', 'Va2', 'Va3', 'Va4', 'Va5', 'Va6'};
    G.OutputName = {'dL1', 'dL2', 'dL3', 'dL4', 'dL5', 'dL6'};

    Gs(i) = {G};
end

d_aligns = [[-0.05,  -0.3,  0];
            [ 0,      0.5,  0];
            [-0.1,   -0.3,  0];
            [ 0,      0.3,  0];
            [-0.05,   0.05, 0];
            [0,       0,    0]]*1e-3;

d_aligns = zeros(6,3);
% d_aligns(1,:) = [-0.05,  -0.3,  0]*1e-3;
d_aligns(2,:) = [ 0,      0.3,  0]*1e-3;

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', 'flexible', ...
                                       'actuator_d_align', d_aligns);

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Encoders

Gdvf = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');   io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dum'],  1, 'openoutput'); io_i = io_i + 1; % Force Sensors

Giff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');   io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Force Sensors

d_aligns = [[-0.05,  -0.3,  0];
            [ 0,      0.5,  0];
            [-0.1,   -0.3,  0];
            [ 0,      0.3,  0];
            [-0.05,   0.05, 0];
            [0,       0,    0]]*1e-3;

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', 'flexible', ...
                                       'actuator_d_align', d_aligns);

Gdvf_struts = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', 'flexible', ...
                                       'actuator_d_align', d_aligns);

Gdvf_plates = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', '2dof');

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');   io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dum'],  1, 'openoutput'); io_i = io_i + 1; % Force Sensors

Giff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% IFF Controller
Kiff_g1 = -(1/(s + 2*pi*40))*...    % LPF: provides integral action above 40Hz
           (s/(s + 2*pi*30))*...    % HPF: limit low frequency gain
           (1/(1 + s/2/pi/500))*... % LPF: more robust to high frequency resonances
           eye(6);                  % Diagonal 6x6 controller

%% IFF controller with Optimal gain
Kiff = 400*Kiff_g1;

save('mat/Kiff.mat', 'Kiff')

%% Initialize the Simscape model in closed loop
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', '2dof', ...
                                       'controller_type', 'iff');

%% Tested IFF gains
iff_gains = [4, 10, 20, 40, 100, 200, 400];

%% Identify the (damped) transfer function from u to dLm for different values of the IFF gain
Gd_iff = {zeros(1, length(iff_gains))};

clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Strut Displacement (encoder)

for i = 1:length(iff_gains)
    Kiff = iff_gains(i)*Kiff_g1*eye(6); % IFF Controller
    Gd_iff(i) = {exp(-s*Ts)*linearize(mdl, io, 0.0, options)};

    isstable(Gd_iff{i})
end

%% Estimated misalignement of the struts
d_aligns = [[-0.05,  -0.3,  0];
            [ 0,      0.5,  0];
            [-0.1,   -0.3,  0];
            [ 0,      0.3,  0];
            [-0.05,   0.05, 0];
            [0,       0,    0]]*1e-3;


%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', 'flexible', ...
                                       'actuator_d_align', d_aligns, ...
                                       'controller_type', 'iff');

%% Optimal IFF controller
load('Kiff.mat', 'Kiff');

%% Linearized inputs/outputs
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'],  1, 'openoutput'); io_i = io_i + 1; % Strut Displacement (encoder)

%% Identification of the plant
Gd_iff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

frf_ol  = load('Measurement_Z_axis.mat'); % Open-Loop
frf_iff = load('Measurement_Z_axis_damped.mat'); % IFF

%% Multiply by 10 (gain in m/s^2/V) and divide by 5 (number of accelerometers)
d_frf_ol = 10/5*(frf_ol.FFT1_H1_4_1_RMS_Y_Mod + frf_ol.FFT1_H1_7_1_RMS_Y_Mod + frf_ol.FFT1_H1_10_1_RMS_Y_Mod + frf_ol.FFT1_H1_13_1_RMS_Y_Mod + frf_ol.FFT1_H1_16_1_RMS_Y_Mod)./(2*pi*frf_ol.FFT1_H1_16_1_RMS_X_Val).^2;
d_frf_iff = 10/5*(frf_iff.FFT1_H1_4_1_RMS_Y_Mod + frf_iff.FFT1_H1_7_1_RMS_Y_Mod + frf_iff.FFT1_H1_10_1_RMS_Y_Mod + frf_iff.FFT1_H1_13_1_RMS_Y_Mod + frf_iff.FFT1_H1_16_1_RMS_Y_Mod)./(2*pi*frf_iff.FFT1_H1_16_1_RMS_X_Val).^2;

%% Load Identification Data
meas_data_lf = {};

for i = 1:6
    meas_data_lf(i) = {load(sprintf('mat/frf_exc_strut_%i_enc_plates_noise.mat', i), 't', 'Va', 'Vs', 'de')};
end

%% Coherence
coh_dvf = zeros(length(f), 6, 6);

for i = 1:6
    coh_dvf(:, :, i) = mscohere(meas_data_lf{i}.Va, meas_data_lf{i}.de, win, [], [], 1/Ts);
end

%% DVF Plant (transfer function from u to dLm)
G_dvf = zeros(length(f), 6, 6);

for i = 1:6
    G_dvf(:,:,i) = tfestimate(meas_data_lf{i}.Va, meas_data_lf{i}.de, win, [], [], 1/Ts);
end

%% Coherence for the IFF plant
coh_iff = zeros(length(f), 6, 6);

for i = 1:6
    coh_iff(:,:,i) = mscohere(meas_data_lf{i}.Va, meas_data_lf{i}.Vs, win, [], [], 1/Ts);
end

%% IFF Plant
G_iff = zeros(length(f), 6, 6);

for i = 1:6
    G_iff(:,:,i) = tfestimate(meas_data_lf{i}.Va, meas_data_lf{i}.Vs, win, [], [], 1/Ts);
end

save('mat/identified_plants_enc_plates.mat', 'f', 'Ts', 'G_iff', 'G_dvf')

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', 'flexible');

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');   io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/Fm'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensors

Giff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Encoders

Gdvf = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '2dof', ...
                                       'flex_top_type', '3dof', ...
                                       'motion_sensor_type', 'struts', ...
                                       'actuator_type', '2dof', ...
                                       'top_plate_type', 'rigid');

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Encoders

Gdvf = linearize(mdl, io, 0.0, options);

size(Gdvf)
isstable(Gdvf)

[sys,g] = balreal(Gdvf);  % Compute balanced realization
elim = (g<1e-4);         % Small entries of g are negligible states
rsys = modred(sys,elim); % Remove negligible states
size(rsys)

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');   io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/Fm'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensors

Giff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Initialize the Simscape model in closed loop
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', 'flexible');

%% Optimal IFF controller
load('Kiff.mat', 'Kiff')

%% Identify the (damped) transfer function from u to dLm for different values of the IFF gain
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Plate Displacement (encoder)

%% Transfer function from u to dL (open-loop)
Gd_ol = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

%% Initialize the Simscape model in closed loop
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', 'flexible', ...
                                       'controller_type', 'iff');

%% Transfer function from u to dL (IFF)
Gd_iff = exp(-s*Ts)*linearize(mdl, io, 0.0, options);

isstable(Gd_iff)

%% Load Identification Data
meas_iff_plates = {};

for i = 1:6
    meas_iff_plates(i) = {load(sprintf('mat/frf_exc_iff_strut_%i_enc_plates_noise.mat', i), 't', 'Va', 'Vs', 'de', 'u')};
end

% Sampling Time [s]
Ts = (meas_iff_plates{1}.t(end) - (meas_iff_plates{1}.t(1)))/(length(meas_iff_plates{1}.t)-1);

% Hannning Windows
win = hanning(ceil(1/Ts));

% And we get the frequency vector
[~, f] = tfestimate(meas_iff_plates{1}.Va, meas_iff_plates{1}.de, win, [], [], 1/Ts);

%% Estimation of the transfer function matrix from u to dL when IFF is applied
G_enc_iff_opt = zeros(length(f), 6, 6);

for i = 1:6
    G_enc_iff_opt(:,:,i) = tfestimate(meas_iff_plates{i}.Va, meas_iff_plates{i}.de, win, [], [], 1/Ts);
end

save('matlab/mat/damped_plant_enc_plates.mat', 'f', 'Ts', 'G_enc_iff_opt')

save('mat/damped_plant_enc_plates.mat', 'f', 'Ts', 'G_enc_iff_opt')

support.type = 1; % On top of vibration table

%% Initialize Nano-Hexapod
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '2dof', ...
                                       'flex_top_type', '3dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', '2dof');

payload.type = 1; % Payload / 1 "mass layer"

%% Identify the DVF Plant (transfer function from u to dLm)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Encoders

%% Identification for all the added payloads
G_dL = {};

for i = i_masses
    fprintf('i = %i\n', i)
    payload.type = i;
    G_dL(i+1) = {exp(-s*frf_ol.Ts)*linearize(mdl, io, 0.0, options)};
end

%% Identify the IFF Plant (transfer function from u to taum)
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/Fm'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensors

%% Identification for all the added payloads
G_tau = {};

for i = 0:3
    fprintf('i = %i\n', i)
    payload.type = i;
    G_tau(i+1) = {exp(-s*frf_ol.Ts)*linearize(mdl, io, 0.0, options)};
end

save('mat/sim_vib_table_m.mat', 'G_tau', 'G_dL')

%% IFF Controller
Kiff_g1 = (1/(s + 2*pi*20))*... % LPF: provides integral action above 20[Hz]
          (s/(s + 2*pi*20))*... % HPF: limit low frequency gain
          (1/(1 + s/2/pi/400)); % LPF: more robust to high frequency resonances

%% Optimal controller
g_opt = -2e2;
Kiff = g_opt*Kiff_g1*eye(6);

save('mat/Kiff_opt.mat', 'Kiff')

support.type = 1; % On top of vibration table

%% Initialize the Simscape model in closed loop
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '2dof', ...
                                       'flex_top_type', '3dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', '2dof', ...
                                       'controller_type', 'iff');

payload.type = 1; % Payload / 1 "mass layer"

%% Identify the (damped) transfer function from u to dLm
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Plate Displacement (encoder)

%% Identify for all add masses
G_dL = {};

for i = i_masses
    payload.type = i;
    G_dL(i+1) = {exp(-s*frf_ol.Ts)*linearize(mdl, io, 0.0, options)};
end

save('mat/sim_iff_vib_table_m.mat', 'G_dL');

sim_iff = load('sim_iff_vib_table_m.mat', 'G_dL');

i_masses = 0:3;

%% Load Identification Data
meas_added_mass = {};

for i_mass = i_masses
    for i_strut = 1:6
        meas_iff_mass(i_strut, i_mass+1) = {load(sprintf('frf_data_exc_strut_%i_iff_vib_table_%im.mat', i_strut, i_mass), 't', 'Va', 'de')};
    end
end

% Sampling Time [s]
Ts = (meas_iff_mass{1,1}.t(end) - (meas_iff_mass{1,1}.t(1)))/(length(meas_iff_mass{1,1}.t)-1);

% Hannning Windows
win = hanning(ceil(1/Ts));

% And we get the frequency vector
[~, f] = tfestimate(meas_iff_mass{1,1}.Va, meas_iff_mass{1,1}.de, win, [], [], 1/Ts);

%% DVF Plant (transfer function from u to dLm)
G_dL = {};

for i_mass = i_masses
    G_dL(i_mass+1) = {zeros(length(f), 6, 6)};
    for i_strut = 1:6
        G_dL{i_mass+1}(:,:,i_strut) = tfestimate(meas_iff_mass{i_strut, i_mass+1}.Va, meas_iff_mass{i_strut, i_mass+1}.de, win, [], [], 1/Ts);
    end
end

save('mat/frf_iff_vib_table_m.mat', 'f', 'Ts', 'G_dL');

%% Load Identification Data
meas_added_mass = {zeros(6,1)};

for i_strut = 1:6
    meas_iff_mass(i_strut) = {load(sprintf('frf_data_exc_strut_%i_iff_vib_table_1m_unbalanced.mat', i_strut), 't', 'Va', 'de')};
end

% Sampling Time [s]
Ts = (meas_iff_mass{1}.t(end) - (meas_iff_mass{1}.t(1)))/(length(meas_iff_mass{1}.t)-1);

% Hannning Windows
win = hanning(ceil(1/Ts));

% And we get the frequency vector
[~, f] = tfestimate(meas_iff_mass{1}.Va, meas_iff_mass{1}.de, win, [], [], 1/Ts);

%% DVF Plant (transfer function from u to dLm)
G_dL = zeros(length(f), 6, 6);
for i_strut = 1:6
    G_dL(:,:,i_strut) = tfestimate(meas_iff_mass{i_strut}.Va, meas_iff_mass{i_strut}.de, win, [], [], 1/Ts);
end

save('mat/frf_iff_unbalanced_vib_table_m.mat', 'f', 'Ts', 'G_dL');

%% Generate the Y-Z trajectory scan
Rx_yz = generateYZScanTrajectory(...
    'y_tot', 4e-6, ... % Length of Y scans [m]
    'z_tot', 4e-6, ... % Total Z distance [m]
    'n', 5, ...     % Number of Y scans
    'Ts', 1e-3, ... % Sampling Time [s]
    'ti', 1, ...    % Time to go to initial position [s]
    'tw', 0, ...    % Waiting time between each points [s]
    'ty', 0.6, ...  % Time for a scan in Y [s]
    'tz', 0.2);     % Time for a scan in Z [s]

%% Compute the reference in the frame of the legs
dL_ref = [J*Rx_yz(:, 2:7)']';

%% Generate the "tilt-spiral" trajectory scan
R_tilt = generateSpiralAngleTrajectory(...
    'R_tot',  20e-6, ... % Total Tilt [ad]
    'n_turn', 5, ...     % Number of scans
    'Ts',     1e-3, ...  % Sampling Time [s]
    't_turn', 1, ...     % Turn time [s]
    't_end',  1);        % End time to go back to zero [s]

%% List of points that draws "NASS"
ref_path = [ ...
    0, 0,0; % Initial Position
    0,0,1; 0,4,1; 3,0,1; 3,4,1; % N
    3,4,0; 4,0,0; % Transition
    4,0,1; 4,3,1; 5,4,1; 6,4,1; 7,3,1; 7,2,1; 4,2,1; 4,3,1; 5,4,1; 6,4,1; 7,3,1; 7,0,1; % A
    7,0,0; 8,0,0; % Transition
    8,0,1; 11,0,1; 11,2,1; 8,2,1; 8,4,1; 11,4,1; % S
    11,4,0; 12,0,0; % Transition
    12,0,1; 15,0,1; 15,2,1; 12,2,1; 12,4,1; 15,4,1; % S
    15,4,0;
           ];

%% Center the trajectory arround zero
ref_path = ref_path - (max(ref_path) - min(ref_path))/2;

%% Define the X-Y-Z cuboid dimensions containing the trajectory
X_max = 10e-6;
Y_max =  4e-6;
Z_max =  2e-6;

ref_path = ([X_max, Y_max, Z_max]./max(ref_path)).*ref_path; % [m]

%% Generating the trajectory
Rx_nass = generateXYZTrajectory('points', ref_path);

%% Lead to increase phase margin
a  = 2;  % Amount of phase lead / width of the phase lead / high frequency gain
wc = 2*pi*100; % Frequency with the maximum phase lead [rad/s]

H_lead = (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a)));

%% Low Pass filter to increase robustness
H_lpf = 1/(1 + s/2/pi/200);

%% Notch at the top-plate resonance
gm = 0.02;
xi = 0.3;
wn = 2*pi*700;

H_notch = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);

%% Decentralized HAC
Khac_iff_struts = -(1/(2.87e-5)) * ... % Gain
                  H_lead * ...       % Lead
                  H_notch * ...      % Notch
                  (2*pi*100/s) * ... % Integrator
                  eye(6);            % 6x6 Diagonal

save('mat/Khac_iff_struts.mat', 'Khac_iff_struts')

L_hac_iff_struts = pagemtimes(permute(frf_iff.G_dL{1}, [2 3 1]), squeeze(freqresp(Khac_iff_struts, frf_iff.f, 'Hz')));

%% Initialize the Simscape model in closed loop
n_hexapod = initializeNanoHexapodFinal('flex_bot_type', '4dof', ...
                                       'flex_top_type', '4dof', ...
                                       'motion_sensor_type', 'plates', ...
                                       'actuator_type', 'flexible', ...
                                       'controller_type', 'hac-iff-struts');

%% Identify the (damped) transfer function from u to dLm
clear io; io_i = 1;
io(io_i) = linio([mdl, '/du'], 1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1; % Plate Displacement (encoder)

%% Identification
Gd_iff_hac_opt = linearize(mdl, io, 0.0, options);

%% Verify the stability
isstable(Gd_iff_hac_opt)

%% Load the experimental data
load('hac_iff_struts_yz_scans.mat', 't', 'de')

%% Pose of the top platform from the encoder values
load('jacobian.mat', 'J');
Xe = [inv(J)*de']';

%% Generate the Y-Z trajectory scan
Rx_yz = generateYZScanTrajectory(...
    'y_tot', 4e-6, ... % Length of Y scans [m]
    'z_tot', 8e-6, ... % Total Z distance [m]
    'n', 5, ...     % Number of Y scans
    'Ts', 1e-3, ... % Sampling Time [s]
    'ti', 1, ...    % Time to go to initial position [s]
    'tw', 0, ...    % Waiting time between each points [s]
    'ty', 0.6, ...  % Time for a scan in Y [s]
    'tz', 0.2);     % Time for a scan in Z [s]

%% Lead
a  = 5;  % Amount of phase lead / width of the phase lead / high frequency gain
wc = 2*pi*110; % Frequency with the maximum phase lead [rad/s]

H_lead = (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a)));

%% Low Pass Filter
H_lpf = 1/(1 + s/2/pi/300);

%% Notch
gm = 0.02;
xi = 0.5;
wn = 2*pi*700;

H_notch = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);

%% HAC Controller
Khac_iff_struts = -2.2e4 * ... % Gain
            H_lead * ...       % Lead
            H_lpf * ...        % Lead
            H_notch * ...      % Notch
            (2*pi*100/s) * ... % Integrator
            eye(6);            % 6x6 Diagonal

%% Load Measurements
load('hac_iff_more_lead_nass_scan.mat', 't', 'de')

%% Compute the pose of the top platform
load('jacobian.mat', 'J');
Xe = [inv(J)*de']';

wc = 100; % Wanted crossover frequency [Hz]
[~, i_wc] = min(abs(frf_iff.f - wc)); % Indice corresponding to wc

%% Plant to be decoupled
frf_coupled = frf_iff.G_dL{2};
G_coupled = sim_iff.G_dL{2};

%% Compute RGA Matrix
RGA_wc = zeros(size(frf_coupled));
for i = 1:length(frf_iff.f)
    RGA_wc(i,:,:) = squeeze(G_dL_wc(i,:,:)).*inv(squeeze(G_dL_wc(i,:,:))).';
end

%% Compute RGA-number
RGA_wc_sum = zeros(size(RGA_wc, 1), 1);
for i = 1:size(RGA_wc, 1)
    RGA_wc_sum(i) = sum(sum(abs(eye(6) - squeeze(RGA_wc(i,:,:)))));
end

%% Compute the RGA matrix for the SVD decoupled plant
RGA_svd = zeros(size(frf_coupled));
for i = 1:length(frf_iff.f)
    RGA_svd(i,:,:) = squeeze(G_dL_svd(i,:,:)).*inv(squeeze(G_dL_svd(i,:,:))).';
end

%% Compute the RGA-number
RGA_svd_sum = zeros(size(RGA_svd, 1), 1);
for i = 1:length(frf_iff.f)
    RGA_svd_sum(i) = sum(sum(abs(eye(6) - squeeze(RGA_svd(i,:,:)))));
end

%% RGA Number for the SVD decoupled plant
figure;
plot(frf_iff.f, RGA_svd_sum, 'k-');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('RGA Number');
xlim([10, 1e3]); ylim([1e-2, 1e2]);

frf_new = frf_iff.G_dL{3};

function [ref] = generateXYZTrajectory(args)
% generateXYZTrajectory -
%
% Syntax: [ref] = generateXYZTrajectory(args)
%
% Inputs:
%    - args
%
% Outputs:
%    - ref - Reference Signal

arguments
    args.points double {mustBeNumeric} = zeros(2, 3) % [m]

    args.ti    (1,1) double {mustBeNumeric, mustBeNonnegative} = 1 % Time to go to first point and after last point [s]
    args.tw    (1,1) double {mustBeNumeric, mustBeNonnegative} = 0.5 % Time wait between each point [s]
    args.tm    (1,1) double {mustBeNumeric, mustBeNonnegative} = 1 % Motion time between points [s]

    args.Ts    (1,1) double {mustBeNumeric, mustBePositive} = 1e-3 % Sampling Time [s]
end

time_i = 0:args.Ts:args.ti;
time_w = 0:args.Ts:args.tw;
time_m = 0:args.Ts:args.tm;

% Go to initial position
xyz = (args.points(1,:))'*(time_i/args.ti);

% Wait
xyz = [xyz, xyz(:,end).*ones(size(time_w))];

% Scans
for i = 2:size(args.points, 1)
    % Go to next point
    xyz = [xyz, xyz(:,end) + (args.points(i,:)' - xyz(:,end))*(time_m/args.tm)];
    % Wait a litle bit
    xyz = [xyz, xyz(:,end).*ones(size(time_w))];
end

% End motion
xyz = [xyz, xyz(:,end) - xyz(:,end)*(time_i/args.ti)];

t = 0:args.Ts:args.Ts*(length(xyz) - 1);

ref = zeros(length(xyz), 7);

ref(:, 1) = t;
ref(:, 2:4) = xyz';

function [ref] = generateYZScanTrajectory(args)
% generateYZScanTrajectory -
%
% Syntax: [ref] = generateYZScanTrajectory(args)
%
% Inputs:
%    - args
%
% Outputs:
%    - ref - Reference Signal

arguments
    args.y_tot (1,1) double {mustBeNumeric, mustBePositive} = 10e-6 % [m]
    args.z_tot (1,1) double {mustBeNumeric, mustBePositive} = 10e-6 % [m]

    args.n     (1,1) double {mustBeInteger, mustBePositive} = 10 % [-]

    args.Ts    (1,1) double {mustBeNumeric, mustBePositive} = 1e-4 % [s]

    args.ti    (1,1) double {mustBeNumeric, mustBeNonnegative} = 1 % [s]
    args.tw    (1,1) double {mustBeNumeric, mustBeNonnegative} = 1 % [s]
    args.ty    (1,1) double {mustBeNumeric, mustBeNonnegative} = 1 % [s]
    args.tz    (1,1) double {mustBeNumeric, mustBeNonnegative} = 1 % [s]
end

time_i = 0:args.Ts:args.ti;
time_w = 0:args.Ts:args.tw;
time_y = 0:args.Ts:args.ty;
time_z = 0:args.Ts:args.tz;

% Go to initial position
y = (time_i/args.ti)*(args.y_tot/2);

% Wait
y = [y, y(end)*ones(size(time_w))];

% Scans
for i = 1:args.n
    if mod(i,2) == 0
        y = [y, -(args.y_tot/2) + (time_y/args.ty)*args.y_tot];
    else
        y = [y,  (args.y_tot/2) - (time_y/args.ty)*args.y_tot];
    end

    if i < args.n
        y = [y, y(end)*ones(size(time_z))];
    end
end

% Wait a litle bit
y = [y, y(end)*ones(size(time_w))];

% End motion
y = [y, y(end) - y(end)*time_i/args.ti];

% Go to initial position
z = (time_i/args.ti)*(args.z_tot/2);

% Wait
z = [z, z(end)*ones(size(time_w))];

% Scans
for i = 1:args.n
    z = [z, z(end)*ones(size(time_y))];

    if i < args.n
        z = [z, z(end) - (time_z/args.tz)*args.z_tot/(args.n-1)];
    end
end

% Wait a litle bit
z = [z, z(end)*ones(size(time_w))];

% End motion
z = [z, z(end) - z(end)*time_i/args.ti];

t = 0:args.Ts:args.Ts*(length(y) - 1);

ref = zeros(length(y), 7);

ref(:, 1) = t;
ref(:, 3) = y;
ref(:, 4) = z;

function [ref] = generateSpiralAngleTrajectory(args)
% generateSpiralAngleTrajectory -
%
% Syntax: [ref] = generateSpiralAngleTrajectory(args)
%
% Inputs:
%    - args
%
% Outputs:
%    - ref - Reference Signal

arguments
    args.R_tot  (1,1) double {mustBeNumeric, mustBePositive} = 10e-6 % [rad]
    args.n_turn (1,1) double {mustBeInteger, mustBePositive} = 5 % [-]
    args.Ts     (1,1) double {mustBeNumeric, mustBePositive} = 1e-3 % [s]
    args.t_turn (1,1) double {mustBeNumeric, mustBePositive} = 1 % [s]
    args.t_end  (1,1) double {mustBeNumeric, mustBePositive} = 1 % [s]
end

time_s = 0:args.Ts:args.n_turn*args.t_turn;
time_e = 0:args.Ts:args.t_end;

Rx = sin(2*pi*time_s/args.t_turn).*(args.R_tot*time_s/(args.n_turn*args.t_turn));
Ry = cos(2*pi*time_s/args.t_turn).*(args.R_tot*time_s/(args.n_turn*args.t_turn));

Rx = [Rx, 0*time_e];
Ry = [Ry, Ry(end) - Ry(end)*time_e/args.t_end];

t = 0:args.Ts:args.Ts*(length(Rx) - 1);

ref = zeros(length(Rx), 7);

ref(:, 1) = t;
ref(:, 5) = Rx;
ref(:, 6) = Ry;

function [M] = getTransformationMatrixAcc(Opm, Osm)
% getTransformationMatrixAcc -
%
% Syntax: [M] = getTransformationMatrixAcc(Opm, Osm)
%
% Inputs:
%    - Opm - Nx3 (N = number of accelerometer measurements) X,Y,Z position of accelerometers
%    - Opm - Nx3 (N = number of accelerometer measurements) Unit vectors representing the accelerometer orientation
%
% Outputs:
%    - M - Transformation Matrix

M = zeros(length(Opm), 6);

for i = 1:length(Opm)
    Ri = [0,         Opm(3,i), -Opm(2,i);
         -Opm(3,i),  0,         Opm(1,i);
          Opm(2,i), -Opm(1,i),  0];
    M(i, 1:3) = Osm(:,i)';
    M(i, 4:6) = Osm(:,i)'*Ri;
end

end

function [J] = getJacobianNanoHexapod(Hbm)
% getJacobianNanoHexapod -
%
% Syntax: [J] = getJacobianNanoHexapod(Hbm)
%
% Inputs:
%    - Hbm - Height of {B} w.r.t. {M} [m]
%
% Outputs:
%    - J - Jacobian Matrix

Fa = [[-86.05,  -74.78, 22.49],
      [ 86.05,  -74.78, 22.49],
      [ 107.79, -37.13, 22.49],
      [ 21.74,  111.91, 22.49],
      [-21.74,  111.91, 22.49],
      [-107.79, -37.13, 22.49]]'*1e-3; % Ai w.r.t. {F} [m]

Mb = [[-28.47, -106.25, -22.50],
      [ 28.47, -106.25, -22.50],
      [ 106.25, 28.47,  -22.50],
      [ 77.78,  77.78,  -22.50],
      [-77.78,  77.78,  -22.50],
      [-106.25, 28.47,  -22.50]]'*1e-3; % Bi w.r.t. {M} [m]

H = 95e-3; % Stewart platform height [m]
Fb = Mb + [0; 0; H]; % Bi w.r.t. {F} [m]

si = Fb - Fa;
si = si./vecnorm(si); % Normalize

Bb = Mb - [0; 0; Hbm];

J = [si', cross(Bb, si)'];