%%
% *By: José Salvador Sánchez Osorio*
%%
% 
% Introduction to optimization E602
% 
%%
% 
% This project is based on the paper MPC for Wind Power
% Gradients-Utilizing Forecasts, Rotor Inertia, and Central Energy Storage
% by T. Hovgaard, L. Larsen, J. Jorgensen and S. Boyd
% 

%%
% 
% The paper bases their results from a NREL 5-MW Wind Turbine, stated on 
% the Technical Report of the National Renewable Energy Laboratory, 
% from which the nominal information to model the dynamics of the system 
% can be obtained.
% According to the paper simulations for the specs of the Wind Turbine the
% peak power coefficient is 0.482 when lambda=7.55 and Betha=0°
%   

 

 


%%
% *Rated information*
wgratedmax=1173.7/60;     %[rps] rated generator max speed
Prated=5000000;     %[W] rated electric power
Tgrated=43094.55;      %[N.m] rated generator torque
wgratedmin=670/60;        %[rps] min generator speed
Betharef=0;            %[degrees] Pitch angle for rated speed
%%
% *Rotor information*
mr=110000;          %[kg]  rotor mass
wr=12.1/60;             %[rps] rotor speed
dr=125.88;           %[m]   rotor diameter
R=dr/2;              %[m]   rotor radius
A=12445.3;           %[m^2] swept area
Ir=11776047;         %[kg.m^2] Rotor inertia
Tg=Tgrated;
%%
% *Generator information*

wg=1173.7/60;           %[rpms]  generator speed
N=97;                %gearbox ratio
ng=0.944;            %generator efficiency
Ig=534.116;          %[kg.m^2] Generator Inertia
Tgref=28116.2;       %[N.m]    Torque 
Taog=0.6;
%%
% *-------Parameters----------*
% 
% In the paper the following parameters are used to obtain the desired
% results
Np=300;             %5 minute horizon
rho=8500;           %[kg/m^3] effective density of steel
t=1:800/Np:800;      %[s]      time 
%%
% To follow the method used in the paper the first value of the speed is the
% rated value 11.4, after that they use a value of 10 m/s for the first 400
% seconds, then it linearly descend down to 8 m/s in the interval from 400s
% to 430s, in the paper there are 4 different tests with different values
% for the speed but for the sake of simplicity only one of the tests will be
% considered for this project

for k=1:Np      
    if t(k)<400
        v(k)=10;
    end
    if t(k)>=400
        v(k)=v(k-1)-2/11;
    end
    if t(k)>430
        v(k)=8;
    end
end
Q=27.9;             %[pu]
gradpu=5.10e-4;     %[pu/s]
Ts=1;
Betha=zeros(1,Np);
%%
% _Betha values depend on wind speed_
for k=1:Np
    if v(k)==12
        Betha(k)=3.83;
    end
    if v(k)==13
        Betha(k)=6.6;
    end
    if v(k)==14
        Betha(k)=8.7;
    end
    if v(k)==15
        Betha(k)=10.45;
    end
end
       
v(1)=11.4;          %m/s wind speed rated value
%%
%
% The values of the coefficient of power [Cp] depend on lambda and Betha,
% the authors of the paper provide the following figure as source for their
% calculations of the Cp, but it is not clear enough to get the values and
% their source material does not state the values for that coefficient so
% the values I am using are taken from [2]
%

%%
% 
% <<Img1.png>>
% 


Tg=1/Taog*(Tgref-Tg);           %generator speed
lambda=R*wr./v;                 %Tip speed Ratio
if lambda>=7.55
    Cp=0.482
end
if (lambda<7.55)&&(lambda>=6)
    Cp=0.35
end
if (lambda<6)&&(lambda>=5)
    Cp=0.3
end
if (lambda<5)&&(lambda>=4)
    Cp=0.25
end    
if (lambda<4)&&(lambda>=3)
    Cp=0.1
end
if (lambda<3)
    Cp=0.02
end

cvx_begin

variable t


% *Wind Turbine Model*
wg=1/(Ig+Ir/(N^2))*(Tr/N-Tg);
Tg=(1/Taog)*(Tgref-Tg);
Tr=1/2*rho*A*Cp*(v^3)/wr;       %rotor torque
Pg=ng*Tg*wg;                    %generator power
Q(t+Ts)=Q(t)+c(t)*Ts;                       %Energy storage
Pgrid=Pg-c(t);                     %Power supplied
%%
% The functions that represent the energy storage Q and the charge rate c
% are not stated in source paper


% *Cost*
mu=Pgrid;
E=int(Pgrid,t,0,800);           %Supplied energy
S=max(mu,T0,Tf)*cstor;          %Storage cost
cstor=1;                        %Tunning parameter whose value is not stated on the source material
%%
% mu is defined as a needed storage capacity, this value depends on
% the kind of batteries used and the value is not stated in the paper, but
% depends on the power supplied so mu is assumed as Pgrid

maximize E-S
subject to 
Tg;
Betha;
c;


% *Constrains*

if v>=25
   Pw=0;
end
Pmin<=Pg<=Prated;
Pg<=Pw;
0<=Betha<=90;
-8<=(Betha-Betha)/Ts<=8;
c<=Pg;
0<=Q<=mu*Prated;

wgratedmin<=wg<=(1+wos)*wgratedmax;
-gradpu<=(Pgrid-Pgrid)/(Prated*Ts)<=gradpu;

cvx_end

%%
% *Results*
%%
% 
% Due to not having the functions that describe Q and c the program cannot
% run.
% 

%%
% *References*
%%
% 
% # [1]Hovgaard, T. G., Larsen, L. F., Jorgensen, J. B., & Boyd, S. (2013). MPC for Wind Power Gradients-Utilizing Forecasts, Rotor Inertia, and Central Storage. European Control Conference, 4071-4076. Retrieved April 26, 2019.
% # [2]J. M. Jonkman, S. Butterfield, W. Musial, and G. Scott, Definition of a 5-MW reference wind turbine for offshore system development. National Renewable Energy Laboratory, Feb 2009.
% # [3]Royal Academy of Engineering. (n.d.). Wind Turbine Power Calculations. Retrieved April 26, 2019, from https://www.raeng.org.uk/publications/other/23-wind-turbine
% 
