Merge pull request #40 from jtgrasb/Controls+PTO-Sim

Controls+pto sim

Controls/Declutching/README.md

@@ -0,0 +1,16 @@
+# WEC-Sim Controls Examples
+
+**Author:**          Jeff Grasberger (Sandia)
+
+**WEC-Sim Version:** v5.0 (or newer)
+
+**Matlab Version:** 2020b (or newer)
+
+**WEC-Sim Model**
+Numerical model for a semi-submerged sphere (diameter = 10 m)  with a declutching controller. "wecSim" can be typed into 
+the command window to run the example with the default setup. "optimalTimeCalc.m" is used to calculate the 
+optimal declutching time. The "mcrBuildTimes.m" script can be run to set up multiple conditions runs, then 
+"wecSimMCR" can be typed into the command window to run the different cases with varying declutching times.
+
+**Questions?**
+* Post all WEC-Sim modeling questions to the [WEC-Sim online forum](https://github.com/WEC-Sim/WEC-Sim/issues).

Controls/Latching/README.md

@@ -0,0 +1,16 @@
+# WEC-Sim Controls Examples
+
+**Author:**          Jeff Grasberger (Sandia)
+
+**WEC-Sim Version:** v5.0 (or newer)
+
+**Matlab Version:** 2020b (or newer)
+
+**WEC-Sim Model**
+Numerical model for a semi-submerged sphere (diameter = 10 m) with a latching controller. "wecSim" can be typed into 
+the command window to run the example with the default setup. "optimalTimeCalc.m" is used to calculate the 
+optimal latching time. The "mcrBuildTimes.m" script can be run to set up multiple conditions runs, then 
+"wecSimMCR" can be typed into the command window to run the different cases with varying latching times.
+
+**Questions?**
+* Post all WEC-Sim modeling questions to the [WEC-Sim online forum](https://github.com/WEC-Sim/WEC-Sim/issues).

Controls/MPC/README.md

@@ -0,0 +1,28 @@
+# WEC-Sim Controls Examples
+
+**Author:**          Jeff Grasberger (Sandia)
+
+**WEC-Sim Version:** v5.0 (or newer)
+
+**Matlab Version:** 2020b (or newer)
+
+**Dependencies:**
+
+   Optimization Toolbox 		            -->	for quadprog command
+   System Identification Toolbox	        -->	for ssdata, tfest, tfdata commands
+   Statistics and Machine Learning Toolbox  -->	for regress command
+   Control System Toolbox		            -->	for frd command
+   Symbolic Math Toolbox		            -->	for subs command
+
+**Description:**
+Numerical model for a semi-submerged sphere (diameter = 10 m) with a model predictive controller (MPC).
+"wecSim" can be typed into the command window to run the example with the default setup. 
+"plotFreqDep.m" solves for and plots the frequency dependent coefficients, 
+which are stored in "coeff.mat". "setupMPC.m" sets the controller up using "makePlantModel.m" 
+and "makePredictiveModel.m" and is called by the input file when "wecSim" 
+is run from the command window. "fexcPrediction.m" and "mpcFcn.m" predict the excitation forces 
+and solve the quadratic programming problem, respectively, and are both called by "sphereMPC.slx" 
+during the simulation. The model predictive controller parameters can be changed in the input file. 
+
+**Questions?**
+* Post all WEC-Sim modeling questions to the [WEC-Sim online forum](https://github.com/WEC-Sim/WEC-Sim/issues).

Controls/Passive (P)/README.md

@@ -0,0 +1,17 @@
+# WEC-Sim Controls Examples
+
+**Author:**          Jeff Grasberger (Sandia)
+
+**WEC-Sim Version:** v5.0 (or newer)
+
+**Matlab Version:** 2020b (or newer)
+
+**WEC-Sim Model**
+Numerical model for a semi-submerged sphere (diameter = 10 m) with a passive controller. 
+"wecSim" can be typed into the command window to run the example with the default setup. "optimalGainCalc.m" 
+is used to calculate the optimal proportional gain. The "mcrBuildGains.m" script can be run to set up 
+multiple conditions runs, then "wecSimMCR" can be typed into the command window to run the different 
+cases with varying proportional gain values.
+
+**Questions?**
+* Post all WEC-Sim modeling questions to the [WEC-Sim online forum](https://github.com/WEC-Sim/WEC-Sim/issues).

Controls/Reactive (PI)/README.md

@@ -0,0 +1,17 @@
+# WEC-Sim Controls Examples
+
+**Author:**          Jeff Grasberger (Sandia)
+
+**WEC-Sim Version:** v5.0 (or newer)
+
+**Matlab Version:** 2020b (or newer)
+
+**WEC-Sim Model**
+Numerical model for a semi-submerged sphere (diameter = 10 m) with a reactive controller. 
+"wecSim" can be typed into the command window to run the example with the default setup. "optimalGainCalc.m" 
+is used to calculate the optimal proportional and integral gains. The "mcrBuildGains.m" script can be run to set up 
+multiple conditions runs, then "wecSimMCR" can be typed into the command window to run the different 
+cases with varying proportional and integral gain values.
+
+**Questions?**
+* Post all WEC-Sim modeling questions to the [WEC-Sim online forum](https://github.com/WEC-Sim/WEC-Sim/issues).

Controls/Reactive (PI)/userDefinedFunctions.m

@@ -34,4 +34,4 @@
 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)))
+mean(controllersOutput.power(startInd:endInd,3))

Controls/Reactive (PI)/userDefinedFunctionsMCR.m

@@ -39,9 +39,9 @@
     ylabel('Integral Gain/Stiffness (N/m)');
     xlabel('Proportional Gain/Damping (Ns/m)');
     % Set color bar and color map
-    C = colorbar('location','EastOutside');
+    %C = colorbar('location','EastOutside');
     colormap(jet);
-    set(get(C,'XLabel'),'String','Power (Watts)')
+    %set(get(C,'XLabel'),'String','Power (Watts)')
     % Create title
     title('Mean Power vs. Proportional and Integral Gains');
 end

Controls/ReactiveWithPTO/README.md

@@ -0,0 +1,19 @@
+# WEC-Sim Controls Examples
+
+**Author:**          Jeff Grasberger (Sandia)
+
+**WEC-Sim Version:** v5.0 (or newer)
+
+**Matlab Version:** 2020b (or newer)
+
+**WEC-Sim Model**
+Numerical model for a semi-submerged sphere (diameter = 10 m) with a reactive controller and 
+simple direct drive power take-off. This example demonstrates the different controller gains required
+for electrical power maximization when compared to mechanical power maximization. "wecSim" can be 
+typed into the command window to run the example with the default setup. The 
+"mcrBuildGains.m" script can be run to set up multiple conditions runs, then "wecSimMCR" can be 
+typed into the command window to run the different cases with varying proportional and integral 
+gain values.
+
+**Questions?**
+* Post all WEC-Sim modeling questions to the [WEC-Sim online forum](https://github.com/WEC-Sim/WEC-Sim/issues).

Controls/ReactiveWithPTO/mcrBuildGains.m

@@ -0,0 +1,19 @@
+close all;clear all;clc;
+%%
+mcr = struct();
+mcr.header = ["controller(1).proportionalIntegral.Kp","controller(1).proportionalIntegral.Ki"];
+mcr.cases = zeros(36,2);
+%%
+Kp = linspace(3e5,5e5,6); Ki = linspace(-2e5,-8e4,6); % good for N = 100
+% Kp = linspace(1e5,1.75e5,6); Ki = linspace(-3e4,-.1e4,6); % good for N = 50
+Kpmat = []; Kimat = [];
+for jj = 1:length(Ki)
+    for ii = 1:length(Kp)
+        Kpmat        = [Kpmat;Kp(jj)];
+        Kimat        = [Kimat;Ki(ii)];
+    end
+end
+%%
+mcr.cases = [Kpmat,Kimat];
+%%
+save mcrCases mcr

Controls/ReactiveWithPTO/userDefinedFunctions.m

@@ -0,0 +1,117 @@
+%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
+
+%% Plots
+figure();
+plot(output.ptoSim(1).time,output.ptoSim(1).velocity)
+set(findall(gcf,'type','axes'),'fontsize',16)
+xlabel('Time (s)')
+ylabel('Velocity (rad/s)')
+title('PTO Velocity')
+grid on
+
+figure();
+plot(output.ptoSim(1).time,output.ptoSim(1).force)
+set(findall(gcf,'type','axes'),'fontsize',16)
+xlabel('Time (s)')
+ylabel('Torque (Nm)')
+title('PTO Shaft Torque')
+grid on
+
+figure();
+plot(output.ptoSim(1).time,output.ptoSim(1).current)
+set(findall(gcf,'type','axes'),'fontsize',16)
+xlabel('Time (s)')
+ylabel('Current (A)')
+title('Generator Current')
+grid on
+
+figure();
+plot(output.ptoSim(1).time,output.ptoSim(1).voltage)
+set(findall(gcf,'type','axes'),'fontsize',16)
+xlabel('Time (s)')
+ylabel('Voltage (V)')
+title('Generator Voltage')
+grid on
+
+figure();
+plot(output.ptoSim(1).time,output.ptoSim(1).I2RLosses)
+set(findall(gcf,'type','axes'),'fontsize',16)
+xlabel('Time (s)')
+ylabel('Losses (W)')
+title('I^2R Losses')
+grid on
+
+figure();
+plot(output.ptoSim(1).time,output.ptoSim(1).elecPower)
+set(findall(gcf,'type','axes'),'fontsize',16)
+xlabel('Time (s)')
+ylabel('Power (W)')
+title('Electrical Power')
+grid on
+
+% Calculate averages to make sure results line up:
+% average last few periods:
+endInd = length(controllersOutput.power);
+startTime = controllersOutput.time(end) - 10*waves.period; % select last 10 periods
+[~,startInd] = min(abs(controllersOutput.time(:) - startTime));
+
+meanControllerPower = mean(controllersOutput.power(startInd:endInd,3))
+meanMechPower = mean(output.ptoSim.absPower(startInd:endInd))
+meanInertiaPower = mean(output.ptoSim.inertiaForce(startInd:endInd).*output.ptoSim.velocity(startInd:endInd))
+meanDampingPower = mean(output.ptoSim.fricForce(startInd:endInd).*output.ptoSim.velocity(startInd:endInd))
+meanGenPower = mean(output.ptoSim.genForce(startInd:endInd).*output.ptoSim.velocity(startInd:endInd))
+meanElecPower = mean(output.ptoSim.elecPower(startInd:endInd))
+meanVIPower = mean(output.ptoSim.current(startInd:endInd).*output.ptoSim.voltage(startInd:endInd))
+meanI2RLosses = mean(output.ptoSim.I2RLosses(startInd:endInd))
+
+figure()
+labels = categorical({'Controller (Ideal)','Mechanical (Drivetrain)','Electrical (Generator)'});
+labels = reordercats(labels,{'Controller (Ideal)','Mechanical (Drivetrain)','Electrical (Generator)'});
+values = [meanControllerPower,meanMechPower,meanElecPower]/1000;
+b = bar(labels, values,'FaceColor',[.4 .8 .2]);
+ylabel('Power (kW)')
+xtickangle(45)
+
+b.FaceColor = 'flat';
+for ii = 1:length(values)
+    if values(ii) > 0
+        b.CData(ii,:) = [.8 .2 .2];
+    end
+end
+
+figure()
+labels = categorical({'Controller (Ideal)','Mechanical (Drivetrain)','Electrical (Generator)'});
+labels = reordercats(labels,{'Controller (Ideal)','Mechanical (Drivetrain)','Electrical (Generator)'});
+values = [-meanControllerPower/waves.power,-meanMechPower/waves.power,-meanElecPower/waves.power];
+b = bar(labels, values,'FaceColor',[.4 .8 .2]);
+ylabel('Capture Width (m)')
+xtickangle(45)
+
+b.FaceColor = 'flat';
+for ii = 1:length(values)
+    if values(ii) < 0
+        b.CData(ii,:) = [.8 .2 .2];
+    end
+end

Controls/ReactiveWithPTO/userDefinedFunctionsMCR.m

@@ -0,0 +1,87 @@
+%% 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.power);
+startTime = controllersOutput.time(end) - 5*waves.period; % select last 10 periods
+[~,startInd] = min(abs(controllersOutput.time(:) - startTime));
+
+mcr.meanControlPower(imcr) = mean(controllersOutput.power(startInd:endInd,3));
+mcr.meanMechPower(imcr) = mean(output.ptoSim.absPower(startInd:endInd));
+mcr.meanElecPower(imcr) = mean(output.ptoSim.elecPower(startInd:endInd));
+mcr.meanShaftTorque(imcr) = mean(output.ptoSim.force(startInd:endInd));
+
+mcr.maxControlPower(imcr) = max(abs(controllersOutput.power(startInd:endInd,3)));
+mcr.maxMechPower(imcr) = max(abs(output.ptoSim.absPower(startInd:endInd)));
+mcr.maxElecPower(imcr) = max(abs(output.ptoSim.elecPower(startInd:endInd)));
+mcr.maxShaftTorque(imcr) = max(abs(output.ptoSim.force(startInd:endInd)));
+
+if imcr == numel(mcr.cases(:,1))
+
+    % Kp and Ki gains
+    kps = unique(mcr.cases(:,1));
+    kis = unique(mcr.cases(:,2));
+
+    i = 1;
+    for kpIdx = 1:length(kps)
+        for kiIdx = 1:length(kis)
+            meanControlPowerMat(kiIdx, kpIdx) = mcr.meanControlPower(i);
+            meanMechPowerMat(kiIdx, kpIdx) = mcr.meanMechPower(i);
+            meanElecPowerMat(kiIdx, kpIdx) = mcr.meanElecPower(i);
+            meanForceMat(kiIdx, kpIdx) = mcr.meanShaftTorque(i);
+
+            maxControlPowerMat(kiIdx, kpIdx) = mcr.maxControlPower(i);
+            maxMechPowerMat(kiIdx, kpIdx) = mcr.maxMechPower(i);
+            maxElecPowerMat(kiIdx, kpIdx) = mcr.maxElecPower(i);
+            maxForceMat(kiIdx, kpIdx) = mcr.maxShaftTorque(i);
+            i = i+1;
+        end
+    end
+
+    % Plot surface for controller power at each gain combination
+    figure()
+    surf(kps,kis,meanControlPowerMat)
+    % Create labels
+    zlabel('Mean Controller (Ideal) Power (W)');
+    ylabel('Integral Gain/Stiffness (N/m)');
+    xlabel('Proportional Gain/Damping (Ns/m)');
+    % Set color bar and color map
+    %C = colorbar('location','EastOutside');
+    colormap(jet);
+    %set(get(C,'XLabel'),'String','Power (Watts)')
+    % Create title
+    title('Mean Controller Power vs. Proportional and Integral Gains');
+
+    % Plot surface for mechanical power at each gain combination
+    figure()
+    surf(kps,kis,meanMechPowerMat)
+    % Create labels
+    zlabel('Mean Mechanical (Drivetrain) Power (W)');
+    ylabel('Integral Gain/Stiffness (N/m)');
+    xlabel('Proportional Gain/Damping (Ns/m)');
+    % Set color bar and color map
+    %C = colorbar('location','EastOutside');
+    colormap(jet);
+    %set(get(C,'XLabel'),'String','Power (Watts)')
+    % Create title
+    title('Mean Mechanical Power vs. Proportional and Integral Gains');
+
+    % Plot surface for electrical power at each gain combination
+    figure()
+    surf(kps,kis,meanElecPowerMat)
+    % Create labels
+    zlabel('Mean Electrical (Generator) Power (W)');
+    ylabel('Integral Gain/Stiffness (N/m)');
+    xlabel('Proportional Gain/Damping (Ns/m)');
+    % Set color bar and color map
+    %C = colorbar('location','EastOutside');
+    colormap(jet);
+    %set(get(C,'XLabel'),'String','Power (Watts)')
+    % Create title
+    title('Mean Electrical Power vs. Proportional and Integral Gains');
+end