Wind turbine controls - MATLAB: S-Function 'FAST_SFunc' does not exist

Wind & Water Controls
1

Dear all,

I used to run these controls using FAST, but due to a configuration error about the version of MATLAB they are outputing an error.
Can you help me to create a block diagram on Simulink that I can run these controls again?

[code]if px==1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optimal control

                Num = 2*r{i}*(u - 2.14*omega_g^2)*(-2*2.14*omega_g) ...
                    +2*q{i}*(Pref - u*omega_g)*(-u);
                DerivCp=geraDerivadaCp(Cp,wind,Pitch,omega_g);
                
                Parte1=r{i}*(u - 2.14*omega_g^2)^2;
                Parte2=q{i}*(Pref - u*omega_g)^2;
                
                Au2=b2*q{i}*omega_g/DerivCp;
                Bu2= -2*b2*r{i}*2.14*omega_g/DerivCp - b2*q{i}*Pref/DerivCp + r{i} + q{i}*omega_g^2;
                Cu2= -( -2*b2*r{i}*2.14^2*omega_g^3/DerivCp + r{i}*2.14*omega_g^2 + q{i}*Pref*omega_g );
                raizes = roots([Au2 Bu2 Cu2]);
                raizMin = min(raizes);
                raizMax = max(raizes);
                if raizMin>0
                    ut3=raizMin;
                else
                    ut3=raizMax;
                end
                u=ut3;
                flag=['markov'];
                if u<0
                    disp('socorro. u negativo')
                end
                flag = ['markov'];
            end
            
            
            
            if px==2
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Jockerman interpretado Yolanda
                if (leer<1129)
                    u = 2.14*omega_g^2;
                elseif (leer>1165)
                    u = Pref/omega_g;
                else
                    alpha = -leer/36 + 1165/36;
                    u = alpha*2.14*omega_g^2 + (1-alpha)*Pref/omega_g;
                end
                flag=['jocker'];
            end
            
            
            if px==3
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Regra chaveado Chino
                if (wind<11.4)
                    u = 2.14*omega_g^2;
                else
                    u = Pref/omega_g;
                end
                flag=['chino'];
            end[/code]

This is the block diagram it was running until the MATLAB version changes.

https://drive.google.com/file/d/1PBhWsKo4T62aR8A4Bj_pLO70dG70ihiq/view

The script now is presenting the following errors:

I’ve searched about these errors I found the following topic:

https://www.mathworks.com/matlabcentral/answers/259203-s-function-does-not-exist

Then I considered that the problem can be the MATLAB version, but still couldn’t install any of them that could run the script.

Can you recommend a block diagram that I can run these already implemented controls?

I’m really grateful for your help.

2

Dear Jonathan,

The error, “S-Function ‘FAST_SFunc’ does not exist” has been discussed on our forum before e.g. see: S-Function 'FAST_SFunc' does not exist - #10 by Jason.Jonkman.

Is there a reason you are running FAST v7 instead of upgrading to FAST v8?

Best regards,

3

Thanks you for your kind consideration.

I followed the instructions to add FAST_SFunc.mexw64 to MATLAB path.
However, now it begins to open the model. Then it displays a message ‘Compiling’ in the bottom and MATLAB suddenly closes.

These files were done during FAST v7, although I could upgrade to FAST v8, just I’m not sure how can I join my files on FAST v8.

Sincerely,
Jonathan

4

Dear Jonathan,

I’m not familiar with the problem you are reporting. What script are you trying to run when MATLAB closes?

The FAST v8 includes documentation and scripts to help you convert a model from FAST v7 to FAST v8.

Best regards,

5

Dear Professor Jonkman,

I was running scriptEnergiesFast.m:

[code]% o arquivo scriptMonteCarlo_WindTurbine.m
% e EnergiesFAST.m estao relacionados com esse aqui
clear all, close all, clc,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% precisa desse arquivo para gerar a funcao numerica de Cp
exemple_DadosWind

% the FAST CertTest folder must be entered below:
CertTest_Dir = ‘C:\Users\Yolanda\CertTest’;

% This script needs the files OpenLoop.mdl, FAST_Sfunc.mexw32,
% and Read_FAST_Input.m located somewhere on the Matlab path so I’ll
% add the Simulink directories in the FAST archive:

addpath(CertTest_Dir);

for px=1:3

%input_fast=['Energies.fst'];
input_fast=['NRELOffshrBsline5MW_YVS.fst'];

% Read FAST input file and set initial conditions
Read_FAST_Input

global leer;
global Pe;
global tiempoGlobal;
global wind;
global DeriCp;
global Pitch;

tiempo = 1150;

lect_usb=0;
open_system('Energies');
%set_param('Energies','StartFcn','localAddEventListener');
set_param('Energies','SimulationCommand','Start');
set_param('Energies', 'StopTime', num2str(tiempo))

pause(2)
h1=add_exec_event_listener('Energies/omega_g','PostOutputs',@datoleer);
h2=add_exec_event_listener('Energies/Pe','PostOutputs',@datoPe);
h3=add_exec_event_listener('Energies/Wind','PostOutputs',@datoWind);
h4=add_exec_event_listener('Energies/DeriCp','PostOutputs',@datoDeriCp);
h5=add_exec_event_listener('Energies/Pitch','PostOutputs',@datoPitch);

u=0; Pref=5000000;

Nit= tiempo*12;
figure(1);clf;
h=plot(nan,nan,'r-','LineWidth',2);
%hold on
%h2=plot(nan,nan,'b-','LineWidth',2);
%grid,xlim([1 Nit]),ylim([0 2.5*Pref/1000]),
grid,xlim([1 Nit]),ylim([1000 2000]),
hold off


