Dear Jason
Thank you. I did MBC and I run mbc3.m but I got this error
**ERROR: the number of colmn vectors in RotTripletIndicesOutput must equal 3, the num of blades
Dear Taha,
I’m not sure why you are receiving this error, but I searched and other people on this forum have received the same error and resolved it, see:
Linearization issue using Fast V7.02 - #15 by Rana.Mamdouh.
Best regards,
Dear Jason,
Sorry, I cannot reach why I have this error. but as written in one of the forum you show me in the previous comment, I have these values.
size(RotTripletIndicesOutput)
ans =
0 0
size(DescOutput)
ans =
13
I am not sure if my linearization process was going right or wrong. So please could you take a look for the attached file if I did it right or wrong. and correct me for the value of (LinTimes). I want to linearize at region 3 with steady wind 18 m/s. rotor speed 12.1 rpm.
Regards
Taha
Linearization Taha.zip (12.2 KB)
Dear Taha,
It looks like you linearization should include 3 outputs for each blade: RootMxb1-3, RootMyb1-3 and RootMzb1-3. Are these outputs in your linearization output file?
I’m not sure I understand your question regarding LinTimes, but you are still not linearizing every 10 degrees for 12.1 rpm.
Best regards,
Dear Jason,
Yes, I have 3 outputs for each blade in my linearized output file.
as the rotor is spinning by 12.1 rpm = 72,6 degree/sec so the time required for every 10 degree azimuth is approximately 0.14 sec.
I adjusted Lintimes to be like this:-
36 NLinTimes
0.14,0.28,0.42,0.56,0.7,0.84,0.98,1.12,1.26,1.4,1.54,1.68,1.82,1.96,2.1,2.24,2.38,2.52,2.66,2.8,2.94,3.08,3.22,3.36,3.5,3.64,3.78,3.92,4.06,4.2,4.34,4.48,4.62,4.76,4.9,5.04 LinTimes
and this gives me 36 linearized output files but I don’t know How can I upload more than 1 linearized file in Getmatf8.m.
actually, I started again with Generator degree of freedom and generator speed output, so I don’t need to make MBC and I added the operating points values to the linearized model. the response of the linearized model becomes nearly to FAST non linear but the nonlinear response always periodic ( oscillate around the generator speed ) I don’t why?
Then I added more degree of freedoms but in this case I have to do MBC and it gives me the error I mentioned before . I cleared RotTripletIndicesOutput and I think this can be wrong but just obtain MBC matrices but in this case I need to do the inverse of MBC for the output but I don’t know which m file can do that?
Regards,
taha
Dear Taha,
Yes, your LinTimes look correct (although I would start outputting the linearized solution after all start-up transients have died out e.g. by setting the first LineTimes to be 60.14 rather than 0.14). On line 17 of GetMats_f8.m is an example showing how MBC3 can be used to process two linearization outputs; you can easily change this to your 36 files.
I’m not sure I can comment on the periodic behavior of the nonlinear model without knowing more about your simulation set-up.
I don’t know why you are getting the error. I suggest stepping through the MATLAB script as it executes to debug.
NREL hasn’t provided the scripts to perform the inverse MBC, but these should not be difficult to implement yourself.
Best regards,
Dear Jason,
Thank you, I did the linearization and I adjusted Lintimes as you said then I changed line 17 in Getmat_f8.m from 2 linearized input to 36 input files but when I run it , it gives me NAzimStep = 12 not 36. I attached here Getmat_f8.m and mbc3.m, Please could you check them and let me know if there is something wrong.
Actually I don’t have mechanical background and it’s my first time to deal with FAST, So please could you give me a hint How can I do the inverse of MBC for the output?
Regards,
Taha
GetMat.rar (104 KB)
Dear Taha,
I briefly reviewed your GetMats_f8.m file and found several typos on line 17, plus several other changes on other lines. I’ve corrected these in the attached update. Running this version produces mostly correct results.
By mostly correct, I noticed that the corrected GetMats_f8.m file generated incorrect values for NRotTripletOutput and RotTripletIndicesOutput. The output was NRotTripletOutput = 0 and RotTripletIndicesOutput = . However, the correct result should be NRotTripletOutput = 3 and RotTripletIndicesOutput = [2,3,4;5,6,7;8,9,10]. I asked Bonnie Jonkman–who originally developed GetMats_f8.m while working at NREL, but now works at Envision Energy–about this, and she traced the problem back down to a bug in the findBladeTriplets.m file included with MBC3. She has now corrected this bug in the updated version of MBC3 available from here: github.com/old-NWTC/MBC.
The inverse MBC transformation is described in Eq. (4) in the User’s Guide to MBC3 that is provided with the source code.
GetMats_f8_JJ.m.txt (7.31 KB)
I hope that helps.
Best regards,
Dear Jason,
Thank you and sorry for bothering you. when I run findBladeTriplets.m, I need to enter these two parameters ( rotFrame,Desc ) which I didn’t get them from Getmatf8.m.
regards,
Taha
Dear Taha,
The “rotFrame” variables (separate variables for States/x, CntrlInpts/u, and Output/y) are contained with the “data” structure. There are separate “Desc” variables also for States, CntrlInpt, and Output.
Best regards,
Dear Jason
I run findBladeTriplets.m but I think also it doesn’t work well. this is my result
findBladeTriplets(data(18).y_rotFrame, data(18).y_desc)
ans =
[]
findBladeTriplets(data(18).x_rotFrame, data(18).x_desc)
ans =
[]
findBladeTriplets(data(18).u_rotFrame, data(18).u_desc)
ans =
4 5 6
this is findBladeTriplets.m code:-
function [Triplets, NTriplets] = findBladeTriplets(rotFrame,Desc )
%% Find the number of, and indices for, triplets in the rotating frame:
chkStr = {‘[Bb]lade \d’, ‘[Bb]lade [Rr]oot \d’, ‘[Bb]\d’, ‘\d’};
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
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: github.com/old-NWTC/MBC/blob/ma … Triplets.m. The use of this correct filed should produce correct results.
Best regards,
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);
OmegaDot = 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
data(iFile) = ReadFASTLinear(FileNames{iFile});
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,:)‘;
OmegaDot( = xdop(i+NActvDOF,:)’;
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).
OmegaDot( = OmegaDot( + xdop(i+NActvDOF,:)’;
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)
[RotTripletIndicesCntrlInpt, NRotTripletCntrlInpt] = findBladeTriplets(data(1).u_rotFrame,DescCntrlInpt );
end
%% Find the number of, and indices for, output measurement triplets in the
% rotating frame:
if (NumOuts > 0 )
[RotTripletIndicesOutput, NRotTripletOutput] = findBladeTriplets(data(1).y_rotFrame,DescOutput );
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:
chkStr = {‘[Bb]lade \d’, ‘[Bb]lade [Rr]oot \d’, ‘[Bb]\d’, ‘\d’};
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
% MBC: Multi-Blade Coordinate Transformation for a turbine with 3-blade rotor
%
% 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_types3)
% RotTripletIndicesCntrlInpt: Control-input triplets in rotating frame (matrix of size rotating_control_input_types3)
% 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)
% OmegaDot : Vector of rotor accelerations at specified azimuths (rad/sec2)
% 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;
%---------- Multi-Blade-Coordinate transformation -------------------------------------------
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-nmnb; % fixed-frame dof
ncfix = NInputs-mcnb; % 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
if ( size(OmegaDot) ~= NAzSteps)
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;
OmgDot = OmegaDot(iaz);
% 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
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,
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’:
FAST_Solution:FAST_AdvanceStates:AD_UpdateStates:BEMT_UpdateStates(node 10, blade 2):BEMT_UnCoupledSolve:DeterminePhiBounds:There is no valid value of phi for these operating conditions! Vx = -7.59392E-02, Vy = 19.774, rlocal = 32.263, theta = 0.97725
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
Dear Taha,
You can avoid receiving errors about “no valid value of phi” by upgrading AeroDyn v15, as reported in the following forum topic: FAST v8 AeroDyn v15 convergence problems for low-speed, high-turbulence - #2 by Jason.Jonkman.
Best regards,
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
Dear Taha,
As described in this forum topic: Learizing Baseline 5MW Wind Turbine with FAST, 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,
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 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
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,