📄 fm_cswt.m
字号:
function fm_cswt(flag)
% FM_CSWT define Constant Speed Wind Turbine
%
% FM_CSWT(FLAG)
% FLAG = 0 initializations
% FLAG = 1 algebraic equations
% FLAG = 2 algebraic Jacobians
% FLAG = 3 differential equations
% FLAG = 4 state matrix
%
%see also FM_WIND
%
%Author: Federico Milano
%Date: 30-Nov-2003
%Version: 1.0.0
%
%E-mail: fmilano@thunderbox.uwaterloo.ca
%Web-site: http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano
global Bus Cswt Wind DAE Settings Varname PV SW
Wn = Settings.rad;
omega_wr = DAE.x(Cswt.omega_wr);
omega_m = DAE.x(Cswt.omega_m);
gamma = DAE.x(Cswt.gamma);
e1r = DAE.x(Cswt.e1r);
e1m = DAE.x(Cswt.e1m);
vw = DAE.x(Wind.vw(Cswt.wind));
rho = Wind.con(Cswt.wind,3);
V = DAE.V(Cswt.bus);
t = DAE.a(Cswt.bus);
st = sin(t);
ct = cos(t);
Vr = -V.*st;
Vm = V.*ct;
r1 = Cswt.con(:,6);
rr = Cswt.con(:,8);
x0 = Cswt.dat(:,1);
x1 = Cswt.dat(:,2);
T10 = Cswt.dat(:,3);
i2Hwr = Cswt.dat(:,4);
i2Hm = Cswt.dat(:,5);
Ks = Cswt.con(:,13);
R = Cswt.dat(:,6);
A = Cswt.dat(:,7);
B = Cswt.dat(:,8);
a = r1.^2+x1.^2;
a13 = r1./a;
a23 = x1./a;
a33 = x0-x1;
Im = -a23.*(e1r-Vr)+a13.*(e1m-Vm);
Ir = a13.*(e1r-Vr)+a23.*(e1m-Vm);
switch flag
case 0
check = 1;
% Constants
% x0
Cswt.dat(:,1) = Cswt.con(:,7) + Cswt.con(:,10);
% x'
Cswt.dat(:,2) = Cswt.con(:,7) + Cswt.con(:,10).*Cswt.con(:,9)./ ...
(Cswt.con(:,10)+Cswt.con(:,9));
% T'0
Cswt.dat(:,3) = (Cswt.con(:,10)+Cswt.con(:,9))./Cswt.con(:,8)/Wn;
% 1/(2*Hwr)
Cswt.dat(:,4) = 1./(2*Cswt.con(:,11));
% 1/(2*Hm)
Cswt.dat(:,5) = 1./(2*Cswt.con(:,12));
% 4*R*pi*f/p
Cswt.dat(:,6) = 2*Wn*Cswt.con(:,14)./Cswt.con(:,15).*Cswt.con(:,17);
% A
Cswt.dat(:,7) = pi*Cswt.con(:,14).*Cswt.con(:,14);
% Initialization of state variables
Pc = Bus.Pg(Cswt.bus);
Qc = Bus.Qg(Cswt.bus);
Vc = DAE.V(Cswt.bus);
ac = DAE.a(Cswt.bus);
vr = -Vc.*sin(ac);
vm = Vc.*cos(ac);
for i = 1:Cswt.n
% parameters
Vr = vr(i);
Vm = vm(i);
Rs = r1(i);
X0 = Cswt.dat(i,1);
X1 = Cswt.dat(i,2);
T10 = Cswt.dat(i,3);
Pg = Pc(i);
Qg = Qc(i);
eqn = ones(5,1);
inc = ones(5,1);
jac = zeros(5,5);
jac(1,1) = Vr;
jac(1,2) = Vm;
jac(2,1) = -Rs;
jac(2,2) = X1;
jac(2,3) = 1;
jac(3,1) = -X1;
jac(3,2) = -Rs;
jac(3,4) = 1;
jac(4,2) = (X0-X1)/T10;
jac(4,3) = -1/T10;
jac(5,1) = -(X0-X1)/T10;
jac(5,4) = -1/T10;
% variables: x(1) = ir
% x(2) = im
% x(3) = e'r
% x(4) = e'm
% x(5) = sigma
% first guess
x = jac([1:5],[1:4])\[Pg;Vr;Vm;0;0];
x(5) = 0;
iter = 0;
while max(abs(inc)) > 1e-5
if iter > 20
fm_disp(['Initialization of constant speed wind turbine #',...
num2str(i),' failed (convergence problem).'])
break
end
eqn(1) = Vr*x(1) + Vm*x(2) - Pg;
eqn(2) = x(3) - Vr - Rs*x(1) + X1*x(2);
eqn(3) = x(4) - Vm - Rs*x(2) - X1*x(1);
eqn(4) = Wn*x(5)*x(4) - (x(3)-(X0-X1)*x(2))/T10;
eqn(5) = -Wn*x(5)*x(3) - (x(4)+(X0-X1)*x(1))/T10;
jac(4,4) = Wn*x(5);
jac(4,5) = Wn*x(4);
jac(5,3) = -Wn*x(5);
jac(5,5) = -Wn*x(3);
inc = -jac\eqn;
x = x + inc;
iter = iter + 1;
end
Qg = Vm*x(1)-Vr*x(2);
Te = x(3)*x(1)+x(4)*x(2);
% shunt capacitor
Cswt.dat(i,8) = (Qc(i)-Qg)/V(i)/V(i);
% wind turbine state variables
DAE.x(Cswt.e1r(i)) = x(3);
DAE.x(Cswt.e1m(i)) = x(4);
DAE.x(Cswt.omega_m(i)) = 1-x(5);
DAE.x(Cswt.gamma(i)) = Te/Cswt.con(i,13);
DAE.x(Cswt.omega_wr(i)) = DAE.x(Cswt.omega_m(i));
% wind power
Pw = Te*DAE.x(Cswt.omega_m(i));
if Pg < 0
fm_disp([' * * Turbine power is negative at bus <',Varname.bus{Cswt.bus(i)},'>.'])
fm_disp([' Wind speed <',num2str(Cswt.wind(i)),'> cannot be initilized.'])
DAE.x(Wind.vw(Cswt.wind(i))) = 1;
continue
end
% wind speed
iter = 0;
incvw = 1;
eqnvw = 1;
R = Cswt.dat(i,6);
A = Cswt.dat(i,7);
vw = 0.9*Wind.con(Cswt.wind(i),2);
while abs(eqnvw) > 1e-8
if iter > 20
wspeed = num2str(Cswt.wind(i));
fm_disp([' * * Initialization of wind speed <', ...
wspeed, '> failed (convergence problem).'])
fm_disp([' Tip: Try increasing the nominal wind speed <',wspeed,'>.'])
check = 0;
break
end
eqnvw = windpower(rho(i),vw, ...
A,R,DAE.x(Cswt.omega_m(i)),1)-Pw*Settings.mva*1e6;
jacvw = windpower(rho(i),vw, ...
A,R,DAE.x(Cswt.omega_m(i)),2);
incvw = -eqnvw/jacvw(2);
%disp([eqnvw+Pw*Cswt.con(i,3),incvw])
vw = vw + incvw;
%disp([vw,iter]);
iter = iter + 1;
end
DAE.x(Wind.vw(Cswt.wind(i))) = vw/Wind.con(Cswt.wind(i),2);
% find & delete static generators
if ~fm_rmgen(Cswt.bus(i)), check = 0; end
end
% random initial phase angle for shadow tower effect
%Cswt.dat(:,9) = 2*pi*rand(Cswt.n,1);
if ~check
fm_disp('Constant speed wind turbine cannot be properly initialized.')
else
fm_disp(['Constant speed wind turbines initialized.'])
end
case 1
DAE.gp = DAE.gp - sparse(Cswt.bus,1,Vr.*Ir+Vm.*Im,Bus.n,1);
DAE.gq = DAE.gq - sparse(Cswt.bus,1,Vm.*Ir-Vr.*Im+B.*V.*V,Bus.n,1);
case 2
Pv = -Ir.*st + Im.*ct;
Pt = -Ir.*Vm + Im.*Vr;
Qv = Ir.*ct + Im.*st;
Qt = Ir.*Vr + Im.*Vm;
IrV = a13.*st-a23.*ct;
ImV = -a23.*st-a13.*ct;
Irt = a13.*Vm-a23.*Vr;
Imt = -a23.*Vm-a13.*Vr;
DAE.J11 = DAE.J11 - sparse(Cswt.bus,Cswt.bus,Pt+Vr.*Irt+Vm.*Imt,Bus.n,Bus.n);
DAE.J12 = DAE.J12 - sparse(Cswt.bus,Cswt.bus,Pv+Vr.*IrV+Vm.*ImV,Bus.n,Bus.n);
DAE.J21 = DAE.J21 - sparse(Cswt.bus,Cswt.bus,Qt+Vm.*Irt-Vr.*Imt,Bus.n,Bus.n);
DAE.J22 = DAE.J22 - sparse(Cswt.bus,Cswt.bus,Qv+Vm.*IrV-Vr.*ImV+2*B.*V,Bus.n,Bus.n);
case 3
slip = 1-omega_m;
% periodic torque pulsation due to "tower shadow" phenomenon
omega_r = 2*Wn*omega_wr.*Cswt.con(:,17)./Cswt.con(:,15)./Cswt.con(:,16);
shadow = 0.025*sin(max(DAE.t,0)*omega_r);
Twr = windpower(rho,vw.*Wind.con(Cswt.wind,2), ...
A,R,omega_wr,1)./omega_wr./Settings.mva/1e6;
DAE.f(Cswt.omega_wr) = (Twr.*(1+shadow)-Ks.*gamma).*i2Hwr;
DAE.f(Cswt.omega_m) = (Ks.*gamma-e1r.*Ir-e1m.*Im).*i2Hm;
DAE.f(Cswt.gamma) = Wn*(omega_wr-omega_m);
DAE.f(Cswt.e1r) = Wn*slip.*e1m-(e1r-a33.*Im)./T10;
DAE.f(Cswt.e1m) = -Wn*slip.*e1r-(e1m+a33.*Ir)./T10;
case 4
dPwdx = windpower(rho,vw.*Wind.con(Cswt.wind,2), ...
A,R,omega_wr,2)./Settings.mva/1e6;
Twr = windpower(rho,vw.*Wind.con(Cswt.wind,2), ...
A,R,omega_wr,1)./omega_wr./Settings.mva/1e6;
slip = 1-omega_m;
IrV = a13.*st-a23.*ct;
ImV = -a23.*st-a13.*ct;
Irt = a13.*Vm-a23.*Vr;
Imt = -a23.*Vm-a13.*Vr;
DAE.Fx = DAE.Fx + sparse(Cswt.omega_wr,Cswt.gamma,-Ks.*i2Hwr,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.omega_m,Cswt.gamma,Ks.*i2Hm,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.gamma,Cswt.omega_wr,Wn,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.gamma,Cswt.omega_m,-Wn,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.omega_wr,Cswt.omega_wr,(dPwdx(:,1)-Twr).*i2Hwr./omega_wr,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.omega_wr,Wind.vw(Cswt.wind),dPwdx(:,2).*i2Hwr./omega_wr,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.omega_m,Cswt.e1r,(-Ir-e1r.*a13+e1m.*a23).*i2Hm,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.omega_m,Cswt.e1m,(-Im-e1r.*a23-e1m.*a13).*i2Hm,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.e1r,Cswt.omega_m,-Wn.*e1m,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.e1r,Cswt.e1r,-(1+a33.*a23)./T10,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.e1r,Cswt.e1m,Wn.*slip+a33.*a13./T10,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.e1m,Cswt.omega_m,Wn.*e1r,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.e1m,Cswt.e1r,-Wn.*slip-a33.*a13./T10,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Cswt.e1m,Cswt.e1m,-(1+a33.*a23)./T10,DAE.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(Cswt.bus,Cswt.e1r, a13.*Vr-a23.*Vm,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(Cswt.bus,Cswt.e1m, a23.*Vr+a13.*Vm,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(Cswt.bus+Bus.n,Cswt.e1r, a23.*Vr+a13.*Vm,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(Cswt.bus+Bus.n,Cswt.e1m,-a13.*Vr+a23.*Vm,2*Bus.n,DAE.n);
DAE.Fy = DAE.Fy + sparse(Cswt.omega_m,Cswt.bus,(-e1r.*Irt-e1m.*Imt).*i2Hm,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Cswt.e1r,Cswt.bus,a33.*Imt./T10,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Cswt.e1m,Cswt.bus,-a33.*Irt./T10,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Cswt.omega_m,Cswt.bus+Bus.n,(-e1r.*IrV-e1m.*ImV).*i2Hm,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Cswt.e1r,Cswt.bus+Bus.n,a33.*ImV./T10,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Cswt.e1m,Cswt.bus+Bus.n,-a33.*IrV./T10,DAE.n,2*Bus.n);
end
%--------------------------------------------------
function output = windpower(rho,vw,Ar,R,omega,type)
%--------------------------------------------------
lambda = omega.*R./vw;
lambdai = 1./(1./lambda+0.002);
switch type
case 1 % Pw
cp = 0.44*(125./lambdai-6.94).*exp(-16.5./lambdai);
output = 0.5*rho.*cp.*Ar.*vw.^3;
case 2
output = zeros(length(omega),2);
a1 = exp(-16.5./lambda-0.0330);
a2 = 125./lambda-6.690;
% d Pw / d omega
output(:,1) = rho.*(-27.5+3.63*a2).*a1.*Ar.*vw.^4./omega./omega./R;
% d Pw / d vw
output(:,2) = rho.*((27.5-3.63*a2)./lambda+0.66*a2).*a1.*Ar.*vw.^2;
% > diff(Pw,vw);
% 16.5 vw 3
%27.50000000 rho exp(- ------- - 0.0330) Ar vw
% omega R
%----------------------------------------------
% omega R
%
% /125 vw \ 16.5 vw 3
% 3.630000000 rho |------- - 6.690| exp(- ------- - 0.0330) Ar vw
% \omega R / omega R
% - ----------------------------------------------------------------
% omega R
%
% /125 vw \ 16.5 vw 2
% + 0.6600000000 rho |------- - 6.690| exp(- ------- - 0.0330) Ar vw
% \omega R / omega R
%
%> diff(Pw,omega);
% 4 16.5 vw
% 27.50000000 rho vw exp(- ------- - 0.0330) Ar
% omega R
%- ----------------------------------------------
% 2
% omega R
%
% /125 vw \ 4 16.5 vw
% 3.630000000 rho |------- - 6.690| vw exp(- ------- - 0.0330) Ar
% \omega R / omega R
% + ----------------------------------------------------------------
% 2
% omega R
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -