## FAST8 Linearization

This forum if for discussing controls. Questions about how to implement controls in FAST are more appropriate to the CAE Tools forum.

Moderator: Bonnie.Jonkman

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Jason

I run findBladeTriplets.m but I think also it doesn't work well. this is my result

ans =

[]

ans =

[]

ans =

4 5 6

function [Triplets, NTriplets] = findBladeTriplets(rotFrame,Desc )

%% Find the number of, and indices for, triplets in the rotating frame:

NTriplets = 0; % first initialize to zero
Triplets = [];
for i = 1:size(rotFrame) % loop through inputs/outputs
if rotFrame(i) % this is in the rotating frame
Tmp = zeros(1,3);
foundIt = false;
for chk = 1:size(chkStr)
BldNoCol = regexp(Desc{i},chkStr{chk},'match');
if ~isempty(BldNoCol)
foundIt = true;

str = regexp(Desc{i},BldNoCol{1},'split'); %this should return the strings before and after the match
FirstStr = [str{1} BldNoCol{1}(1:end-1) '.' ];
checkThisStr = [FirstStr str{2}];
k = str2double(BldNoCol{1}(end));
Tmp(k) = i;
break;
end % check the next one;
end

% find the other match values
if foundIt
for j = (i+1):size(rotFrame) % loop through all remaining control inputs
if rotFrame(j) % this is in the rotating frame
BldNoCol = regexp(Desc{j},checkThisStr,'match'); % match all but the blade number
if ~isempty(BldNoCol)
Num = regexp(Desc{j},FirstStr,'match'); % match all but the blade number
k = str2double(Num{1}(end));
Tmp(k) = j; % save the indices for the remaining blades
if ( all(Tmp ~=0) ) % true if all the elements of Tmp are nonzero; thus, we found a triplet of rotating indices
NTriplets = NTriplets + 1; % this is the number of control input triplets in the rotating frame
Triplets(NTriplets,:) = Tmp; % these are the indices for control input triplets in the rotating frame
break;
end
end
end
end % j - all remaining active control inputs
end

end % in the rotating frame
end % for i

return;
end

So, I just took the values that you wrote before ( NRotTripletOutput = 3 and RotTripletIndicesOutput = [2,3,4;5,6,7;8,9,10] ) and I run mbc3.m by it.

Regards,
Taha

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Taha,

I don't see that you've used the corrected findBladeTriplets.m file that I reported above. You can find the corrected file at this deeplink: https://github.com/old-NWTC/MBC/blob/ma ... Triplets.m. The use of this correct filed should produce correct results.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Jason,

Thank you, please could you check Getmatf8.m mbc3.m and findBladeTriples.m and let me know what is wrong in my case because I used the new code for findBladeTriples.m but still gives me wrong results.

I have a question regarding to MBC . If I do the MBC for the linearized model, it should give me 2 components of moments ( Cos, Sin) . In my case I have 3 outputs for each blade( blade root moments in 3 directions)( RootMxb1, RootMxb2, RootMxb3,RootMyb1, RootMyb2, RootMyb3 ,RootMzb1, RootMzb2, RootMzb3 ). So I will have ( RootMxb Cos, RootMxb sin, RootMyb cos, RootMyb sin, RootMzb Cos, RootMyb Sin) as outputs from MBC transformation which mean the outputs reduce from 9 outputs to 6 outputs but I am still having 9 outputs after I did MBC transformation and A,B,C,D matrices are in the same dimensions ( I think the columns of C and D matrices should be reduced by 3 but I still have the same dimensions like before ). then If I did the inverse MBC, it will give me 9 outputs again. Please, correct me if I am wrong.

