From 02c6d1a9c8beb81d303ed4d8b55f3d4a85961ada Mon Sep 17 00:00:00 2001 From: Lucas Antonio Pelike Date: Wed, 15 May 2024 16:38:38 -0300 Subject: [PATCH] Add scripts related to ARX estimation of generalized plant Generalized plant is obtained from sysid experiments considering correctors to BPMs path. Co-authored-by: Guilherme Ricioli --- sysid/generate_plant_models.m | 45 ++++++++++++++++++++++ sysid/plant_arx_fit.m | 72 +++++++++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+) create mode 100644 sysid/generate_plant_models.m create mode 100644 sysid/plant_arx_fit.m diff --git a/sysid/generate_plant_models.m b/sysid/generate_plant_models.m new file mode 100644 index 0000000..6c834b6 --- /dev/null +++ b/sysid/generate_plant_models.m @@ -0,0 +1,45 @@ +% Copyright (C) 2024 CNPEM (cnpem.br) +% Author: Lucas Pelike +% Modified by: Guilherme Ricioli + +clc; clear; + +prbs_ol_acq_fpath = '/current/'; +ps_names_fpath = 'ps_names.txt'; + +ps_names_table = readtable(ps_names_fpath, Delimiter=',', ... + ReadVariableNames=false, TextType='string'); +ncorr = size(ps_names_table, 1); + +excluded_corr = [1, 80, 81, 160]; + +sys = cell(ncorr, 1); +fit = zeros(ncorr, 1); + +figure(); +for i=1:ncorr + tic; + ps_name = string(ps_names_table{i, 1}); + fprintf('Corrector %d: %s\n', i, ps_name); + if ismember(i, excluded_corr) + sys{i} = NaN; + else + fpath = strcat(prbs_ol_acq_fpath, ps_name, '.mat'); + sys{i} = plant_arx_fit(fpath, [6 6 2]); + fit(i) = sys{i}.Report.Fit.FitPercent; + + opts = bodeoptions; + opts.FreqUnits = 'Hz'; + bode(sys{i}, opts); + hold on; + end + fprintf('Elapsed time: %f s\n', toc); +end + +save('sysid_res', 'sys'); + +figure(); +plot(fit); +title('Fit Percent'); +xlabel('Corrector index'); +ylabel('Fit [%]'); diff --git a/sysid/plant_arx_fit.m b/sysid/plant_arx_fit.m new file mode 100644 index 0000000..1669a3a --- /dev/null +++ b/sysid/plant_arx_fit.m @@ -0,0 +1,72 @@ +% Copyright (C) 2024 CNPEM (cnpem.br) +% Author: Lucas Pelike +% Modified by: Guilherme Ricioli + +function sys = plant_arx_fit(fpath, arx_params, n_prbs_T_to_use) +% plant_arx_fit +% +% Provides an ARX fit to the studied system with PRBS excitation +% +% sys = plant_arx_fit(fpath, arx_params, n_prbs_T_to_use) +% +% INPUTS: +% fpath: Filepath to the plant PRBS acquistion MATLAB object +% generated by prbs_ol_acq_h5_to_mat.m script +% arx_params: ARX estimation parameters in the format [na nb nk] (see +% 'arx' documentation) +% n_prbs_T_to_use: (optional parameter) Number of PRBS periods to be +% considered +% +% OUTPUTS: +% sys: Resulting system in idpoly format + +% Loads .mat fpath +prbs_ol_acq = load(fpath); +Ts = 1/prbs_ol_acq.data.sampling_frequency; +prbs_u = prbs_ol_acq.data.prbs_data; +prbs_lfsr_len = prbs_ol_acq.data.prbs_lfsr_len(1); +prbs_step_duration = prbs_ol_acq.data.prbs_step_duration(1); +prbs_mov_avg_taps = double(prbs_ol_acq.data.prbs_mov_avg_taps); +bpm_y = prbs_ol_acq.data.orb(prbs_ol_acq.data.bpm_idx_max_response, :); + +% PRBS period length +prbs_T = (2^(prbs_lfsr_len) - 1)*prbs_step_duration; + +% Removes transient +n_prbs_T_transient = 4; +assert(length(prbs_u) > prbs_T*n_prbs_T_transient, ... + "Could not remove transient: dataset is too small"); +prbs_u = prbs_u(prbs_T*n_prbs_T_transient + 1:end); +bpm_y = bpm_y(prbs_T*n_prbs_T_transient + 1:end); + +% Guarantees that there's an integer amount of PRBS periods in datasets +rem = mod(length(prbs_u), prbs_T); +prbs_u = prbs_u(1:end - rem); +bpm_y = bpm_y(1:end - rem); + +% Reduces dataset (optional parameter) +if exist('n_prbs_T_to_use', 'var') + n_p = length(prbs_u)/prbs_T; + assert(n_prbs_T_to_use <= n_p, ... + "Could not reduce dataset: not enough PRBS periods"); + prbs_u = prbs_u(1:n_prbs_T_to_use*prbs_T); + bpm_y = bpm_y(1:n_prbs_T_to_use*prbs_T); +end + +% Mimics the moving average filter that's applied to the PRBS excitation in +% gateware so we don't end up identifying it +mov_avg_filt = dfilt.dffir(ones(1, 2^prbs_mov_avg_taps)/2^prbs_mov_avg_taps); +prbs_u = filter(mov_avg_filt, prbs_u); + +% Averaging +prbs_u_avg = mean(reshape(prbs_u, prbs_T, []), 2); +bpm_y_avg = mean(reshape(bpm_y, prbs_T, []), 2); + +% Removes the mean value +prbs_u_avg = prbs_u_avg - mean(prbs_u_avg); +bpm_y_avg = bpm_y_avg - mean(bpm_y_avg); + +% ARX fit +plant_iddata = iddata(bpm_y_avg, prbs_u_avg, Ts); +sys = arx(plant_iddata, arx_params, ... + arxOptions('Focus', 'Simulation', 'EnforceStability', true));