Adding load mitigating controls and OWC

.github/workflows/run-tests-dev.yml

@@ -31,13 +31,15 @@ jobs:
                  End_Stops,
                  Free_Decay,
                  Generalized_Body_Modes,
+                 Load_Mitigating_Controls,
                  Mean_Drift,
                  Morison_Element,
                  MOST,
                  Multiple_Condition_Runs,
                  Multiple_Wave_Headings,
                  Nonhydro_Body,
                  Nonlinear_Hydro,
+                 OWC,
                  Passive_Yaw,
                  PTO-Sim,
                  Radiation_Force_Options,
@@ -50,8 +52,12 @@ jobs:
             products: Simulink Simscape Simscape_Multibody Simscape_Fluids
           - folder: Controls
             products: Simulink Simscape Simscape_Multibody Control_System_Toolbox Optimization_Toolbox System_Identification_Toolbox Statistics_and_Machine_Learning_Toolbox Symbolic_Math_Toolbox
+          - folder: Load_Mitigating_Controls
+            products: Simulink Simscape Simscape_Multibody DSP_System_Toolbox
+          - folder: OWC
+            products: Simulink Simscape Simscape_Multibody Control_System_Toolbox
           - folder: WECCCOMP
-            products: Simulink Simscape Simscape_Multibody Control_System_Toolbox Optimization_Toolbox System_Identification_Toolbox Statistics_and_Machine_Learning_Toolbox 
+            products: Simulink Simscape Simscape_Multibody Control_System_Toolbox Optimization_Toolbox System_Identification_Toolbox Statistics_and_Machine_Learning_Toolbox
     name: ${{ matrix.folder }} MATLAB ${{ matrix.release }} on Linux
     steps:
       - name: Check out repository (repository dispatch)

.gitignore

@@ -88,6 +88,7 @@ yaw*.mat
 !/Generalized_Body_Modes/hydroData/WAMIT/*
 !/Mean_Drift/WAMIT/inputFiles/*
 !/Mean_Drift/WAMIT/outputFiles/*
+!/OWC/OrificeModel/hydroData/*.h5
 
 # Other
 .DS_Store

Load_Mitigating_Controls/CalcImpedance/Identify_Impedance.m

@@ -0,0 +1,158 @@
+% calculate impedance force/velocity from WEC-Sim output
+close all;
+baseDir = pwd;
+
+%% describe signal properties
+
+% multisine has 300 s period: grabbing a middle 300 s to avoid startup
+% transients
+t_start = 200; 
+t_out = 500; 
+load('mcrOut'); % loads wecsim outputs from each tested time series
+
+% define frequency range of interest, should correspond to excited band in
+% multisine
+f_min = 0.015; % min frequency in Hz
+f_max = 1.5;
+T_rep = 300;
+
+% find indices of start/stop times in outputs
+[~,idx_s] = min(abs(mcr.OutputLog1.time - t_start));
+[~,idx_e] = min(abs(mcr.OutputLog1.time - t_out));
+
+% define a time vector and calculate time step of data based on indices
+tvec = mcr.OutputLog1.time(idx_s:idx_e-1);
+dt = mean(diff(tvec));
+
+%% build data matrices
+
+% for each multisine run, append the force and velocity data (time domain)
+for k = 1:6
+    nameField = strcat(['OutputLog' num2str(k)]);
+    force(1,k,:) = mcr.(nameField).forces(t_start/dt : t_out/dt -1,1); % surge
+    force(2,k,:) = mcr.(nameField).forces(t_start/dt : t_out/dt -1,2); % heave
+    force(3,k,:) = mcr.(nameField).forces(t_start/dt : t_out/dt -1,3); % pitch
+    vel(1,k,:) = mcr.(nameField).vel(t_start/dt : t_out/dt -1,1); % surge 
+    vel(2,k,:) = mcr.(nameField).vel(t_start/dt : t_out/dt -1,2); % heave
+    vel(3,k,:) = mcr.(nameField).vel(t_start/dt : t_out/dt -1,3); % pitch
+end
+
+%% Move to frequency domain
+
+% take fft along 3rd dimension 
+[freq, FORCE] = ampSpectraForZcalc(force,tvec,T_rep);
+[~, VEL]= ampSpectraForZcalc(vel,tvec,T_rep);
+[~, ind_f_min] = min(abs(f_min-freq));
+[~, ind_f_max] = min(abs(f_max-freq));
+
+% restrict to frequency range of interest
+f_len = (ind_f_max-ind_f_min)+1;
+f_vec = freq(ind_f_min:ind_f_max);
+% restrict to of-interest frequency range
+FORCE=FORCE(:,:,[ind_f_min:ind_f_max]);
+VEL=VEL(:,:,[ind_f_min:ind_f_max]);
+
+%% Calculate impedance and admittance
+
+Z = zeros(3,3,f_len); % initialize impdedance
+Y = zeros(3,3,f_len); % initialize admittance
+
+for k = 1: f_len % solve impedance, frequency-by-frequency
+    Y(:,:,k) = VEL(:,:,k)/FORCE(:,:,k);
+    Z(:,:,k) = FORCE(:,:,k)/VEL(:,:,k);
+end
+
+cd(baseDir)
+
+% depending on parameters, there may be a negative real part of impedance.
+% This is not physically meaningful, and is a result of some numerical
+% errors. 
+
+% compare to BEM form from *h5 data
+iw = permute(1i * body.hydroData.simulation_parameters.w,[3,1,2]); 
+static_mass = diag([body.mass(1); body.mass(1); body.inertia(1)]);
+added_mass = simu.rho .* body.hydroData.hydro_coeffs.added_mass.all([1:2:5],[1:2:5],:);
+rad_damp = simu.rho .*repmat(permute(body.hydroData.simulation_parameters.w,[3,1,2]),3,3,1) .* body.hydroData.hydro_coeffs.radiation_damping.all([1:2:5],[1:2:5],:);
+hyd_stff = repmat(diag([24000; simu.rho .* simu.gravity .* body.hydroData.hydro_coeffs.linear_restoring_stiffness(3,3); simu.rho .* simu.gravity .* body.hydroData.hydro_coeffs.linear_restoring_stiffness(5,5)])...
+    ,1,1,length(iw)); % added spring in surge
+Z_bem = iw .* (static_mass + added_mass) + rad_damp + hyd_stff./iw;  
+
+%% create impedance plot comparing BEM-calc'd to WEC-Sim calc'd
+
+% bode plot, diagonals
+figure; clf;
+subplot(2,1,1)
+semilogx(2*pi*f_vec,mag2db(squeeze(abs(Z(1,1,:)))),'r','LineWidth',1.4); % surge
+hold on; grid on;
+semilogx(2*pi*f_vec,mag2db(squeeze(abs(Z(2,2,:)))),'b','LineWidth',1.4); % heave
+semilogx(2*pi*f_vec,mag2db(squeeze(abs(Z(3,3,:)))),'g','LineWidth',1.4); % pitch
+semilogx(body.hydroData.simulation_parameters.w,mag2db(squeeze(abs(Z_bem(1,1,:)))),'--r','LineWidth',1.4); % surge BEM
+semilogx(body.hydroData.simulation_parameters.w,mag2db(squeeze(abs(Z_bem(2,2,:)))),'--b','LineWidth',1.4); % heave BEM
+semilogx(body.hydroData.simulation_parameters.w,mag2db(squeeze(abs(Z_bem(3,3,:)))),'--g','LineWidth',1.4); % pitch BEM
+ylabel('|Z| db')
+xlim([0.5 10])
+s(1) = semilogx([0.5 10],[NaN NaN],'-k','LineWidth',1.2);
+s(2) = semilogx([0.5 10],[NaN NaN],'--k','LineWidth',1.2);
+s(3) = semilogx([0.5 10],[NaN NaN],'ob','MarkerFaceColor','b');
+s(4) = semilogx([0.5 10],[NaN NaN],'or','MarkerFaceColor','r');
+s(5) = semilogx([0.5 10],[NaN NaN],'og','MarkerFaceColor','g');
+legend(s,{'Z_i','Z_{BEM}', 'Heave', 'Surge', 'Pitch'},'AutoUpdate','off','Location','SouthOutside','Orientation','Horizontal');
+subplot(2,1,2)
+semilogx(2*pi*f_vec,(180/pi).*atan2(squeeze(imag(Z(1,1,:))),squeeze(real(Z(1,1,:)))),'r','LineWidth',1.4);
+hold on; grid on;
+semilogx(2*pi*f_vec,(180/pi).*atan2(squeeze(imag(Z(2,2,:))),squeeze(real(Z(2,2,:)))),'b','LineWidth',1.4);
+semilogx(2*pi*f_vec,(180/pi).*atan2(squeeze(imag(Z(3,3,:))),squeeze(real(Z(3,3,:)))),'g','LineWidth',1.4);
+semilogx(body.hydroData.simulation_parameters.w,(180/pi).*atan2(squeeze(imag(Z_bem(1,1,:))),squeeze(real(Z_bem(1,1,:)))),'--r','LineWidth',1.4);
+semilogx(body.hydroData.simulation_parameters.w,(180/pi).*atan2(squeeze(imag(Z_bem(2,2,:))),squeeze(real(Z_bem(2,2,:)))),'--b','LineWidth',1.4);
+semilogx(body.hydroData.simulation_parameters.w,(180/pi).*atan2(squeeze(imag(Z_bem(3,3,:))),squeeze(real(Z_bem(3,3,:)))),'--g','LineWidth',1.4);
+xlim([0.5 10])
+xlabel('Freq. (rad/s)')
+ylabel('Phase(^o)')
+
+% bode plot, off-diagonals
+figure; clf;
+subplot(2,1,1)
+semilogx(2*pi*f_vec,mag2db(squeeze(abs(Z(1,3,:)))),'m','LineWidth',1.4); % pitch-to-surge coupling
+hold on; grid on
+semilogx(2*pi*f_vec,mag2db(squeeze(abs(Z(3,1,:)))),'c','LineWidth',1.4); % surge-to-pitch coupling
+semilogx(body.hydroData.simulation_parameters.w,mag2db(squeeze(abs(Z_bem(1,3,:)))),'--m','LineWidth',1.4); % pitch-to-surge coupling, BEM
+semilogx(body.hydroData.simulation_parameters.w,mag2db(squeeze(abs(Z_bem(3,1,:)))),'--c','LineWidth',1.4); % surge-to-pitch coupling, BEM
+ylabel('|Z| db')
+xlim([0.5 10])
+s2(1) = semilogx([0.5 10],[NaN NaN],'-k','LineWidth',1.2);
+s2(2) = semilogx([0.5 10],[NaN NaN],'--k','LineWidth',1.2);
+s2(3) = semilogx([0.5 10],[NaN NaN],'om','MarkerFaceColor','m');
+s2(4) = semilogx([0.5 10],[NaN NaN],'oc','MarkerFaceColor','c');
+legend(s2,{'Z_i','Z_{BEM}', 'Pitch-to-surge coupling', 'Surge-to-pitch coupling'},'AutoUpdate','off','Location','SouthOutside','Orientation','Horizontal');
+subplot(2,1,2)
+semilogx(2*pi*f_vec,(180/pi).*atan2(squeeze(imag(Z(1,3,:))),squeeze(real(Z(1,3,:)))),'m','LineWidth',1.4);
+hold on; grid on;
+semilogx(2*pi*f_vec,(180/pi).*atan2(squeeze(imag(Z(3,1,:))),squeeze(real(Z(3,1,:)))),'c','LineWidth',1.4);
+semilogx(body.hydroData.simulation_parameters.w,(180/pi).*atan2(squeeze(imag(Z_bem(1,3,:))),squeeze(real(Z_bem(1,3,:)))),'--m','LineWidth',1.4);
+semilogx(body.hydroData.simulation_parameters.w,(180/pi).*atan2(squeeze(imag(Z_bem(3,1,:))),squeeze(real(Z_bem(3,1,:)))),'--c','LineWidth',1.4);
+xlim([0.5 10])
+xlabel('Freq. (rad/s)')
+ylabel('Phase(^o)')
+
+figure; clf;
+subplot(2,1,1)
+semilogx(2*pi*f_vec,real(squeeze(Z(1,1,:))),'b','LineWidth',1.4); % full model
+hold on; grid on;
+semilogx(2*pi*f_vec,real(squeeze(Z(2,2,:))),'r','LineWidth',1.4); % 3DOF model
+semilogx(2*pi*f_vec,real(squeeze(Z(3,3,:))),'g','LineWidth',1.4); % 1DOF model
+%semilogx(body.hydroData.simulation_parameters.w,real(squeeze(Z_bem(1,1,:))),'--b','LineWidth',1.4);
+%semilogx(body.hydroData.simulation_parameters.w,real(squeeze(Z_bem(2,2,:))),'--r','LineWidth',1.4);
+%semilogx(body.hydroData.simulation_parameters.w,real(squeeze(Z_bem(3,3,:))),'--g','LineWidth',1.4);
+ylabel('Re')
+xlim([0.5 10])
+subplot(2,1,2)
+semilogx(2*pi*f_vec,imag(squeeze(Z(1,1,:))),'b','LineWidth',1.4);
+hold on; grid on;
+semilogx(2*pi*f_vec,imag(squeeze(Z(2,2,:))),'r','LineWidth',1.4);
+semilogx(2*pi*f_vec,imag(squeeze(Z(3,3,:))),'g','LineWidth',1.4);
+%semilogx(body.hydroData.simulation_parameters.w,imag(squeeze(Z_bem(1,1,:))),'--b','LineWidth',1.4);
+%semilogx(body.hydroData.simulation_parameters.w,imag(squeeze(Z_bem(2,2,:))),'--r','LineWidth',1.4);
+%semilogx(body.hydroData.simulation_parameters.w,imag(squeeze(Z_bem(3,3,:))),'--g','LineWidth',1.4);
+xlim([0.5 10])
+xlabel('Freq. (rad/s)')
+ylabel('Im')

Load_Mitigating_Controls/CalcImpedance/ampSpectraForZcalc.m

@@ -0,0 +1,29 @@
+function [freq,compSpec]=ampSpectraForZcalc(Data,time,rep_len)
+
+%% INPUTS:
+% Data: array of 3D matrices. Index 1 is the DOF, index 2 is the test number,
+% and index 3 is the time series (along which to take the fft). Ensure it
+% is periodic or you'll have spectral leakage, zero padded to reduce
+% time: the time series (in s) related to the Data
+%% OUTPUTS:
+% freq: the vector of frequencies for the one-sided amplitude spectrum
+% ampSpec: the one-sided amplitude spectrum
+
+% detrend, zero pad, and fft
+Data = Data - mean(Data,3); % detrend
+L=length(Data(1,1,:)); 
+%Ln=2^nextpow2(L);
+compSpec=fft(Data,[],3)./L;
+compSpec=compSpec(:,:,(1:L/2+1));
+
+% handles frequency vector
+dt=mean(diff(time));
+if var(diff(time)) > 1e-5;
+    warning('Non-uniform sample time. Must resample data and try again')
+end
+Fs=1/dt;
+df = 1/rep_len; 
+freq=0:df:Fs/2;
+
+end
+

Load_Mitigating_Controls/CalcImpedance/userDefinedFunctions.m

@@ -0,0 +1,21 @@
+% userDefinedFunction for impedance calculation, saves the necessary
+% structure. 
+if exist('imcr','var') 
+    nameField = strcat(['OutputLog' num2str(imcr)]);
+
+    % save velocities
+    mcr.(nameField).vel =[output.bodies(1).velocity(:,1), output.bodies(1).velocity(:,3),output.bodies(1).velocity(:,5)];
+
+    % save forces 
+    mcr.(nameField).forces =[F_actuation.Data(:,1),F_actuation.Data(:,2),-1.*F_actuation.Data(:,3)];
+
+    % save time
+    mcr.(nameField).time(:,1) = F_actuation.Time;
+
+    if imcr==6; % save fully-populated structure
+        save('mcrOut','mcr')
+        Identify_Impedance;
+    end
+end
+
+

Load_Mitigating_Controls/CalcImpedance/userDefinedFunctionsMCR.m

@@ -0,0 +1,9 @@
+% This is an intentionally blank script. Chances are, you have another
+% userDefinedFunctionsMCR somewhere in your path and that can potentially
+% error out or cause undersirable behavior for this example if it tries to 
+% run. 
+
+% To that end, this blank function was included to do nothing and avoid
+% path issues.
+
+return

Load_Mitigating_Controls/CalcImpedance/wecSimInputFile.m

@@ -0,0 +1,59 @@
+%% Simulation Data
+simu = simulationClass();                 % Initialize Simulation Class
+simu.simMechanicsFile = 'WaveBot3Dof_ID'; % Specify Simulink Model File
+simu.mode = 'normal';                     % Specify Simulation Mode ('normal','accelerator','rapid-accelerator')
+simu.explorer = 'on';                     % Turn SimMechanics Explorer (on/off)
+simu.startTime = 0;                       % Simulation Start Time [s]
+simu.rampTime = 20;                       % Wave Ramp Time [s]
+simu.endTime = 600;                       % Simulation End Time [s]
+simu.solver = 'ode45';                    % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step 
+simu.dt = 0.01; 						  % Simulation time-step [s]
+simu.rho = 1025;                          % water density
+simu.domainSize = 100;                    % domain size for visualization
+simu.cicEndTime = 10;
+simu.mcrMatFile = 'mcrImpedanceRuns.mat'; 
+
+%% load multisine
+% This is handled by the 'LoadFile' header in mcrCaseFile, so keep
+% commented out for wecSimMCR runs. For single wecSim runs, uncomment:
+load('multisine3DOFA'); 
+
+%% Wave Information 
+% noWaveCIC, no waves with radiation CIC  
+waves = waveClass('noWaveCIC');       % Initialize Wave Class and Specify Type  
+
+%% System Identification %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% usually only use one at a time.
+FBEnable = 0; % one to enable feedback controller
+FFEnable = 1; % one to enable feedforward actuation
+
+%% define forcing multisines amplitude, by degree of freedom (surge heave pitch)
+OpenLoopGains = [50;100;4]; 
+
+%% Body Data
+% Float
+body(1) = bodyClass('../hydroData/waveBotBuoy.h5');      
+body(1).geometryFile = '../geometry/float.stl' ;  % Location of Geomtry File
+body(1).mass =1156.5; % heave mass is 893, surge mass is 1420, averaged used 
+body(1).inertia = [84 84 84];  %Moment of Inertia [kg*m^2], roll and yaw will be unused
+body(1).linearDamping(1:2:5)=[1000 1000 100]; 
+body(1).quadDrag.cd = [1.15 1.15 1 0.5 0.5 0]; %
+body(1).quadDrag.area = [2.9568 2.9568 5.4739 5.4739 5.4739 0]; %
+
+%% PTO and Constraint Parameters
+% Floating (3DOF) Joint
+constraint(1) = constraintClass('Constraint1'); % Initialize Constraint Class for Constraint1
+constraint(1).location = [0 0 0];                    % Constraint Location [m]
+
+%% Mooring Parameters
+mooring(1) = mooringClass('mooring');
+mooring(1).location = [0,0,0];
+mooring(1).matrix.stiffness = zeros(6,6);
+mooring(1).matrix.stiffness(1,1) = 24000;
+mooring(1).matrix.stiffness(2,2) = 24000;
+mooring(1).matrix.stiffness(3,3) = 5000;
+mooring(1).matrix.stiffness(5,5) = 1000;
+mooring(1).matrix.damping = zeros(6,6);
+mooring(1).matrix.damping(1,1) = 5;
+mooring(1).matrix.damping(3,3) = 5;
+mooring(1).matrix.preTension = zeros(1,6); 

Load_Mitigating_Controls/ControlTests/interp_impedance.m

@@ -0,0 +1,31 @@
+function Z_interp = interp_impedance(Z,f,f_interp);
+
+% INPUTS
+% Z: the impedance FRF to interpolate, with the frequency content on the
+%   third dimensions (i.e., size is [DOF,DOF,nFreq].
+% f: the frequency associated with the third dimension of the input
+%   impedance model Z.
+% f_interp: the frequency to which you should interpolate. f should bound
+%   f_vec and be in common units.
+% OUTPUTS
+% Z_interp: the interpolated Z 
+
+dims = size(Z(:,:,1));
+
+rPart = zeros(dims(1),dims(2),length(f_interp));
+iPart = zeros(dims(1),dims(2),length(f_interp));
+
+for k = 1:dims(1)
+    for kk = 1:dims(2) 
+        rPart(k,kk,:) = interp1(f, squeeze(real(Z(k,kk,:))),f_interp); 
+        iPart(k,kk,:) = interp1(f, squeeze(imag(Z(k,kk,:))),f_interp);
+    end
+end
+
+Z_interp = rPart + 1i * iPart; 
+
+end
+
+
+
+

Load_Mitigating_Controls/ControlTests/userDefinedFunctions.m

@@ -0,0 +1,23 @@
+% userDefinedFunction for impedance calculation, saves the necessary
+% structure.
+if exist('imcr','var');
+    nameField = strcat(['OutputLog' num2str(imcr)]);
+
+    % save velocities
+    mcr.(nameField).F_PTO =F_PTO.data;
+
+    % save forces
+    mcr.(nameField).F_TOT =F_TOT.data;
+
+    % save gains
+    mcr.(nameField).Gains = gains.data;
+
+    % save time
+    mcr.(nameField).time(:,1) = F_PTO.time;
+
+    if imcr==51 % save fully-populated structure
+        save('mcrOut','mcr')
+    end
+end
+
+

Load_Mitigating_Controls/ControlTests/userDefinedFunctionsMCR.m

@@ -0,0 +1,9 @@
+% This is an intentionally blank script. Chances are, you have another
+% userDefinedFunctionsMCR somewhere in your path and that can potentially
+% error out or cause undersirable behavior for this example if it tries to 
+% run. 
+
+% To that end, this blank function was included to do nothing and avoid
+% path issues.
+
+return