Hi,
I am currently trying to design a PI pitch controller with openFAST and Simulink. Therefore, I am following the “Advanced Control Design for Wind Turbines”-paper ([url]https://www.nrel.gov/docs/fy08osti/42437.pdf[/url]).
To obtain the values for Ki and Kp I linearized my model at a desired operating point in region 3, linearizing every 10-degrees of azimuth. For the linearization I enabled the FlapDOF1, GenDOF and TwFADOF1 DOFs. Afterwards I successfully calculated the state-space and output matrices as well as the azimuth-averaged matrices, using the fx_mbc3.m script from the openFAST Matlab toolbox.
In the paper (p. 45/46) it states that one uses the azimuth averaged state matrices (AvgA, AvgB, etc.) to solve equation 4.1 and I also assume equation 3.8 (p. 29), is that correct?
The linearization in the paper results in the constants A, B, Bd (p.29). However, the fx_mbc.m-script outputs the azimuth averaged state 10x10 matrices. How do I process these, to obtain constants as in the paper?
Thanks and kind regards
Jannik
Dear Jason,
Thank you for your quick response.
As you suggested I reduced the the DOFs to the GenDOF only. However, I still encounter some issues.
- The MBC now fails with the error “fx_mbc3 WARNING: Number of blades is 0, not 3. MBC transformation was not performed.” I can still obtain the AvgA, AvgB, etc. matrices.
- The AvgA matrix is now resulting in [0,1;67.126,-5.777], if I understand you correctly the first column is supposed to be all zero (which is not the case) and can then be eliminated. What could be the issue here? And even if I could eliminate the first column, I would still have a vector with two values, how would I use this vector in eq. 3.8 from the paper?
- The AvgB matrix is a 2x9 matrix, what would be the procedure here?
Thanks for your help. I also attached my linearization files and a screenshot of my linearization-settings.
Kind regards
Jannik
model-windturbine-lin.zip (86.2 KB)
Dear Jannik,
I haven’t reviewed all your files, but just a few comments:
- When you don’t have any states in the rotating frame, then MBC3 does not apply, but you can still azimuth-average the matrices.
- As the forum post linked above discussed, you should eliminate the column associated with the generator-azimuth state (assuming it is virtually zero) and you can then eliminate the first row of A and B because it has nothing but zeros and the trivial relationship that generator speed equals itself (A(1,2) = 1). The A matrix will then be a scalar and B will be a row vector.
- If A(2,1) is nonzero (before eliminating rows and columns), I would guess the value is nonzero because of numerical roundoff. I would suggest plotting A(2,1) across azimuth angle (which is reported in x_op(1)) to see if the mean is close to zero, as I would expect it to be. Ensuring the solution is in periodic steady state before linearizing (e.g., virtually constant speed), ensuring that the linearization covers one full rotation of the rotor, and ensuring the output precision is high enough (e.g., OutFmt = “ES17.9E3”) should help.
- You can ignore any column in B that is not needed.
Best regards,
Dear Jason,
thanks for your timely response und very helpful comments. I have a few follow up questions though.
- As you suggested, I made sure that the rotor speed was constant before linearizing, that the linearization covers a full rotation and increased the output precision. I also reduced the time step. My azimuth-averaged A(2,1) is still nonzero. Do you have another suggestion? I plotted A(2,1) vs. azimuth angles (see attachment). Does it seem reasonable to you? How can I estimate how big the effect of a nonzero average A(2,1) would be?
- Which column of B corresponds to the generator speed? I could not find any resources which describe this.
Thank you and kind regards
Jannik
Dear Jannik,
From what I can see, the average value does look close to zero except for the presence of “spikes” at 60, 180, and 300 degrees azimuth, which I presume are associated with with the tower influence model. My guess is that the mean will get close to zero if the tower influence model is disabled. It looks likes enabling the tower influence model will result in some nonzero mean. I’m not sure how important tower influence is in your case, but if it is very important to you, I would probably suggest linearizing at a finer azimuth resolution to better resolve the “spike”.
The generator speed is not a column of B. The list of inputs making up the columns of B is described at the top of the linearization (*.lin) file in the section titled, “Order of inputs:”.
Best regards,
Dear Jason,
- I reduced the Azimuth to 2 degrees and was indeed able to obtain a value for A(2,1) which is close to zero. The “spike” did actually resolve in a downward and upward spike (see attachment).
- Thanks for the hint regarding the description of input ports. Is my assumption correct that I will have to use B(2,9), which is “collective blade-pitch command”, as input for B?
- As my baseline pitch controller model so far does not react as expected, I am trying to factor out possible mistakes. Regarding the linearization, would I have to linearize with ElastoDyn enabled or disabled? If so, I assume I would have to set the rated torque as constant?
Thanks for your help.
Kind regards
Jannik
Dear Jannik,
I’m glad you better resolved the “spike” which now results in a zero mean for A(2,1).
I agree that column 9 in the standard inputs is collective blade-pitch, which sounds like the input you need.
I’m not sure understand enough about what you are doing to be able to answer your last question. OpenFAST does not allow you to disable ElastoDyn altogether (although you can disable all of the ElastoDyn DOFs).
Best regards,
Dear Jason,
sorry for the confusion, I meant ServoDyn not ElastoDyn. So far, I linearized my model with ServoDyn disabled. To linearize at an operating point with rated generator speed I must set a high pitch angle due to the lacking generator torque. My reasoning was, that I might need to include the generator torque through ServoDyn. Is that correct?
Kind regards,
Jannik
Dear Jannik,
Yes, that sounds correct.
Best regards,
Dear Jason,
I am currently using the OpenFAST linearization tool to find linear models for control design. The models I am generating have no rotating states, and three rotating inputs (blade pitch 1/2/3).
In your post dated Fri Jun 18, 2021, you mention that MBC transforms do not apply when there are no rotating states. This makes sense when I look at the A matrix transform (for x_dot = Ax + Bu), however, I am confused by the B matrix transformation.
Equation (30) of the MBC3 user manual is:
B_NR = [T1^-1 0; 0 T1^-1] * B * T1c
My understanding is that the first matrix [T1^-1 0; 0 T1^-1] relates to the states - if there are no rotating states this is an identity matrix. The third matrix T1c is related to the inputs, again if there are no rotating inputs this is an identity matrix, in my case with three rotating inputs this should not be an identity matrix. When I look at the Matlab function fx_mbc3 it seems that if there are no rotating states but there are rotating inputs, the transformation relating to T1c won’t be applied. Is this the intended function, have I misunderstood some aspect of the MBC process?
Thanks,
Brendan
Dear Brendan,
I misspoke in my prior post. I should have said that MBC3 does not apply if you don’t have any rotating states, rotating inputs, or rotating outputs. However, MBC3 does apply if you have rotating states, rotating inputs, or rotating outputs.
Best regards,
Dear Jason,
Thank you for clarifying.
The Matlab function fx_mbc3 from the Matlab toolbox on GitHub will only perform the MBC transformations if it finds the number of blades equals 3 (this check is performed on line 60). The function solves for the number of blades by looking at the size of the matrix that contains the indices of the state triplets (line 53). Thus, it seems to me that the function will not perform any transformation if there are no rotating states even if there are rotating inputs and/or outputs. I think that line 60 should actually be checking the size of states, inputs, and outputs to determine which (if any) transformations should be applied. Does this make sense, or have I missed something elsewhere in the code?
Thanks,
Brendan
Dear Brendan,
I agree, it looks like you’ve found a bug in fx_mbc3.m. To be honest, I’m not sure I’ve ever ran MBC for a model without states. Can you raise this as an issue on github: github.com/OpenFAST/matlab-toolbox/issues? Or better yet, submit a pull request with the fix?
Best regards,
Dear Jason,
Sorry for the late reply, thank you for looking into it. It looks like it has been addressed by [url]https://github.com/OpenFAST/matlab-toolbox/pull/29[/url]. My apologies for not getting to this sooner.
Thanks for your help,
Brendan
Dear Brendan,
Yes, Bonnie Jonkman of Envision Energy kindly fixed the issues you raised through pull requests to the MATLAB Toolbox.
Best regards,
Dear Jason,
FIrstly, thank you for helping me run openfast on linux!
Recently, I am trying to learn design controller and reading “Advanced Control Design for Wind Turbines Part I” . I am also finding way to work out the constants A, B, Bd (p.29) in this report. But I don’t know whether to obtian the openfast file for CART. I want to use it to verify whether I can calculate these three constants with the method you said.
Best regards
Fuhang
Dear FuHang,
Unfortunately, the CART models are not available in the public domain, as discussed a bit in the following forum topic: Help running FAST v8 for the first time - #6 by Jason.Jonkman.
Best regards,
Oh, ok.
Thank you Jason.
Best regards,
Dear Jonkman,
I encountered the same problem with Jannik.Luehning. First, I followed your suggestion and switch off TwrAero, and set TwrPotent and TwrShadow=0, but it still has the spike in at 60, 180, and 300 degrees azimuth, and then I set NLinTimes=180 as Jannik.Luehening which he got good answer, but the spike existed for me. I am a litter puzzled. Could you please give me some advice? My OPENFAST Version is 3.0. After then, I set NLinTimes=3600, the spike still exsit.
Best regards!
