windturbine2_caseboh_control_matlabsimulink_windturbine_源码

% function to compute the average phase shift, "theta" function [sys, x0]=test_def_clk(t,x,u,flag,npts,ts) global tact; if flag == 0 % Initialization % Number of discrete states for storage: nstates = 2*npts + 6; sys = [0;nstates;2;2;0;0]; % Sizes of system (see SFUNTMPL) x0 = zeros(nstates/2,2); % Initial conditions x0(1,1) = 1; % Set counter to 1. x0(nstates/2,2) = 0; x0 = x0(:); elseif abs(flag) == 2 % Flag equals 2 on a real time hit % Is it a sample hit ? nstates = npts + 3; sys = zeros(nstates,2); sys(:) = x(:); % each 'nes' samples compute FFTs for tip speed ratio and power coefficient % extract the phase information for the two sequences % substact and average the result if rem(t + 1000*eps, ts) < 10000*eps % t is multiple of ts x(1,1) = x(1,1) + 1; sys(x(1,1),:) = u(:)'; div = x(1,1)/npts; if (div == round(div)) % blc = ones(npts,1); blc = blackman (npts); % avoid Gibbs effect % FFT ya = [sys(x(1,1)+1:npts+1,1);sys(2:x(1,1),1)]; ya = ya.*blc; yb = [sys(x(1,1)+1:npts+1,2);sys(2:x(1,1),2)]; yb= yb.*blc; ga = fft(ya,npts); gb = fft(yb,npts); a1 = real(ga); b1 = imag(ga); a2 = real(gb); b2 = imag(gb); % the phase spectral lines % arrange for [0,pi] for i=1:length(a1) at = atan(b1(i)/a1(i)); % at1= at; if (a1(i)>0), if (b1(i)>0) ar1(i)=at; elseif (b1(i)<0) ar1(i)=2*pi+at; else ar1(i)=0; end; elseif (a1(i)==0), if (b1(i)>0) ar1(i)=pi/2; elseif (b1(i)<0) ar1(i)=3*pi/2; else ar1(i)=0; end; elseif (b1(i)~=0), ar1(i)=at+pi; else ar1(i)=pi; end; at = atan(b2(i)/a2(i)); % at2= at; if (a2(i)>0), if (b2(i)>0) ar2(i)=at; elseif (b2(i)<0) ar2(i)=2*pi+at; else ar2(i)=0; end; elseif (a2(i)==0), if (b2(i)>0) ar2(i)=pi/2; elseif (b2(i)<0) ar2(i)=3*pi/2; else ar2(i)=0; end; elseif (b2(i)~=0), ar2(i)=at+pi; else ar2(i)=pi; end; % if % substract and average dfi(i) = ar1(i) - ar2(i); % 0 <= dfi <= pi % dfi(i) = at2 - at1 ; if ((dfi(i) < -pi)|(dfi(i) > pi)), dfi(i) = 2*pi - abs(dfi(i)); end; if ((dfi(i) >= -pi)&(dfi(i) < 0)), dfi(i) = abs(dfi(i)); end; end; % for dfi sys(nstates,2)=mean(dfi); tact=~tact; end; if sys(1,1) == npts x(1,1) = 1; end sys(1,1) = x(1,1); end; sys = sys(:); elseif flag == 3 sys(1) =x(length(x)) ; sys(2) =tact ; % outputs also the "theta" computing moment elseif flag == 4 % Return next sample hit ns = t/ts; % Number of samples sys = (1 + floor(ns + 1e-13*(1+ns)))*ts; else % Other flag options ignored sys = []; end % sfuntf