1-
% GetMats_f8.m
% Written by J. Jonkman, NREL
% 19-July-2016: Updated to convert FAST v8.16 linearization files into
% format expected by mbc3.m
% NOTE that we assume all the files in FileNames contain the same data
% structures (same state matrices; same number of inputs, outputs, states,
% etc.; same rotor speed; but should have different azimuth angles.

% This m-file is used to read in the data written to multiple FAST linear output
% (.lin) files, compute the state matrix, [A], at each of the equally-spaced
% azimuth steps and their azimuth-average, along with their eigenvalues and
% eigenvectors.

format short g;

FileNames = {'Test18.1.lin','Test18.2.lin','Test18.3.lin','Test18.4.lin','Test18.5.lin','Test18.6.lin','Test18.7.lin','Test18.8.lin','Test18.9.lin','Test18.10.lin','Test18.11.lin','Test18.12.lin','Test18.13.lin','Test18.14.lin','Test18.15.lin','Test18.16.lin','Test18.17.lin','Test18.18.lin','Test18.19.lin','Test18.20.lin','Test18.21.lin','Test18.22.lin','Test18.23.lin','Test18.24.lin','Test18.25.lin','Test18.26.lin','Test18.27.lin','Test18.28.lin','Test18.29.lin','Test18.30.lin','Test18.31.lin','Test18.32.lin','Test18.33.lin','Test18.34.lin','Test18.35.lin','Test18.36.lin'};

if isempty(FileNames)
FileNames = {'Test18.1.lin','Test18.2.lin','Test18.3.lin','Test18.4.lin','Test18.5.lin','Test18.6.lin','Test18.7.lin','Test18.8.lin','Test18.9.lin','Test18.10.lin','Test18.11.lin','Test18.12.lin','Test18.13.lin','Test18.14.lin','Test18.15.lin','Test18.16.lin','Test18.17.lin','Test18.18.lin','Test18.19.lin','Test18.20.lin','Test18.21.lin','Test18.22.lin','Test18.23.lin','Test18.24.lin','Test18.25.lin','Test18.26.lin','Test18.27.lin','Test18.28.lin','Test18.29.lin','Test18.30.lin','Test18.31.lin','Test18.32.lin','Test18.33.lin','Test18.34.lin','Test18.35.lin','Test18.36.lin'};
end

% Input data:

MdlOrder = 1;
NAzimStep = length(FileNames);
data(NAzimStep) = ReadFASTLinear(FileNames{1}); %we'll read this twice so we can allocate space first; putting it at NAzimStep saves some reallocation later
N = data(NAzimStep).n_x;
NActvDOF = N / 2;
NInputs = data(NAzimStep).n_u;
NumOuts = data(NAzimStep).n_y;

clear xdop xop AMat DescCntrlInpt BMat DescOutput OutName CMat DMat AvgBMat AvgCMat AvgDMat;
clear RotTripletIndicesStates RotTripletIndicesCntrlInpt RotTripletIndicesOutput;

%% .................................
% allocate space for these variables
% ..................................
Azimuth = zeros(NAzimStep, 1);
Omega = zeros(NAzimStep, 1);
NDisturbs = 0;
if ( N > 0 )
DescStates = data(NAzimStep).x_desc;
xdop = zeros(N, NAzimStep);
xop = zeros(N, NAzimStep);
AMat = zeros(N, N, NAzimStep);
end

if ( NInputs > 0 )
DescCntrlInpt = data(NAzimStep).u_desc;
if (N>0)
BMat = zeros(N, NInputs, NAzimStep);
end
end

if ( NumOuts > 0 )
DescOutput = data(NAzimStep).y_desc;
if ( N > 0 )
CMat = zeros(NumOuts, N, NAzimStep);
end
if ( NInputs > 0 )
DMat = zeros(NumOuts, NInputs, NAzimStep);
end
end

%% .................................
% get data into variables expected by GetMats (concatenate data from
% different azimuths into matrices)
% ..................................

for iFile = 1:NAzimStep

Omega(iFile) = data(iFile).RotSpeed;
Azimuth(iFile) = data(iFile).Azimuth*180/pi;

if (isfield(data(iFile), 'A'))
AMat(:,:,iFile) = data(iFile).A;
end
if (isfield(data(iFile), 'B'))
BMat(:,:,iFile) = data(iFile).B;
end
if (isfield(data(iFile), 'C'))
CMat(:,:,iFile) = data(iFile).C;
end
if (isfield(data(iFile), 'D'))
DMat(:,:,iFile) = data(iFile).D;
end

if (isfield(data(iFile), 'x_op'))
xop(:,iFile) = cell2mat(data(iFile).x_op);
end
if (isfield(data(iFile), 'xdot_op'))
xdop(:,iFile) = cell2mat(data(iFile).xdot_op);
end

end

%% Find the azimuth-averaged linearized 1st order state matrices:
if ( isfield(data(1),'A') )
Avgxdop = mean(xdop,2);
Avgxop = mean(xop, 2);
AvgAMat = mean(AMat, 3);
end
if ( isfield(data(1),'B') )
AvgBMat = mean(BMat, 3);
end
if ( isfield(data(1),'C') )
AvgCMat = mean(CMat, 3);
end
if ( isfield(data(1),'D') )
AvgDMat = mean(DMat, 3);
end

%%
for i=1:NActvDOF
col = strfind(DescStates{i},'DOF_GeAz'); % find the starting index of the string 'DOF_GeAz'
if ( ~isempty(col) ) % true if the DescStates contains the string 'DOF_GeAz'
Omega(:) = xdop(i,:)';
break;
end
end
for i=1:NActvDOF
col = strfind(DescStates{i},'DOF_DrTr'); % find the starting index of the string 'DOF_DrTr'
if ( ~isempty(col) ) % true if the DescStates contains the string 'DOF_GeAz'
Omega(:) = Omega(:) + xdop(i,:)'; %This always comes after DOF_GeAz so let's just add it here (it won't get written over later).
break;
end
end

% ----------- Find multi-blade coordinate (MBC) transformation indices ----

%% Find the number of, and indices for, state triplets in the rotating
% frame:
NRotTripletStates = 0; % first initialize to zero
for i = 1:NActvDOF % loop through all active (enabled) DOFs
if data(1).x_rotFrame(i) % this is a state in the rotating frame

col = strfind(DescStates{i},'blade'); % find the starting index of the string 'blade'
if ~isempty(col) % true if the DescStates{I} contains the string 'blade'
k = str2double(DescStates{i}(col+6)); % save the blade number for the initial blade
Tmp = zeros(1,3); % first initialize to zero
Tmp(k) = i; % save the index for the initial blade
for j = (i+1):NActvDOF % loop through all remaining active (enabled) DOFs
if strncmp(DescStates{j},DescStates{i},col) % true if we have the same state from a different blade
k = str2double(DescStates{j}(col+6)); % save the blade numbers for the remaining blades
Tmp(k) = j; % save the indices for the remaining blades

if ( all(Tmp) ) % true if all the elements of Tmp are nonzero; thus, we found a triplet of rotating indices
NRotTripletStates = NRotTripletStates + 1; % this is the number of state triplets in the rotating frame
RotTripletIndicesStates(NRotTripletStates,:) = Tmp; % these are the indices for state triplets in the rotating frame
break;
end
end
end %for j
end

end
end % i - all active (enabled) DOFs

%% Find the number of, and indices for, control input triplets in the
% rotating frame:
if (NInputs > 0)
end

%% Find the number of, and indices for, output measurement triplets in the
% rotating frame:
if (NumOuts > 0 )
end

% ----------- Clear some unneeded variables -------------------------------
clear iFile i j k col Tmp;

% ----------- We're finished ----------------------------------------------

% Tell the user that we are finished:
disp( ' ' );
disp( 'GetMats_f8.m completed ' );
disp( 'Type "who" to list available variables' );

2-
function [Triplets, NTriplets] = findBladeTriplets(rotFrame,Desc )

%% Find the number of, and indices for, triplets in the rotating frame:

NTriplets = 0; % first initialize to zero
Triplets = [];
for i = 1:size(rotFrame) % loop through inputs/outputs
if rotFrame(i) % this is in the rotating frame
Tmp = zeros(1,3);
foundIt = false;
for chk = 1:size(chkStr)
BldNoCol = regexp(Desc{i},chkStr{chk},'match');
if ~isempty(BldNoCol)
foundIt = true;

str = regexp(Desc{i},BldNoCol{1},'split'); %this should return the strings before and after the match
FirstStr = [str{1} BldNoCol{1}(1:end-1) '.' ];
checkThisStr = [FirstStr str{2}];
checkThisStr = strrep(strrep(strrep(checkThisStr,')','\)'), '(', '\('),'^','\^'); %we need to get rid of the special characters
k = str2double(BldNoCol{1}(end));
Tmp(k) = i;
break;
end % check the next one;
end

% find the other match values
if foundIt
for j = (i+1):size(rotFrame) % loop through all remaining control inputs
if rotFrame(j) % this is in the rotating frame
BldNoCol = regexp(Desc{j},checkThisStr,'match'); % match all but the blade number
if ~isempty(BldNoCol)
Num = regexp(Desc{j},FirstStr,'match'); % match all but the blade number
k = str2double(Num{1}(end));
Tmp(k) = j; % save the indices for the remaining blades
if ( all(Tmp ~=0) ) % true if all the elements of Tmp are nonzero; thus, we found a triplet of rotating indices
NTriplets = NTriplets + 1; % this is the number of control input triplets in the rotating frame
Triplets(NTriplets,:) = Tmp; % these are the indices for control input triplets in the rotating frame
break;
end
end
end
end % j - all remaining active control inputs
end

end % in the rotating frame
end % for i

return;
end

3-
% function [...] = mbc(...) % later: convert this script into a function script
%
% Developed by Gunjit Bir, NREL (303-384-6953, gunjit_bir@nrel.gov)

% Last update: 08/30/2006
%----------------------------------------------------------------------------------------

% Objectives:
% 1. Given state-space matrices (A,B,Bd) and output matrices (C,D,Dd), defined partly in the
% rotoating frame and partly in the fixed frame, transform these matrices to the fixed
% coordinate frame using multi-blade coordinate transformation (MBC). The transformned
% matrices are MBC_A, MBC_B, MBC_Bd, MBC_C, MBC_D and MBC_Dd.
%
% 2. Given second-order matrices (M,Dmp,K), control input matrix (F), wind input disturbance
% matrix (Fd), displacement output matrix (Dc), and velocity output matrix (Vc), transform
% these to the fixed coordinate frame using multi-blade coordinate transformation (MBC).
% The transformned matrices are MBC_M, MBC_Dmp, MBC_K, MBC_F, MBC_Fd, MBC_Dc, and MBC_Vc.
%
% 3. Azimuth-average the MBC_A matrix and compute its eigenvalues and eigenvectors. The
% eigenvectors, referred to the fixed coordinates, are presented in both complex and
% amplitude-phase forms. The eigenvalues, also referred to the fixed coordinates, are
% presented in complex and damping-ratio/damped-frequency forms.

% ***WARNING: Do not copy or modify this code, or transfer algoritms (in original or modified form)
% from this code into another code without permission from NREL (Gunjit Bir, 303-384-6953).
%
% ***Disclaimer: This code is still in the developmental stage and no guarantee is given
% as to its proper functioning or accuracy of its results.

% ------------ INPUTS (these will be accessed later via the fuction arguments) ---
% ndof : total number of degrees of freedom
% NInputs : total num of control inputs
% NumOuts : total num of outputs
% ns : number of states
% A,B,Bd : matrices associated with the state-spece equation
% C,D,Dd : matrices associated with the ouputs equation
% M,Dmp,K : 2nd-order mass, stiffness, and damping/gyroscopic matrices
% F : control input matrix (associated with 2nd-order equation)
% Fd : wind input disturbance matrix (associated with 2nd-order equation)
% Dc, Vc : displacement and velocity output matrices (associated with 2nd-order equation)
% azm : vector of rotor azimuth positions wrt the the user-specified reference position [see AzimB1Up, FAST User's guide] (deg)
% omg : rotor speed (rad/sec)
% OmgDot : rotor acceleraton (rad/sec^2)

% RotTripletIndicesStates : State triplets in rotating frame (matrix of size rotating_dof_types*3)
% RotTripletIndicesCntrlInpt: Control-input triplets in rotating frame (matrix of size rotating_control_input_types*3)
% RotTripletIndicesOutput : Output triplets in rotating frame (matrix of size rotating_output_types*3)

% DescStates : description of states associated with input matrices (FAST-specific) %%

% --------------------- OUTPUTS ----------------------------------------------------
% MBC_A,MBC_B,MBC_Bd : state-space matrices transformed to the fixed frame
% MBC_C,MBC_D, MBC_Dd : output matrices transformed to the fixed frame
% MBC_M,MBC_Dmp,MBC_K : second-order mass, damping/gyroscopic and stiffness matrices transformed to the fixed frame
% MBC_F : control input matrix transformed to the fixed frame
% MBC_Fd : wind input disturbance matrix transformed to the fixed frame
% MBC_Dc,MBC_Vc : displacement and velocity output matrices (Vc) transformed to the fixed frame

% MBC_States : description of mbc-transformed states (FAST-specific) %%

% -----------------------------------------------------------------------------------------
%**NOTE 1: All inputs and output matrices may not be present. For example, user may supply or may be interested
% in multi-blade coordinate transformation of A matrix only. In this case, only MBC_A matrix along with
% fixed-frame eigensolution will be genertaed as outputs. The code checks for consistency and completness
% of selected inputs and generates associated outputs.

%**NOTE 2: The code currently assumes that the following inputs are available in the Matlab worspace
% (for example those generated by the FAST post-processor Eigenanalysis.m)

% AMat, BMat, BdMat, CMat, DMat, DdMat: 1st-order input matrices
% MassMat, DampMat, StffMat, FMat, FdMat, VelCMat, DspCMat: 2nd-order input matrices
% Omega : Vector of rotor speeds at specified azimuths (rad/sec)
% NAzimStep : number of azimuth steps at which i/o matrices are specified/generated
% NActvDOF : total number of degrees of freedom
% N : number of states
% Azimuth : vector of azimuth positions (in deg)
% -----------------------------------------------------------------------------------------
%% TBD: MBC flags/error messages / status messages / convert this script to function/

% error(nargchk(2, 3, nargin)) : implement later
%
% function [s, varargout] = mysize(x)
% msg = nargoutchk(1, 3, nargout);
% if isempty(msg)
% nout = max(nargout, 1) - 1;
% s = size(x);
% for k = 1:nout, varargout(k) = {s(k)}; end
% else
% disp(msg) & or, error('msg') ??
% end --------------------------------------------------------------------------------------

format short g;

ProgName = 'mbc3 (v1.00.00a-gbir, 15-Oct-2008)';

fprintf( '\n Running %s\n\n', ProgName );

% --- The following lines will be removed later ---
NAzSteps = NAzimStep;
ndof = NActvDOF;
ns = N;
azm = Azimuth;

new_seq_cont = 0;
new_seq_out = 0;

id_exist = exist('RotTripletIndicesStates');
if(id_exist == 0)
display('*** There are no rotating states. MBC transformation, therefore, cannot be performed.');
return
end

[nm,nb] = size(RotTripletIndicesStates);
if(nb ~= 3)
display('**ERROR: the number of colmn vectors in RotTripletIndicesStates must equal 3, the num of blades');
return
end
if(nm*nb > ndof)
display('**ERROR: the rotating dof exceeds the total num of dof');
return
end
new_seq_dof = get_new_seq(RotTripletIndicesStates,ndof);
new_seq_states = [new_seq_dof new_seq_dof+ndof];

mc=0; % number of rotating-frame control triplets
if (exist('RotTripletIndicesCntrlInpt'))
[mc,nb] = size(RotTripletIndicesCntrlInpt);
if(nb ~= 3)
display('**ERROR: the number of colmn vectors in RotTripletIndicesCntrlInpt must equal 3, the num of blades');
return
end
new_seq_cont = get_new_seq(RotTripletIndicesCntrlInpt,NInputs);
end

mo=0; % number of rotating-frame output triplets
if (exist('RotTripletIndicesOutput'))
[mo,nb] = size(RotTripletIndicesOutput);
if(nb ~= 3)
display('**ERROR: the number of colmn vectors in RotTripletIndicesOutput must equal 3, the num of blades');
return
end
new_seq_out = get_new_seq(RotTripletIndicesOutput,NumOuts);
end

ndfix = ndof-nm*nb; % fixed-frame dof
ncfix = NInputs-mc*nb; % fixed-frame controls
nofix = NumOuts-mo*nb; % fixed-frame outputs

if (exist('AMat'))
MBC_AvgA = zeros(ns);
end

if ( size(Omega) ~= NAzSteps)
display('**ERROR: the size of Omega vector must equal NAzSteps, the num of azimuth steps');
return
end
display('**ERROR: the size of OmegaDot vector must equal NAzSteps, the num of azimuth steps');
return
end

for iaz = 1:NAzSteps % begin azimuth loop

% compute azimuth positions of blades:
del_az = 2*pi/nb;
az1 = azm(iaz)*pi/180.0; % in radian
az2 = az1+del_az;
az3 = az2+del_az;

% get rotor speed and acceleration
omg = Omega(iaz);
omg2 = omg*omg;

% compute transformation matrices
c1 = cos(az1);
c2 = cos(az2);
c3 = cos(az3);
s1 = sin(az1);
s2 = sin(az2);
s3 = sin(az3);

cos_col = [c1;c2;c3];
sin_col = [s1;s2;s3];

tt = [ones(3,1), cos_col, sin_col];
ett = 1.5*sqrt(3);

cos231_col = [c2;c3;c1];
cos312_col = [c3;c1;c2];
sin231_col = [s2;s3;s1];
sin312_col = [s3;s1;s2];

tc1 = cos231_col.*sin312_col-sin231_col.*cos312_col;
ttv = [tc1, (sin231_col-sin312_col), (cos312_col-cos231_col)]'/ett;
tt2 = [zeros(3,1), -sin_col, cos_col];
tt3 = [zeros(3,1), -cos_col, -sin_col];

%---
T1 = [eye(ndfix)];
for ii = 1:nm
T1 = blkdiag(T1, tt);
end

T1v = [eye(ndfix)];
for ii = 1:nm
T1v = blkdiag(T1v, ttv);
end

T2 = [zeros(ndfix)];
for ii = 1:nm
T2 = blkdiag(T2, tt2);
end

T3 = [zeros(ndfix)];
for ii = 1:nm
T3 = blkdiag(T3, tt3);
end
%---
T1c = [eye(ncfix)];
for ii = 1:mc;
T1c = blkdiag(T1c, tt);
end

T1o = [eye(nofix)];
for ii = 1:mo
T1o = blkdiag(T1o, tt);
end

T1ov = [eye(nofix)];
for ii = 1:mo
T1ov = blkdiag(T1ov, ttv);
end

% mbc transformation of first-order matrices
% if ( MBC_EqnsOrder == 1 ) % activate later

if (exist('AMat'))
xAMat = row_col_xtion(AMat(:,:,iaz), new_seq_states, 1, 0); %--
xAMat = row_col_xtion(xAMat, new_seq_states, 2, 0); %--

AK = xAMat(ndof+1:ns,1:ndof);
AC = xAMat(ndof+1:ns,ndof+1:ns);
MBC_A(:,:,iaz) = [zeros(ndof), eye(ndof);
T1v*(AK*T1+omg*AC*T2-omg2*T3-OmgDot*T2), T1v*(AC*T1-2*omg*T2)];

MBC_A(:,:,iaz) = row_col_xtion(MBC_A(:,:,iaz), new_seq_states, 2, 1); %--
MBC_A(:,:,iaz) = row_col_xtion(MBC_A(:,:,iaz), new_seq_states, 1, 1); %--
MBC_AvgA = MBC_AvgA + MBC_A(:,:,iaz);

clear xAMat; %--
end

if (exist('BMat'))
xBMat = row_col_xtion(BMat(:,:,iaz), new_seq_states, 1, 0); %--
if (new_seq_cont ~= 0); xBMat = row_col_xtion(xBMat, new_seq_cont, 2, 0); end %--

B1 = xBMat(1:ndof,:);
B2 = xBMat(ndof+1:ns,:);
MBC_B(:,:,iaz) = [T1v*B1;T1v*B2]*T1c;

if (new_seq_cont ~= 0); MBC_B(:,:,iaz) = row_col_xtion(MBC_B(:,:,iaz), new_seq_cont, 2, 1); end %--
MBC_B(:,:,iaz) = row_col_xtion(MBC_B(:,:,iaz), new_seq_states, 1, 1); %--

clear xBMat; %--
end

if (exist('BdMat'))
xBdMat = row_col_xtion(BdMat(:,:,iaz), new_seq_states, 1, 0); %--

Bd1 = xBdMat(1:ndof,:);
Bd2 = xBdMat(ndof+1:ns,:);
MBC_Bd(:,:,iaz) = [T1v*Bd1;T1v*Bd2];

MBC_Bd(:,:,iaz) = row_col_xtion(MBC_Bd(:,:,iaz), new_seq_states, 1, 1); %--

clear xBdMat; %--
end

if (exist('CMat'))
xCMat = CMat(:,:,iaz);
if (new_seq_out ~= 0); xCMat = row_col_xtion(CMat(:,:,iaz), new_seq_out, 1, 0); end %--
xCMat = row_col_xtion(xCMat, new_seq_states, 2, 0); %--

MBC_C(:,:,iaz) = T1ov*xCMat*[T1, zeros(ndof); omg*T2, T1];

MBC_C(:,:,iaz) = row_col_xtion(MBC_C(:,:,iaz), new_seq_states, 2, 1); %--
if (new_seq_out ~= 0); MBC_C(:,:,iaz) = row_col_xtion(MBC_C(:,:,iaz), new_seq_out, 1, 1); end %--

clear xCMat; %--
end

if (exist('DMat'))
xDMat = DMat(:,:,iaz);
if (new_seq_out ~= 0); xDMat = row_col_xtion(DMat(:,:,iaz), new_seq_out, 1, 0); end %--
if (new_seq_cont ~= 0); xDMat = row_col_xtion(xDMat, new_seq_cont, 2, 0); end %--

MBC_D(:,:,iaz) = T1ov*xDMat*T1c;

if (new_seq_cont ~= 0); MBC_D(:,:,iaz) = row_col_xtion(MBC_D(:,:,iaz), new_seq_cont, 2, 1); end %--
if (new_seq_out ~= 0); MBC_D(:,:,iaz) = row_col_xtion(MBC_D(:,:,iaz), new_seq_out, 1, 1); end %--

clear xDMat; %--
end

if (exist('DdMat'))
xDdMat = DdMat(:,:,iaz);
if (new_seq_out ~= 0); xDdMat = row_col_xtion(DdMat(:,:,iaz), new_seq_out, 1, 0); end %--
MBC_Dd(:,:,iaz) = T1ov*xDdMat;

if (new_seq_out ~= 0); MBC_Dd(:,:,iaz) = row_col_xtion(MBC_Dd(:,:,iaz), new_seq_out, 1, 1); end %--

clear xDdMat; %--
end

% end

% mbc transformation of second-order matrices
% if ( MBC_EqnsOrder == 2 ) %% activate later

if (exist('MassMat'))
xMassMat = row_col_xtion(MassMat(:,:,iaz), new_seq_dof, 1, 0); %--
xMassMat = row_col_xtion(xMassMat, new_seq_dof, 2, 0); %--
MBC_M(:,:,iaz) = xMassMat*T1;
MBC_M(:,:,iaz) = row_col_xtion(MBC_M(:,:,iaz), new_seq_dof, 2, 1); %--
MBC_M(:,:,iaz) = row_col_xtion(MBC_M(:,:,iaz), new_seq_dof, 1, 1); %--

xDampMat = row_col_xtion(DampMat(:,:,iaz), new_seq_dof, 1, 0); %--
xDampMat = row_col_xtion(xDampMat, new_seq_dof, 2, 0); %--
MBC_Dmp(:,:,iaz) = 2*omg*xMassMat*T2+xDampMat*T1;
MBC_Dmp(:,:,iaz) = row_col_xtion(MBC_Dmp(:,:,iaz), new_seq_dof, 2, 1); %--
MBC_Dmp(:,:,iaz) = row_col_xtion(MBC_Dmp(:,:,iaz), new_seq_dof, 1, 1); %--

xStffMat = row_col_xtion(StffMat(:,:,iaz), new_seq_dof, 1, 0); %--
xStffMat = row_col_xtion(xStffMat, new_seq_dof, 2, 0); %--
MBC_K(:,:,iaz) = omg2*xMassMat*T3+OmgDot*xMassMat*T2+omg*xDampMat*T2+xStffMat*T1;
MBC_K(:,:,iaz) = row_col_xtion(MBC_K(:,:,iaz), new_seq_dof, 2, 1); %--
MBC_K(:,:,iaz) = row_col_xtion(MBC_K(:,:,iaz), new_seq_dof, 1, 1); %--

clear xMassMat xDampMat xStffMat; %--
end

if (exist('FMat'))
xFMat = row_col_xtion(FMat(:,:,iaz), new_seq_dof, 1, 0); %--
if (new_seq_cont ~= 0); xFMat = row_col_xtion(xFMat, new_seq_cont, 2, 0); end %--

MBC_F(:,:,iaz) = xFMat*T1c;

if (new_seq_cont ~= 0); MBC_F(:,:,iaz) = row_col_xtion(MBC_F(:,:,iaz), new_seq_cont, 2, 1); end %--
MBC_F(:,:,iaz) = row_col_xtion(MBC_F(:,:,iaz), new_seq_dof, 1, 1); %--

clear xFMat; %--
end

if (exist('FdMat'))
xFdMat = row_col_xtion(FdMat(:,:,iaz), new_seq_dof, 1, 0); %--

MBC_Fd(:,:,iaz) = xFdMat;

MBC_Fd(:,:,iaz) = row_col_xtion(MBC_Fd(:,:,iaz), new_seq_dof, 1, 1); %--

clear xFdMat; %--
end

if (exist('DspCMat'))
xVelCMat=VelCMat(:,:,iaz);
if (new_seq_out ~= 0); xVelCMat = row_col_xtion(VelCMat(:,:,iaz), new_seq_out, 1, 0); end %--
xVelCMat = row_col_xtion(xVelCMat, new_seq_dof, 2, 0); %--
MBC_Vc(:,:,iaz) = T1ov*xVelCMat*T1;
MBC_Vc(:,:,iaz) = row_col_xtion(MBC_Vc(:,:,iaz), new_seq_dof, 2, 1); %--
if (new_seq_out ~= 0); MBC_Vc(:,:,iaz) = row_col_xtion(MBC_Vc(:,:,iaz), new_seq_out, 1, 1); end %--

xDspCMat=DspCMat(:,:,iaz);
if (new_seq_out ~= 0); xDspCMat = row_col_xtion(DspCMat(:,:,iaz), new_seq_out, 1, 0); end %--
xDspCMat = row_col_xtion(xDspCMat, new_seq_dof, 2, 0); %--
MBC_Dc(:,:,iaz) = T1ov*(omg*xVelCMat*T2+xDspCMat*T1);
MBC_Dc(:,:,iaz) = row_col_xtion(MBC_Dc(:,:,iaz), new_seq_dof, 2, 1); %--
if (new_seq_out ~= 0); MBC_Dc(:,:,iaz) = row_col_xtion(MBC_Dc(:,:,iaz), new_seq_out, 1, 1); end %--

clear xDspCMat xVelCMat; %--
end

% end

end % end of azimuth loop
%---------------------------- Eigensolution ---------------------------------------

if (exist('AMat'))
MBC_AvgA = MBC_AvgA/NAzSteps; % azimuth-average of azimuth-dependent MBC_A matrices

% --- Eigensolution of the Azimuth-averaged MBC_AvgA matrix (in the fixed frame)
% call function 'eianalysis' to obtain the following
% nss : number of re-sequenced states
% EigenVects : nsXns
% Evals : nsX1
% eigenVects : ndofXnss
% eigenVals : nssX1
% DampedFrequencies : nssX1
% DampedFreqs_Hz : nssX1
% NaturalFrequencies : nssX1
% NaturalFreqs_Hz : nssX1
% DecrementRate : nssX1
% DampRatios : nssX1
% MagnitudeModes : ndofXnss
% PhaseModes_deg : ndofXnss

[EigenVects,eigenVects,Evals,eigenVals,DampedFrequencies,DampedFreqs_Hz,NaturalFrequencies,...
NaturalFreqs_Hz,DecrementRate,DampRatios,MagnitudeModes,PhaseModes_deg,nss] = eianalysis(MBC_AvgA);

MBC_EigenVects = EigenVects;
MBC_Evals = Evals;
MBC_eigenVects = eigenVects;
MBC_eigenVals = eigenVals;

MBC_DampedFrequency = DampedFrequencies;
MBC_DampedFrequencyHz = DampedFreqs_Hz;
MBC_NaturalFrequency = NaturalFrequencies;
MBC_NaturalFrequencyHz = NaturalFreqs_Hz;
MBC_DecrementRate = DecrementRate;
MBC_DampingRatio = DampRatios;
MBC_ModeShapeMagnitude = MagnitudeModes;
MBC_ModeShapePhaseDeg = PhaseModes_deg;

end
% ---- Label fixed-coordinate states following mbc transformation
% if (exist('DescStates'))
% MBC_States = DescStates;
% for is = ndfix+1:3:ndof;
% k = strfind(MBC_States{is},'DOF');
% MBC_States{is} = [MBC_States{is}(1:k-1) 'collective DOF'];
% k = strfind(MBC_States{is+1},'DOF');
% MBC_States{is+1} = [MBC_States{is+1}(1:k-1) 'cosine DOF'];
% k = strfind(MBC_States{is+2},'DOF');
% MBC_States{is+2} = [MBC_States{is+2}(1:k-1) 'sine DOF'];
% end
% end

% --- Clear non-required data and display successful completion of MBC
clear omg NAzSteps nb ndof ndof_rot ns azm iaz ndfix ncfix nofix nm mc mo omg2 OmgDot;
clear T1 T2 T3 T1v T1c T1o T1ov;
clear cos_col sin_col tt ett c1 c2 c3 s1 s2 s3 del_az az1 az2 az3 is k;
clear tc1 cos231_col ttv cos312_col tt2 sin231_col tt3 sin312_col ii;
if (exist('Bd1'))
clear Bd1;
end
if (exist('Bd2'))
clear Bd2;
end
if (exist('B1'))
clear B1;
end
if (exist('B2'))
clear B2;
end
if (exist('AK'))
clear AK;
end
if (exist('AC'))
clear AC;
end

% ----------- Clear unneeded variables -------------------------------
clear EigenVects eigenVects Evals eigenVals DampedFrequencies DampedFreqs_Hz NaturalFrequencies;
clear NaturalFreqs_Hz DecrementRate DampRatios MagnitudeModes PhaseModes_deg;
clear new_seq_dof new_seq_states;

if (exist('RotTripletIndicesCntrlInpt'))
clear new_seq_cont;
end

if (exist('RotTripletIndicesOutput'))
clear new_seq_out;
end

disp(' ');
disp(' Multi-Blade Coordinate transformation completed ');
%-----------------------------------------------------------
%%%xmbc =[MBC_DecrementRate MBC_DampedFrequencyHz MBC_DampingRatio*100 MBC_NaturalFrequencyHz]'

Best Regards
Taha

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Taha,

I cannot review your code in that format. Besides, I've sent you all of the files I've used, which should work to obtain proper results.

MBC3 does not reduce the number of terms from 3 (one for each blade) to 2 (cos, sin). Instead, MBC3 converts 3 (one for each blade in the rotating frame) to 3 (collective, cos, sin in the fixed frame).

I hope that helps.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Jason,

I would like to thank you for your help.

I have just a small question, After I did the linearization. I estimated the states using an estimator and it works good for my linearized model and the nonlinear FAST turbine. I have implemented LQR but I find that LQR is working good with the linear model where the nonlinear FAST turbine it gives me this error
Error reported by S-function 'FAST_SFunc' in 'OpenLoop/FAST Nonlinear Wind Turbine/S-Function':

I understand that this error comes from large change in the pitch input. I reduced my step pitch input as I can but I am still having the same error.

I am thinking that this error comes from high feedback gain of LQR, I tried to change in Q and R matrices of LQR but it doesn't help.

So Please Can you tell me from where this error comes?

Best Regards
Taha Fouda

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Taha,

You can avoid receiving errors about "no valid value of phi" by upgrading AeroDyn v15, as reported in the following forum topic: viewtopic.php?f=4&t=1752&p=8379.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Jason,

I have read the mentioned post and I updated AeroDyn v15 but I still face the same problem. I think instability may be the reason because I have checked the poles of the linear model and they are

AvgeigenVals =

-0.117 + 2.0848i
-0.00079865 + 0i
-1.0902 + 0i
-2.6852 + 4.0457i
-2.6969 + 3.7327i
-2.7048 + 3.7429i

and after I have done MBC, they are
MBC_eigenVals =

-2.6911 + 4.9911i
-2.6852 + 4.0468i
-2.6832 + 2.4638i
-0.14414 + 2.0887i
6.3093e-06 + 0i
-1.0915 + 0i

I am wondering why there is one postivie pole in the RHS after I did MBC?

these are my DOF
True FlapDOF1 - First flapwise blade mode DOF
True GenDOF - Generator DOF
True TwFADOF1 - First fore-aft tower bending-mode DOF

Regards
Taha

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Taha,

As described in this forum topic: viewtopic.php?f=13&t=621, rigid-body modes (i.e. modes without stiffness) show up in MBC3 as a pair of zero-valued (or near-zero-valued) frequencies with +/- inf damping (i.e., eigenvalues with real values only). That is, each rigid-body mode will introduce an additional mode beyond the number of enabled DOFs and the damping is unphysical. The rigid-body mode should not cause a problem in FAST.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Jason,

Sorry, I was off last days, now I am working on solving this error ( FAST_Solution:FAST_AdvanceStates:AD_UpdateStates:BEMT_UpdateStates(node 10, blade 3):BEMT_UnCoupledSolve:DeterminePhiBounds:There is no valid value of phi for these operating conditions! Vx = -3.1687, Vy = 51.327, rlocal = 32.36, theta = 1.0028) I have upgraded the Aerodyn V15 input file from this link https://github.com/old-NWTC/AeroDyn/find/master as you refer me before, but I am still facing the same issue.

in my case, I have estimated the nonlinear fast turbine states and I am using LQR as a full state feedback controller. this error appears when I enable the feedback gain. So please could you help me to solve this issue. I tried many possibilities, reducing the pitch command input, Changing in the LQR parameters but no change.

Regards
Taha

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Taha,

It is not the input files that change; instead, the AeroDyn v15 source code needs to be updated and FAST v8 needs to be recompiled to fix the issue regarding no valid values of phi.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Wan.Y
Posts: 1
Joined: Tue Oct 10, 2017 7:31 am
Organization: University of Shanghai for Sicence and Technology
Location: China

### Re: FAST8 Linearization

Dear Taha.Fouda

Sorry to bother you! I have a question for you. In the linearization, how do you control the number of input in the FAST V8 similar to fastv7

Yours sincerely
Wan.Y

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Wan,

In FAST 8, you cannot choose which number of inputs as you want but you can partion the B and D matrix for which inputs you want. e.g choose the columns and rows that represent the inputs which you want and delete the other.

Regards
Taha

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Wan,

I agree with Taha.

The input parameter LinInputs from the FAST v8 primary input file controls which system inputs are included in the linear state-space model. See the FAST v8 ReadMe file for more information: https://wind.nrel.gov/nwtc/docs/README_FAST8.pdf.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Taha.Fouda
Posts: 44
Joined: Thu Jun 22, 2017 1:51 am
Organization: DLR Institiue of Flight Systems
Location: Germany

### Re: FAST8 Linearization

Dear Jason,

I need to calculate the fatigue damage for the blade and the tower, so Could you please suggest me a matlab script to calculate the random vector A from the time simulation to be used for rainflow analysis.

the random vector A comes from calculating the peaks and valleys. It is the first time for me to deal with this issue. I was going to use LIFE 2 that has been supported by FAST, but I have used my own turbulence model that differs from FAST turbulence.

Regards,
Taha

Jason.Jonkman
Posts: 5734
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: FAST8 Linearization

Dear Taha,

I'm sorry, but I don't really understand your question. I don't know what you mean by the "random vector A" or by "LIFE 2".

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov