Merge pull request #30 from jtgrasb/Controls

Controls Applications Cases

Cable/userDefinedFunctions.m

@@ -30,7 +30,7 @@
 figure()
 subplot(2,1,1)
 initialLength = norm(cable(1).base.location  - cable(1).follower.location );
-plot([0 output.cables(1).time(end)], [cable(1).length cable(1).length],'k--',...
+plot([0 output.cables(1).time(end)], [cable(1).cableLength cable(1).cableLength],'k--',...
     output.cables(1).time, output.cables(1).position(:,3) + initialLength) 
 ylabel('Displacement [m]');
 title('Cable Displacement');

Cable/wecSimInputFile.m

@@ -75,7 +75,7 @@
 cable(1) = cableClass('Cable','constraint(2)','constraint(3)');
 cable(1).stiffness = 1000000;
 cable(1).damping = 100;
-cable(1).length = 17.8; % Cable equilibrium length [m] 
+cable(1).cableLength = 17.8; % Cable equilibrium length [m] 
 % cable(1).preTension = 5100000; % Cable equilibrium pre-tension [N]
 cable(1).quadDrag.cd = [1.4 1.4 1.4 0 0 0];
 cable(1).quadDrag.area = [10 10 10 0 0 0];

Controls/Declutching/declutching.m

@@ -0,0 +1,15 @@
+function [Force, tDeclutch, tNormal, vel] = declutching(vel, tDeclutch, tNormal, velPrev, declutchTime, dt, Kp)
+
+% Declutch when velocity is zero and reengage when declutching time is
+% reached
+if ((diff(sign([vel,velPrev]))~=0 && tDeclutch==0 && tNormal>.2) || (tDeclutch>0 && tDeclutch<declutchTime))
+    Force = 0;
+    tDeclutch = tDeclutch + dt;
+    tNormal = 0;
+else
+    Force = -Kp*vel;
+    tNormal = tNormal + dt;
+    tDeclutch = 0;
+end
+
+end

Controls/Declutching/mcrBuildTimes.m

@@ -0,0 +1,11 @@
+close all;clear all;clc;
+%%
+mcr = struct();
+mcr.header = ["controller(1).declutching.declutchTime"];
+mcr.cases = zeros(9,1);
+%%
+time = linspace(0,1.2,13);
+%%
+mcr.cases = time';
+%%
+save mcrCases mcr

Controls/Declutching/optimalTimeCalc.m

@@ -0,0 +1,73 @@
+% This script identifies the dynamics of the float in the respective wave 
+% conditions and determines the optimal proportional gain value for a 
+% passive controller (for regular waves)
+
+close all; clear all; clc;
+
+% Inputs (from wecSimInputFile)
+simu = simulationClass();
+body(1) = bodyClass('../hydroData/sphere.h5');
+waves.height = 1;
+waves.period = 3.5;
+
+% Load hydrodynamic data for float from BEM
+hydro = readBEMIOH5(body.h5File, 1, body.meanDrift);
+
+% Define wave conditions
+H = waves.height;
+A = H/2;
+T = waves.period;
+omega = (1/T)*(2*pi);
+
+% Define excitation force based on wave conditions
+ampSpect = zeros(length(hydro.simulation_parameters.w),1);
+[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w));
+ampSpect(closestIndOmega) = A;
+FeRao = squeeze(hydro.hydro_coeffs.excitation.re(3,:))'*simu.rho*simu.gravity + squeeze(hydro.hydro_coeffs.excitation.im(3,:))'*simu.rho*simu.gravity*1j;
+Fexc = ampSpect.*FeRao;
+
+% Define the intrinsic mechanical impedance for the device
+mass = simu.rho*hydro.properties.volume;
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(3,3,:))*simu.rho;
+radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(3,3,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
+hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(3,3)*simu.rho*simu.gravity;
+Gi = -((hydro.simulation_parameters.w)'.^2.*(mass+addedMass)) + 1j*hydro.simulation_parameters.w'.*radiationDamping + hydrostaticStiffness;
+Zi = Gi./(1j*hydro.simulation_parameters.w');
+
+% Calculate magnitude and phase for bode plot
+Mag = 20*log10(abs(Zi));
+Phase = (angle(Zi))*(180/pi);
+
+% Determine natural frequency based on the phase of the impedance
+[~,closestIndNat] = min(abs(imag(Zi)));
+natFreq = hydro.simulation_parameters.w(closestIndNat);
+natT = (2*pi)/natFreq;
+
+% Create bode plot for impedance
+figure()
+subplot(2,1,1)
+semilogx((hydro.simulation_parameters.w)/(2*pi),Mag)
+xlabel('freq (hz)')
+ylabel('mag (dB)')
+grid on
+xline(natFreq/(2*pi))
+xline(1/T,'--')
+legend('','Natural Frequency','Wave Frequency','Location','southwest')
+
+subplot(2,1,2)
+semilogx((hydro.simulation_parameters.w)/(2*pi),Phase)
+xlabel('freq (hz)')
+ylabel('phase (deg)')
+grid on
+xline(natFreq/(2*pi))
+xline(1/T,'--')
+legend('','Natural Frequency','Wave Frequency','Location','northwest')
+
+% Determine optimal latching time
+optDeclutchTime = 0.5*(natT - T)
+KpOpt = sqrt(radiationDamping(closestIndOmega)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(closestIndOmega)))^2)
+
+% Calculate the maximum potential power
+P_max = -sum(abs(Fexc).^2./(8*real(Zi)))
+
+

Controls/Declutching/userDefinedFunctions.m

@@ -0,0 +1,50 @@
+%Example of user input MATLAB file for post processing
+close all
+
+%Plot waves
+waves.plotElevation(simu.rampTime);
+try 
+    waves.plotSpectrum();
+catch
+end
+
+%Plot heave response for body 1
+output.plotResponse(1,3);
+
+%Plot heave forces for body 1
+output.plotForces(1,3);
+
+controllersOutput = controller1_out;
+signals = {'force','power'};
+for ii = 1:length(controllersOutput)
+    for jj = 1:length(signals)
+        controllersOutput(ii).(signals{jj}) =  controllersOutput(ii).signals.values(:,(jj-1)*6+1:(jj-1)*6+6);
+    end
+end
+
+% Plot controller power
+figure()
+plot(controllersOutput.time,controllersOutput.power(:,3))
+title('Controller Power')
+ylabel('Power (W)')
+xlabel('Time (s)')
+
+% Calculate average power for last 10 wave periods
+endInd = length(controllersOutput.time);
+startTime = controllersOutput.time(end) - 10*waves.period; % select last 10 periods
+[~,startInd] = min(abs(controllersOutput.time - startTime));
+disp('Controller Power:')
+mean(controllersOutput.power(startInd:endInd,3))
+
+% Plot declutching
+figure()
+yyaxis left
+plot(output.bodies.time,output.bodies.velocity(:,3))
+hold on
+xlabel('Time (s)')
+xlim([200 210])
+ylabel('Velocity (m/s)')
+ylim([-.4 .4])
+yyaxis right
+plot(output.bodies.time,output.bodies.forceExcitation(:,3)/1000)
+ylabel('Excitation Force (kN)')

Controls/Declutching/userDefinedFunctionsMCR.m

@@ -0,0 +1,28 @@
+%% User Defined Functions for MCR run
+controllersOutput = controller1_out;
+signals = {'force','power'};
+for ii = 1:length(controllersOutput)
+    for jj = 1:length(signals)
+        controllersOutput(ii).(signals{jj}) =  controllersOutput(ii).signals.values(:,(jj-1)*6+1:(jj-1)*6+6);
+    end
+end
+
+endInd = length(controllersOutput.time);
+startTime = controllersOutput.time(end) - 10*waves.period; % select last 10 periods
+[~,startInd] = min(abs(controllersOutput.time - startTime));
+
+mcr.meanPower(imcr) = mean(controllersOutput.power(startInd:endInd,3));
+mcr.meanForce(imcr) = mean(controllersOutput.force(startInd:endInd,3));
+
+mcr.maxPower(imcr) = max(abs(controllersOutput.power(startInd:endInd,3)));
+mcr.maxForce(imcr) = max(abs(controllersOutput.force(startInd:endInd,3)));
+
+if imcr == length(mcr.cases)
+    figure()
+    plot(mcr.cases,mcr.meanPower)
+    title('Mean Power vs. Declutching Time')
+    xlabel('Declutching Time (s)')
+    ylabel('Mean Power (W)')
+    xline(.4166, '--')
+    legend('','Theoretical Optimal')
+end

Controls/Declutching/wecSimInputFile.m

@@ -0,0 +1,37 @@
+%% Simulation Data
+simu = simulationClass();               % Initialize Simulation Class
+simu.simMechanicsFile = 'declutch.slx';      % 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 = 100;                    % Wave Ramp Time [s]
+simu.endTime = 400;                     % Simulation End Time [s]
+simu.solver = 'ode4';                   % simu.solver = 'ode4' for fixed step & simu.solver = 'ode45' for variable step 
+simu.dt = 0.01; 					    % Simulation time-step [s]
+simu.mcrMatFile = 'mcrCases.mat';
+
+%% Wave Information 
+
+% % Regular Waves  
+waves = waveClass('regular');           % Initialize Wave Class and Specify Type                                 
+waves.height = 2.5;                     % Wave Height [m]
+waves.period = 9.6664;                       % Wave Period [s]
+
+%% Body Data
+% Sphere
+body(1) = bodyClass('../hydroData/sphere.h5');          % Create the body(1) Variable
+body(1).geometryFile = '../geometry/sphere.stl';        % Location of Geomtry File
+body(1).mass = 'equilibrium';                           % Body Mass
+body(1).inertia = [20907301 21306090.66 37085481.11];   % Moment of Inertia [kg*m^2]     
+
+%% PTO and Constraint Parameters
+
+% Translational PTO
+pto(1) = ptoClass('PTO1');                      % Initialize PTO Class for PTO1
+pto(1).stiffness = 0;                           % PTO Stiffness [N/m]
+pto(1).damping = 0;                       % PTO Damping [N/(m/s)]
+pto(1).location = [0 0 0];                      % PTO Location [m]
+
+% Declutching Controller
+controller(1).declutching.declutchTime = 0.8;
+controller(1).declutching.Kp = 2.3202e5;

Controls/Latching/latchingTime.m

@@ -0,0 +1,14 @@
+function [Force, tLatch, tNormal, vel] = latchingTime(vel, tLatch, tNormal, velPrev, latchTime, dt, Kp, forceCoeff)
+
+% Latch when velocity is zero and release when latching time is reached
+if ((diff(sign([vel,velPrev]))~=0 && tLatch==0 && tNormal>.2) || (tLatch>0 && tLatch<latchTime))
+    Force = -forceCoeff*vel; %-forceOnBody;
+    tLatch = tLatch + dt;
+    tNormal = 0;
+else
+    Force = -Kp*vel;
+    tNormal = tNormal + dt;
+    tLatch = 0;
+end
+
+end

Controls/Latching/mcrBuildTimes.m

@@ -0,0 +1,11 @@
+close all;clear all;clc;
+%%
+mcr = struct();
+mcr.header = ["controller(1).latching.latchTime"];
+mcr.cases = zeros(9,1);
+%%
+time = linspace(2.0666,2.9666,25);
+%%
+mcr.cases = time';
+%%
+save mcrCases mcr

Controls/Latching/optimalTimeCalc.m

@@ -0,0 +1,74 @@
+% This script identifies the dynamics of the float in the respective wave 
+% conditions and determines the optimal proportional gain value for a 
+% passive controller (for regular waves)
+
+close all; clear all; clc;
+
+% Inputs (from wecSimInputFile)
+simu = simulationClass();
+body(1) = bodyClass('../hydroData/sphere.h5');
+waves.height = 2.5;
+waves.period = 9.6664;
+
+% Load hydrodynamic data for float from BEM
+hydro = readBEMIOH5(body.h5File, 1, body.meanDrift);
+
+% Define wave conditions
+H = waves.height;
+A = H/2;
+T = waves.period;
+omega = (1/T)*(2*pi);
+
+% Define excitation force based on wave conditions
+ampSpect = zeros(length(hydro.simulation_parameters.w),1);
+[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w));
+ampSpect(closestIndOmega) = A;
+FeRao = squeeze(hydro.hydro_coeffs.excitation.re(3,:))'*simu.rho*simu.gravity + squeeze(hydro.hydro_coeffs.excitation.im(3,:))'*simu.rho*simu.gravity*1j;
+Fexc = ampSpect.*FeRao;
+
+% Define the intrinsic mechanical impedance for the device
+mass = simu.rho*hydro.properties.volume;
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(3,3,:))*simu.rho;
+radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(3,3,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
+hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(3,3)*simu.rho*simu.gravity;
+Gi = -((hydro.simulation_parameters.w)'.^2.*(mass+addedMass)) + 1j*hydro.simulation_parameters.w'.*radiationDamping + hydrostaticStiffness;
+Zi = Gi./(1j*hydro.simulation_parameters.w');
+
+% Calculate magnitude and phase for bode plot
+Mag = 20*log10(abs(Zi));
+Phase = (angle(Zi))*(180/pi);
+
+% Determine natural frequency based on the phase of the impedance
+[~,closestIndNat] = min(abs(imag(Zi)));
+natFreq = hydro.simulation_parameters.w(closestIndNat);
+natT = (2*pi)/natFreq;
+
+% Create bode plot for impedance
+figure()
+subplot(2,1,1)
+semilogx((hydro.simulation_parameters.w)/(2*pi),Mag)
+xlabel('freq (hz)')
+ylabel('mag (dB)')
+grid on
+xline(natFreq/(2*pi))
+xline(1/T,'--')
+legend('','Natural Frequency','Wave Frequency','Location','southwest')
+
+subplot(2,1,2)
+semilogx((hydro.simulation_parameters.w)/(2*pi),Phase)
+xlabel('freq (hz)')
+ylabel('phase (deg)')
+grid on
+xline(natFreq/(2*pi))
+xline(1/T,'--')
+legend('','Natural Frequency','Wave Frequency','Location','northwest')
+
+% Determine optimal latching time
+optLatchTime = 0.5*(T - natT)
+KpOpt = radiationDamping(closestIndOmega)
+force = 80*(mass+addedMass(closestIndOmega))
+
+% Calculate the maximum potential power
+P_max = -sum(abs(Fexc).^2./(8*real(Zi)))
+
+

Controls/Latching/userDefinedFunctions.m

@@ -0,0 +1,50 @@
+%Example of user input MATLAB file for post processing
+close all
+
+%Plot waves
+waves.plotElevation(simu.rampTime);
+try 
+    waves.plotSpectrum();
+catch
+end
+
+%Plot heave response for body 1
+output.plotResponse(1,3);
+
+%Plot heave forces for body 1
+output.plotForces(1,3);
+
+controllersOutput = controller1_out;
+signals = {'force','power'};
+for ii = 1:length(controllersOutput)
+    for jj = 1:length(signals)
+        controllersOutput(ii).(signals{jj}) =  controllersOutput(ii).signals.values(:,(jj-1)*6+1:(jj-1)*6+6);
+    end
+end
+
+% Plot controller power
+figure()
+plot(controllersOutput.time,controllersOutput.power(:,3))
+title('Controller Power')
+ylabel('Power (W)')
+xlabel('Time (s)')
+
+% Calculate average power for last 10 wave periods
+endInd = length(controllersOutput.time);
+startTime = controllersOutput.time(end) - 10*waves.period; % select last 10 periods
+[~,startInd] = min(abs(controllersOutput.time - startTime));
+disp('Controller Power:')
+mean( mean(controllersOutput.power(startInd:endInd,3)))
+
+% Plot latching
+figure()
+yyaxis left
+plot(output.bodies.time,output.bodies.velocity(:,3))
+xlabel('Time (s)')
+xlim([200 220])
+ylabel('Velocity (m/s)')
+ylim([-10 10])
+hold on
+yyaxis right
+plot(output.bodies.time,output.bodies.forceExcitation(:,3)/1000)
+ylabel('Excitation Force (kN)')