Question about PI control
Moderator: Bonnie.Jonkman

 Posts: 3
 Joined: Tue Jul 04, 2017 11:12 am
 Organization: self
 Location: Italy
Question about PI control
hello to everyone
Im trying to implement in matlab the PI control of a Nrel 5 MW wind turbine and i wanted to ask if anyone know which are the equations to find Kp and Ki with the pitch angle beta
Or at least some values for some beta so that i can find an equation myself
Thank you
Im trying to implement in matlab the PI control of a Nrel 5 MW wind turbine and i wanted to ask if anyone know which are the equations to find Kp and Ki with the pitch angle beta
Or at least some values for some beta so that i can find an equation myself
Thank you

 Posts: 3857
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Question about PI control
Dear Vlad,
The proportional and integral gains in the NREL 5MW baseline controller are defined by Eq. (718)  (7.20) in the specifications document for the NREL 5MW baseline turbine: http://www.nrel.gov/docs/fy09osti/38060.pdf.
Best regards,
The proportional and integral gains in the NREL 5MW baseline controller are defined by Eq. (718)  (7.20) in the specifications document for the NREL 5MW baseline turbine: http://www.nrel.gov/docs/fy09osti/38060.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
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

 Posts: 3
 Joined: Tue Jul 04, 2017 11:12 am
 Organization: self
 Location: Italy
Re: Question about PI control
Yes those i already saw
I was asking if someone had already them tabulated or calculated
Also i saw online that instead of these 2 equations there is the chance to use the equation:
Kp(i)=Kp(theta 0)* GK
with GK=1/(1+theta(i)/theta_k) with thetak=6.30°
Is that valid?
and if it is valid,how do i write the code for the new theta knowing kp and ki?
Best regards
I was asking if someone had already them tabulated or calculated
Also i saw online that instead of these 2 equations there is the chance to use the equation:
Kp(i)=Kp(theta 0)* GK
with GK=1/(1+theta(i)/theta_k) with thetak=6.30°
Is that valid?
and if it is valid,how do i write the code for the new theta knowing kp and ki?
Best regards

 Posts: 3857
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Question about PI control
Dear Vlad,
Yes, the equation you state is equivalent to what is in the paper I linked to, with:
Kp(theta 0) = 2*I_Drivetrain*Omega_0*zeta_phi*omega_phin/(N_Gear*dP/dtheta(theta 0))
But I'm not sure I understand your last question about how to calculate theta.
Best regards,
Yes, the equation you state is equivalent to what is in the paper I linked to, with:
Kp(theta 0) = 2*I_Drivetrain*Omega_0*zeta_phi*omega_phin/(N_Gear*dP/dtheta(theta 0))
But I'm not sure I understand your last question about how to calculate theta.
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
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

 Posts: 3
 Joined: Tue Jul 04, 2017 11:12 am
 Organization: self
 Location: Italy
Re: Question about PI control
if i have
n=50 with a step=1
the part of the code i'm interested in is this:
if w<wmax then e=0 %where e is the error i need to correct with pi controller
else
e=wwmax
Kp(1)=0.01882681
Ki(1)=0.008068634
bk=6.302336
Gk(i)=1/(1+b(i1)/bk)
Kp(i) =Kp(1)*Gk(i)
Ki(i) =Ki(1)*Gk(i)
so i'll have a :
b=kp*e+ki* integral(e)
now how do i code this last part to find b at every iteration with the pi control active?
Best regards
n=50 with a step=1
the part of the code i'm interested in is this:
if w<wmax then e=0 %where e is the error i need to correct with pi controller
else
e=wwmax
Kp(1)=0.01882681
Ki(1)=0.008068634
bk=6.302336
Gk(i)=1/(1+b(i1)/bk)
Kp(i) =Kp(1)*Gk(i)
Ki(i) =Ki(1)*Gk(i)
so i'll have a :
b=kp*e+ki* integral(e)
now how do i code this last part to find b at every iteration with the pi control active?
Best regards

 Posts: 3857
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Question about PI control
Dear Vlad,
The DISCON.f90 file in the FAST archive contains the Fortran implementation of the control logic for the NREL 5MW baseline controller. You should be able to review that source file to see the exact coding details. Starting from the Fortran implementation, it should not be hard to recode the controller in MATLAB.
Best regards,
The DISCON.f90 file in the FAST archive contains the Fortran implementation of the control logic for the NREL 5MW baseline controller. You should be able to review that source file to see the exact coding details. Starting from the Fortran implementation, it should not be hard to recode the controller in MATLAB.
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
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

 Posts: 11
 Joined: Tue Feb 21, 2017 12:41 pm
 Organization: Middle East Technical University
 Location: Ankara/Turkey
Re: Question about PI control
Dear Jonkman,
In my wind turbine model prepared in matlab/simulink, i am trying to implemet Gain Scheduled PI Blade Pitch Controller which you applied before.
Following is what i wrote in a Simulink Block for the Gain Schudeled PI Control for Region 3, i.e blade pitch control in order to calculate your case first . I used the values of parameters according to your explanations as much as i understand fro 5MW turbine. But it seems giving little bit different results when compared to your plot, ThetaK, Kp and Ki. I do not know which value selected below causing this difference. I think, Omg0 should be in rpm not in radian. I=Irotor+Ngear^2*Igen. Damping is selected as 0.7. There seems no problem. I do not know little difference in Ki and Kp effects the results much either.
Omg0=12.1
I=(115926+(97^2)*534.116)
wn=0.6
dping=0.7
Ngear=97
Ps=25.52e+06
thetak=6.302336
GK=1/(1+theta/thetak)
payi=I*Omg0*(wn^2)
payp=2*I*Omg0*dping*wn
payda=(Ngear*Ps)
Ki=(payi/payda)*GK
Kp=(payp/payda)*GK
After seeing similar calculation results in your plot. I did the above same things for my case, NREL 5MW turbine with no tilt and precone angles. Followings are according to my case. I calculated a different pitch sensitivity, Ps using my model and ThetaK and Ki and Kp gains.
Omg0=12.1
I=(115926+(97^2)*534.116)
wn=0.6
dping=0.7
Ngear=97
Ps=22507000
thetak=5.1217
GK=1/(1+theta/thetak)
payi=I*Omg0*(wn^2)
payp=2*I*Omg0*dping*wn
payda=(Ngear*Ps)
Ki=(payi/payda)*GK
Kp=(payp/payda)*GK
In the attachement 1, the best fit line of blade pitch sensitivity is given. My results and your results are on the same plot. I assumed frozen wake assumption. But at the rated speed, 11,4 m/s, my pitch angle value which gives the rated power (5.296610 MW) is not zero degree unlike yours, it is 1.9455 degrees. But for now, i took the pitch angle as zero, which does not give the rated power. Is this angle always should be zero? or is this the pitch angle that gives the rated power at the rated wind speed? just like others giving the rated power? If so, they may change from a wind turbine simulation model to model.
Attachement 2 shows your GK and the PI gains and mine together on the same plot. I divided GK by 10 to let you see all results clearly on same plot.
My model does not include any blade actuator dynamics, no saturation, or no rate limiter etc. only Gain Scheduled PI Controller. The pitch angle, theta, to calculate gain correction factor, GK, is obtained by a lookup table which includes the wind speed information and their respective pitch angles giving the rated power.
The PI controlller output, which is the pitch angle, is in attachement 3. The power output of the turbine, on the other hand,is seen in attachement 4.
Last thing, i applied a PI control with no gain scheduling. It works but i realized that the pitch angles which gives the rated power is negative in all above rated wind speeds. I think, the pitch angles sould be larger and positive when the wind speed increase in region 3. When i check the model, in fact positive pitch angles also give the rated power but the controller gives the negative ones since negative pitch angles also give rated power between 6 and 8 degree pitch angles. Any idea? or do i have to limit something in model?
i would be very grateful if you could help me
Yours sincerely.
Mustafa SAHIN.
METU Aerospace Engineering Department
In my wind turbine model prepared in matlab/simulink, i am trying to implemet Gain Scheduled PI Blade Pitch Controller which you applied before.
Following is what i wrote in a Simulink Block for the Gain Schudeled PI Control for Region 3, i.e blade pitch control in order to calculate your case first . I used the values of parameters according to your explanations as much as i understand fro 5MW turbine. But it seems giving little bit different results when compared to your plot, ThetaK, Kp and Ki. I do not know which value selected below causing this difference. I think, Omg0 should be in rpm not in radian. I=Irotor+Ngear^2*Igen. Damping is selected as 0.7. There seems no problem. I do not know little difference in Ki and Kp effects the results much either.
Omg0=12.1
I=(115926+(97^2)*534.116)
wn=0.6
dping=0.7
Ngear=97
Ps=25.52e+06
thetak=6.302336
GK=1/(1+theta/thetak)
payi=I*Omg0*(wn^2)
payp=2*I*Omg0*dping*wn
payda=(Ngear*Ps)
Ki=(payi/payda)*GK
Kp=(payp/payda)*GK
After seeing similar calculation results in your plot. I did the above same things for my case, NREL 5MW turbine with no tilt and precone angles. Followings are according to my case. I calculated a different pitch sensitivity, Ps using my model and ThetaK and Ki and Kp gains.
Omg0=12.1
I=(115926+(97^2)*534.116)
wn=0.6
dping=0.7
Ngear=97
Ps=22507000
thetak=5.1217
GK=1/(1+theta/thetak)
payi=I*Omg0*(wn^2)
payp=2*I*Omg0*dping*wn
payda=(Ngear*Ps)
Ki=(payi/payda)*GK
Kp=(payp/payda)*GK
In the attachement 1, the best fit line of blade pitch sensitivity is given. My results and your results are on the same plot. I assumed frozen wake assumption. But at the rated speed, 11,4 m/s, my pitch angle value which gives the rated power (5.296610 MW) is not zero degree unlike yours, it is 1.9455 degrees. But for now, i took the pitch angle as zero, which does not give the rated power. Is this angle always should be zero? or is this the pitch angle that gives the rated power at the rated wind speed? just like others giving the rated power? If so, they may change from a wind turbine simulation model to model.
Attachement 2 shows your GK and the PI gains and mine together on the same plot. I divided GK by 10 to let you see all results clearly on same plot.
My model does not include any blade actuator dynamics, no saturation, or no rate limiter etc. only Gain Scheduled PI Controller. The pitch angle, theta, to calculate gain correction factor, GK, is obtained by a lookup table which includes the wind speed information and their respective pitch angles giving the rated power.
The PI controlller output, which is the pitch angle, is in attachement 3. The power output of the turbine, on the other hand,is seen in attachement 4.
Last thing, i applied a PI control with no gain scheduling. It works but i realized that the pitch angles which gives the rated power is negative in all above rated wind speeds. I think, the pitch angles sould be larger and positive when the wind speed increase in region 3. When i check the model, in fact positive pitch angles also give the rated power but the controller gives the negative ones since negative pitch angles also give rated power between 6 and 8 degree pitch angles. Any idea? or do i have to limit something in model?
i would be very grateful if you could help me
Yours sincerely.
Mustafa SAHIN.
METU Aerospace Engineering Department

 Posts: 11
 Joined: Tue Feb 21, 2017 12:41 pm
 Organization: Middle East Technical University
 Location: Ankara/Turkey
Re: Question about PI control
Attachements are given below.
Mustafa SAHIN
Mustafa SAHIN
 Attachments

 3pitch.png (7.67 KiB) Viewed 2201 times

 2Gains.jpg (29.34 KiB) Viewed 2201 times

 1Pitch sensitivity.jpg (39.55 KiB) Viewed 2201 times

 Posts: 11
 Joined: Tue Feb 21, 2017 12:41 pm
 Organization: Middle East Technical University
 Location: Ankara/Turkey
Re: Question about PI control
Last Attachement (4)
Mustafa SAHIN
Mustafa SAHIN
 Attachments

 4power.png (7.37 KiB) Viewed 2200 times

 Posts: 11
 Joined: Tue Feb 21, 2017 12:41 pm
 Organization: Middle East Technical University
 Location: Ankara/Turkey
Re: Question about PI control
Attachement 2 with explantions
 Attachments

 2gains.jpg (60.2 KiB) Viewed 2199 times

 Posts: 3857
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Question about PI control
Dear Mustafa,
I can't follow everything you are saying, but here are a few comments:
Best regards,
I can't follow everything you are saying, but here are a few comments:
 The rated rotor speed (Omg0) should have the units of rad/s, not rpm i.e. you should change 12.1 to 1.2671 rad/s.
 You are using the hub inertia of the NREL 5MW turbine (115926 kg*m^2) instead of the rotor inertia (38759228 kg*m^2).
 The bladepitch angle at rated need not be zero degrees, but it is commonly close to zero.
 Pitching the blades negatively will result in pitchtostall control rather than pitchtofeather control.
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
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

 Posts: 11
 Joined: Tue Feb 21, 2017 12:41 pm
 Organization: Middle East Technical University
 Location: Ankara/Turkey
Re: Question about PI control
Dear Jonkman,
Thank for your quick reply. I will apply what you have said. when i have further questions, i write them here on Forum.
Yours sincerely.
Mustafa SAHIN
Thank for your quick reply. I will apply what you have said. when i have further questions, i write them here on Forum.
Yours sincerely.
Mustafa SAHIN

 Posts: 11
 Joined: Tue Feb 21, 2017 12:41 pm
 Organization: Middle East Technical University
 Location: Ankara/Turkey
Re: Question about PI control
Dear Jonkman,
First of all, thank you for sharing your valuable ideas, experinces with us. I have applied the Gain Scheduling method to the turbine's above rated region. It seems working! but not as good as i expect. When i check simulation results, It takes long time (150200 second) for the signals to converge the steady state with wn=0.6 damping ratio=0.7, no rate limiter, no limitation like 090 degrees in blade pitch etc. Usually in how much second the system will convenge to steady state? Since the 5MW is quite large, it should take long time?
Another thing is about the value of rotor inertia. You said it is 38759228 kg*m^2. I could not find this value in the article ''Definition of a 5MW Reference Wind Turbine for Offshore System Development,'' which is published in February 2009. Could you please recommend me a published documentation for that.
For the application of my Gain scheduling, although my pitch angle at 11.4 m/s rated wind speed that gives the rated power is different than zero, 1.9455 which is not even so close to zero, i obtained the pitch sensitivity at zero pitch angle from the best fit line to my data. i obtained thetak which is the blade pitch angle corresponding to two times the pitch sensitivity at zero blade pitch angle. When i do it as explained above, i get similar results to yours in terms of pitch sensitivity and thetak. My pitch sensitivity and thetak are 24164000 watt/rad and 5.6232 degrees. Yours are, respectively, 25520000 and 6.302336 degrees. This may be caused due to no precone and tilt or some adjustments in my model compared to yours.
if i calculate the pitch sensitivity from the best fit line corresponding to 1.9544 blade pitch angle (instead of zero) which gives the rated power at 11.4 m/s and then obtain thetak as the two times of pitch sensitiviy at 1.9544. The pitch sensitivity and thetak are quite different than yours, respectively as 32524000 watt/rad and 9.5141 degrees. I think this is wrong! (?) Therefore, while calculating pitch sensitivity, i should take the pitch sentivity at zero degree blade pitch whatever the blade pitch angle is calculated from the simulation model at rated wind speed? Because my all results are almost same as yours but not the blade pitch angle at rated wind speed.
Thank you.
Yours sincerely.
Mustafa SAHIN
First of all, thank you for sharing your valuable ideas, experinces with us. I have applied the Gain Scheduling method to the turbine's above rated region. It seems working! but not as good as i expect. When i check simulation results, It takes long time (150200 second) for the signals to converge the steady state with wn=0.6 damping ratio=0.7, no rate limiter, no limitation like 090 degrees in blade pitch etc. Usually in how much second the system will convenge to steady state? Since the 5MW is quite large, it should take long time?
Another thing is about the value of rotor inertia. You said it is 38759228 kg*m^2. I could not find this value in the article ''Definition of a 5MW Reference Wind Turbine for Offshore System Development,'' which is published in February 2009. Could you please recommend me a published documentation for that.
For the application of my Gain scheduling, although my pitch angle at 11.4 m/s rated wind speed that gives the rated power is different than zero, 1.9455 which is not even so close to zero, i obtained the pitch sensitivity at zero pitch angle from the best fit line to my data. i obtained thetak which is the blade pitch angle corresponding to two times the pitch sensitivity at zero blade pitch angle. When i do it as explained above, i get similar results to yours in terms of pitch sensitivity and thetak. My pitch sensitivity and thetak are 24164000 watt/rad and 5.6232 degrees. Yours are, respectively, 25520000 and 6.302336 degrees. This may be caused due to no precone and tilt or some adjustments in my model compared to yours.
if i calculate the pitch sensitivity from the best fit line corresponding to 1.9544 blade pitch angle (instead of zero) which gives the rated power at 11.4 m/s and then obtain thetak as the two times of pitch sensitiviy at 1.9544. The pitch sensitivity and thetak are quite different than yours, respectively as 32524000 watt/rad and 9.5141 degrees. I think this is wrong! (?) Therefore, while calculating pitch sensitivity, i should take the pitch sentivity at zero degree blade pitch whatever the blade pitch angle is calculated from the simulation model at rated wind speed? Because my all results are almost same as yours but not the blade pitch angle at rated wind speed.
Thank you.
Yours sincerely.
Mustafa SAHIN

 Posts: 3857
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Question about PI control
Dear Mustafa,
How long the solution will take to converge to steady state depends on the turbine you modeling, the degrees of freedom enabled, the initial conditions you've set, and the wind conditions you are operating with. I can't really comment in general.
The rotor inertia is not documented in the NREL 5MW baseline turbine specifications report, but it can be calculated by running the FAST model (the rotor inertia is written to the FAST summary (*.fsm) file).
Regarding the gain scheduling, all you really need is the slope of the curve of ∂P/∂θ versus θ. Selecting θk as the value of the pitch that gives ∂P/∂θ = 2*∂P/∂θ(θ=0) is just a convenient way of calculating the slope. But I don't really know anything about your turbine, so, can't really comment on whether the values you are reporting make sense.
Best regards,
How long the solution will take to converge to steady state depends on the turbine you modeling, the degrees of freedom enabled, the initial conditions you've set, and the wind conditions you are operating with. I can't really comment in general.
The rotor inertia is not documented in the NREL 5MW baseline turbine specifications report, but it can be calculated by running the FAST model (the rotor inertia is written to the FAST summary (*.fsm) file).
Regarding the gain scheduling, all you really need is the slope of the curve of ∂P/∂θ versus θ. Selecting θk as the value of the pitch that gives ∂P/∂θ = 2*∂P/∂θ(θ=0) is just a convenient way of calculating the slope. But I don't really know anything about your turbine, so, can't really comment on whether the values you are reporting make sense.
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
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

 Posts: 43
 Joined: Tue Nov 14, 2017 6:50 am
 Organization: ECN
 Location: france
Re: Question about PI control
Dear Jason,
I try to use SIMULINK to repeat to the Baseline PI control, but the result is very different from the .DISCON.dll you proposed.
Here is my SIMULINK model, could you please tell me where I was wrong? and I'm quiet cannot understand this line IntSpdErr = MIN( MAX( IntSpdErr, PC_MinPit/( GK*PC_KI ) ), PC_MaxPit/( GK*PC_KI )) . I don't know its effect on the system.
Thank you in advance.
Best regards,
Cheng
I try to use SIMULINK to repeat to the Baseline PI control, but the result is very different from the .DISCON.dll you proposed.
Here is my SIMULINK model, could you please tell me where I was wrong? and I'm quiet cannot understand this line IntSpdErr = MIN( MAX( IntSpdErr, PC_MinPit/( GK*PC_KI ) ), PC_MaxPit/( GK*PC_KI )) . I don't know its effect on the system.
Thank you in advance.
Best regards,
Cheng
 Attachments

 Capture.JPG (91.18 KiB) Viewed 1662 times
Who is online
Users browsing this forum: No registered users and 1 guest