omega_g=0.105*leer;
while isempty(omega_g)
    omega_g=0.105*leer;
    pause(0.1);
end

X{1}=[0;11]; X{2}=[11;11.2]; X{3}=[11.2;11.4]; X{4}=[11.4;11.6]; X{5}=[11.6;12]; X{6}=[12;Inf];
% %r{1}=2.5; r{2}=1.8; r{3}=1.6; r{4}=0.4; r{5}=0.2; r{6}=1e-9; %funciona
% %bastante pero no completo
%r{1}=6; r{2}=4;    r{3}=3;   r{4}=2;   r{5}=0.2;  r{6}=1e-9;
% q{1}=0; q{2}=0.01; q{3}=0.1; q{4}=0.2; q{5}=0.95; q{6}=1;
%r{1}=2.5; r{2}=1.8; r{3}=1.6; r{4}=1; r{5}=0.05; r{6}=1e-9;
%q{1}=0;  q{2}=0.05; q{3}=0.1; q{4}=0.2; q{5}=0.95; q{6}=1;

%% prova nova metodo randomico
q{1} =  0;   q{2} = 0.6463;  q{3} =  0.7094;  q{4} =  0.7547;  q{5} =  0.2760; q{6} =  1;
r{1} = 1.5904; r{2} = 0.1869; r{3} =0.4898; r{4} =0.4456; r{5} = 0.05; r{6} = 1.0000e-09;

% esse funciona mais ou menos ateh certo ponto
%r{1}=2.5; r{2}=1.8; r{3}=1.6; r{4}=0.1; r{5}=0.05; r{6}=1e-9;
%q{1}=0;  q{2}=0.05; q{3}=0.1; q{4}=0.2; q{5}=0.95; q{6}=1;


% cut-in; cut-out; 3 m/s, 11.4 m/s, 25 m/s
% area 12,445.3 m2
b2=(1/115926000);

u=0;
vecTheta=[];
for k=1:Nit
    omega_g=0.105*leer;
    
    Pref_int=Pref/1000;
    
    for i=1:max(size(X))
        if (X{i}(1)<=wind) && (wind<X{i}(2))
            vecTheta(k)=i;
            
            if px==1
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Optimal control
                
                Num = 2*r{i}*(u - 2.14*omega_g^2)*(-2*2.14*omega_g) ...
                    +2*q{i}*(Pref - u*omega_g)*(-u);
                DerivCp=geraDerivadaCp(Cp,wind,Pitch,omega_g);
                
                Parte1=r{i}*(u - 2.14*omega_g^2)^2;
                Parte2=q{i}*(Pref - u*omega_g)^2;
                
                Au2=b2*q{i}*omega_g/DerivCp;
                Bu2= -2*b2*r{i}*2.14*omega_g/DerivCp - b2*q{i}*Pref/DerivCp + r{i} + q{i}*omega_g^2;
                Cu2= -( -2*b2*r{i}*2.14^2*omega_g^3/DerivCp + r{i}*2.14*omega_g^2 + q{i}*Pref*omega_g );
                raizes = roots([Au2 Bu2 Cu2]);
                raizMin = min(raizes);
                raizMax = max(raizes);
                if raizMin>0
                    ut3=raizMin;
                else
                    ut3=raizMax;
                end
                u=ut3;
                flag=['markov'];
                if u<0
                    disp('socorro. u negativo')
                end
                flag = ['markov'];
            end
            
            
            
            if px==2
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Jockerman interpretado Yolanda
                if (leer<1129)
                    u = 2.14*omega_g^2;
                elseif (leer>1165)
                    u = Pref/omega_g;
                else
                    alpha = -leer/36 + 1165/36;
                    u = alpha*2.14*omega_g^2 + (1-alpha)*Pref/omega_g;
                end
                flag=['jocker'];
            end
            
            
            if px==3
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Regra chaveado Chino
                if (wind<11.4)
                    u = 2.14*omega_g^2;
                else
                    u = Pref/omega_g;
                end
                flag=['chino'];
            end
            
            break
        end
    end
    
    u_int= u;
    
    vecOmega_g(k)=inv(0.105)*omega_g;
    vecU(k)=u_int;
    vecPe(k)=Pe;
    vecPref(k)=Pref_int;
    vecWind(k)=wind;
    set_param('Energies/tau','Gain',num2str(u_int));
    set_param('Energies/Pref','Gain',num2str(Pref_int));
    vectiempoGlobal(k)=tiempoGlobal;
    pause (1e-6);
    set(h,'xdata',[1:k],'ydata',[vecOmega_g([1:k])]);
    pause(0.1);
    
end

set_param('Energies','SimulationCommand','Stop');
TS = sum(vectiempoGlobal(2:end)-vectiempoGlobal(1:end-1))/Nit
totalTime = vectiempoGlobal(end)-vectiempoGlobal(1)
 
vecPower=OutData(:,1);
mean_Power = mean(vecPower(1:end-5))
std_Power= std(vecPower(1:end-5))

vecClock=clock;
save(['dataMarkov_FAST_5MW' date sprintf('_%ih%im.mat',[vecClock(4:5)])]);

end

% mean_Power =
%
% 4.3693e+03
%
%
% std_Power =
%
% 150.0559

%multiplica r_i por 10 pra obter isso abaixo
%mean_Power = 4.371513071698347e+03
%std_Power = 1.374370323910432e+02
[/code]

I had read the documentation to convert the models, I’ll follow the procedures to run the code in the newer version.

Wish you have a good day,
Jonatan