pmec-osu / fostwin Goto Github PK

@bretbosma The code that calculates the excitation force time series for SID, based on S(w), seems a little convoluted and may contain some errors (not entirely sure yet).

A potentially simpler version that should yield the same result is this (to be checked)

  % generate spectrum with random phases
    phi = 2*pi*rand(size(S)); % random phases
    A = sqrt(2*S*df); %NOT complex A - phase is dealt with further below .*exp(1i*phi); % complex A
    
    % interpolate wamit results to fit spectral frequency
    wamitExAft = interp1(wamit.w,wamit.FexAftPitch,w,'spline','extrap');
    wamitExBow = interp1(wamit.w,wamit.FexBowPitch,w,'spline','extrap');
    
    % multiply JONSWAP amplitude spectrum and wamit excitation spectrum
    FexAft_w = A .* wamitExAft;
    FexBow_w = A .* wamitExBow;

    t = 0:TsTwin:duration; % around 100 hz for tsTwin

    % Note: WecSim uses Re*cos(x) - Im*sin(x), which corresponds to
    % real(A*exp(-ix))
    FexAft = zeros(size(t));
    FexBow = zeros(size(t));
    for n=1:length(w)
        arg = -1i*(w(n)*t + phi(n));
        FexAft = FexAft + real(FexAft_w(n)*exp(arg));  
        FexBow = FexBow + real(FexBow_w(n)*exp(arg));  
    end

When running the simulation in realtime, there are often issues when running the pTg.stop line to stop the simulation on the speedgoat.

Symptoms:

• run pTg.stop or stoptarget (when coming from the UI)
• Get Cannot stop application on target 'EGIBaseline2' : Timed-out waiting for application to stop
• Simulation is no longer transmitting data - also the current step time stops increasing if looking at the speedgoat screen/ using the UI - no data coming in anymore
• About 30 seconds later, the speedgoat ends up in the "DONE" state and the logsout variable is transferred to the host machine

Adding Simulink.sdi.setRecordData(false); didn't seem to fix the issue

@bretbosma I believe I found a first issue as to why SID and wecSim give different results. There are likely more differences, but this particular issue should probably be resolved first.

I compared the excitation force from the "excitationWAMIT.mat" file with that from the h5 file used in wecSim. The Matlab code used for the comparison is pasted at the end of this issue. This code will need wecSim on the path for the readBEMIOH5 function.

The results:

Observations:

• the comparison looks good for the abs
• for the real and imag, it seems the phase is "scrambled" (for lack of a better word) for the aft flap for the SID case. The behavior of this excitation force does not look physical. I suspect something went wrong when the "excitationWAMIT.mat" file was created

Possible actions:

• Replace the "excitationWAMIT.mat" file with a corrected version (if found to be incorrect)
• Remove the "excitationWAMIT.mat" file, and read from the h5 file as in the code below. If going this route, we should extract the required h5 read code from the wecSim example, and have a little helper function in our repo (to avoid needing a wecSim install on the path for the SID option).

@bretbosma - any views?

% check excitation force between WAMIT mat file for SID and h5 file for
% WecSim

clearvars; close all; clc;

% SID excitation force from WAMIT
sid = load('excitationWAMIT.mat');

dofPitch = 5;
rho = 1020;
g = 9.81;

% use wecSim script to read H5 file
hydroWecSimFlapBow = readBEMIOH5('hydroData/foswec.h5',1,0); % bow flap is body 1
hydroWecSimFlapAft = readBEMIOH5('hydroData/foswec.h5',3,0); % aft flap is body 3

wecSimFexBowPitch =  rho*g* (squeeze(hydroWecSimFlapBow.hydro_coeffs.excitation.re(dofPitch,1,:)) ...
                        + 1i*squeeze(hydroWecSimFlapBow.hydro_coeffs.excitation.im(dofPitch,1,:)));

wecSimFexAftPitch =  rho*g* (squeeze(hydroWecSimFlapAft.hydro_coeffs.excitation.re(dofPitch,1,:)) ...
                        + 1i*squeeze(hydroWecSimFlapAft.hydro_coeffs.excitation.im(dofPitch,1,:)));

figure
subplot(2,1,1)
plot(sid.w, abs(sid.FexBowPitch))
hold on
plot(hydroWecSimFlapBow.simulation_parameters.w, abs(wecSimFexBowPitch))
ylabel('Fex Pitch Bow (Nm/m)')
xlabel('t (s)')
title('Abs')
legend('SID','wecSim')

subplot(2,1,2)
title('Abs')
plot(sid.w, abs(sid.FexAftPitch))
hold on
plot(hydroWecSimFlapAft.simulation_parameters.w, abs(wecSimFexAftPitch))
ylabel('Fex Pitch Aft (Nm/m)')
xlabel('t (s)')
legend('SID','wecSim')

figure
subplot(2,1,1)
plot(sid.w, real(sid.FexBowPitch))
hold on
plot(hydroWecSimFlapBow.simulation_parameters.w, real(wecSimFexBowPitch))
ylabel('Fex Pitch Bow (Nm/m)')
xlabel('t (s)')
title('Real')
legend('SID','wecSim')

subplot(2,1,2)
plot(sid.w, real(sid.FexAftPitch))
hold on
plot(hydroWecSimFlapAft.simulation_parameters.w, real(wecSimFexAftPitch))
ylabel('Fex Pitch Aft (Nm/m)')
xlabel('t (s)')
legend('SID','wecSim')

figure
subplot(2,1,1)
plot(sid.w, imag(sid.FexBowPitch))
hold on
plot(hydroWecSimFlapBow.simulation_parameters.w, imag(wecSimFexBowPitch))
ylabel('Fex Pitch Bow (Nm/m)')
xlabel('t (s)')
title('Imag')
legend('SID','wecSim')

subplot(2,1,2)
plot(sid.w, imag(sid.FexAftPitch))
hold on
plot(hydroWecSimFlapAft.simulation_parameters.w, imag(wecSimFexAftPitch))
ylabel('Fex Pitch Aft (Nm/m)')
xlabel('t (s)')
legend('SID','wecSim')

@bretbosma The JONSWAP spectrum is calculated here:

This appears to be missing a scaling factor that adjusts the overall spectrum energy to account for gamma. Based on IEC TC114 -2 this factor should be of the form:

For gamma = 3.3 this introduces ~25 error in the forcing / surface elevation.

@Graham-EGI - Can you take care of this?

Hello @johannes-spinneken, I'm a PhD student working on this competition and I have a question about the frequency response in your paper (Fig. 7).

I'm not sure what is the input and the output here, I guess this is based on the flap frame, and shows the frequency response from flap torque to flap velocity, is that correct? I tried to compare this to the bode plot I got from your SystemID model. The result shows that the two systems are totally different. Are we supposed to get the same system in this competition?

Besides, I tried to find out the frequency response of your wecSim model, and this seems to be different from the SystemID model (but somewhat similar). So I am also wondering how much difference is expected between the two models.

Thank you for your time and response.

From today's information session: Can you add a flag that the controller can see to indicate the state machine state?

@bretbosma @ryancoe

The linear and quadratic damping used for the WecSim FOSTWIN both seem a little high.

Is there any documentation / reporting available for the data analysis that led to those coefficients?

I made a rough start on the description for the docs for the ctrl state machine, but know it needs some work.

Section begins here:

From today's information session: Can we change the time step?

Not sure if this is the behavior we want, but when restarting a simulation through the UI, where the model isn't re-compiled, the number of faults/ unstable states that the controller is allowed to into before landing in the fault state isn't re-set.

Might become an issue if someone is using the realtime system to test optimizing parameters, where some of the tests land the controller in an unstable state.

From today's information session: Can we change/affect the drive train dynamics?

From today's information session: What do you use to detect instability and trigger the state machine?

@bretbosma The way the FOSWTIN code is currently setup, gamma is set to 3.3 for SID, and not set explicitly for the WecSim simulations.

WecSim will calculate gamma based on Tp and Hs using IEC TC114 -2:

For the standard case (Hs = 0.136m and Tp=2.61s), WecSim will use gamma of 1.0, whereas gamma = 3.3 (hard coded) for SID. This means WecSim and SID may often use different gamma and results are difficult to compare.

Solutions:

• Could bring gamma into the main init script, and then explicitly set gamma for both SID and WecSim
• Could calculate gamma in the SID wave generator script, based on the same equations as WecSim

The purpose of this issue is to track comparisons between the FOSTWIN SID and WecSim results. This issue covers Case A: Hs=0.136m and Tp = 2.61m