📄 fm_ddsg.m
字号:
function fm_ddsg(flag)
% FM_DDSG define a direct drive synchronous generator wind turbine
%
% FM_DDSG(FLAG)
% FLAG = 0 initializations
% FLAG = 1 algebraic equations
% FLAG = 2 algebraic Jacobians
% FLAG = 3 differential equations
% FLAG = 4 state matrix
% FLAG = 5 non-windup limiters
%
%see also FM_WIND, FM_CSWT and FM_DFIG
%
%Author: Federico Milano
%Date: 12-May-2004
%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 Ddsg Wind DAE Settings Varname PV SW
omega_m = DAE.x(Ddsg.omega_m);
theta_p = DAE.x(Ddsg.theta_p);
ids = DAE.x(Ddsg.ids);
iqs = DAE.x(Ddsg.iqs);
idc = DAE.x(Ddsg.idc);
vw = DAE.x(Wind.vw(Ddsg.wind));
rho = Wind.con(Ddsg.wind,3);
V = DAE.V(Ddsg.bus);
th = DAE.a(Ddsg.bus);
st = sin(th);
ct = cos(th);
rs = Ddsg.con(:,6);
xd = Ddsg.con(:,7);
xq = Ddsg.con(:,8);
psip = Ddsg.con(:,9);
Hm = Ddsg.con(:,10);
Kp = Ddsg.con(:,11);
Tp = Ddsg.con(:,12);
Kv = Ddsg.con(:,13);
Tv = Ddsg.con(:,14);
Tep = Ddsg.con(:,15);
Teq = Ddsg.con(:,16);
R = Ddsg.dat(:,4);
A = Ddsg.dat(:,5);
Vref = Ddsg.dat(:,1);
Pref = Ddsg.dat(:,6);
Qref = Ddsg.dat(:,7);
ids_max = -Ddsg.con(:,24);
ids_min = -Ddsg.con(:,23);
idc_max = -Ddsg.con(:,24);
idc_min = -Ddsg.con(:,23);
iqs_max = -Ddsg.con(:,22);
iqs_min = -Ddsg.con(:,21);
switch flag
case 0
check = 1;
%check time constants
idx = find(Hm == 0);
if idx
ddsgwarn(idx, 'Inertia Hm cannot be zero. Hm = 2.5 s will be used.'),
end
Ddsg.con(idx,5) = 2.5;
idx = find(Tp == 0);
if idx
ddsgwarn(idx, 'Time constant Tp cannot be zero. Tp = 0.001 s will be used.'),
end
Ddsg.con(idx,6) = 0.001;
idx = find(Tv == 0);
if idx
ddsgwarn(idx, 'Time constant Tv cannot be zero. Tv = 0.001 s will be used.'),
end
Ddsg.con(idx,9) = 0.001;
idx = find(Teq == 0);
if idx
ddsgwarn(idx, 'Time constant Teq cannot be zero. Teq = 0.001 s will be used.'),
end
Ddsg.con(idx,15) = 0.001;
idx = find(Tep == 0);
if idx
ddsgwarn(idx, 'Time constant Tep cannot be zero. Tep = 0.001 s will be used.'),
end
Ddsg.con(idx,17) = 0.001;
% Constants
% etaGB*4*R*pi*f/p
Ddsg.dat(:,4) = 2*Settings.rad*Ddsg.con(:,17)./ ...
Ddsg.con(:,18).*Ddsg.con(:,20);
% A
Ddsg.dat(:,5) = pi*Ddsg.con(:,17).*Ddsg.con(:,17);
% Initialization of state variables
Pc = Bus.Pg(Ddsg.bus);
Qc = Bus.Qg(Ddsg.bus);
Vc = DAE.V(Ddsg.bus);
ac = DAE.a(Ddsg.bus);
vdc = -Vc.*sin(ac);
vqc = Vc.*cos(ac);
for i = 1:Ddsg.n
% idc
DAE.x(Ddsg.idc(i)) = cos(ac(i))*(Qc(i)-Pc(i)*tan(ac(i)))/Vc(i);
% vref
Ddsg.dat(i,1) = DAE.x(Ddsg.idc(i))/Kv(i)+Vc(i);
jac = zeros(4,4);
jac(3,1) = -1;
jac(4,2) = 1;
jac(3,3) = -rs(i);
jac(3,4) = xq(i);
jac(4,3) = xd(i);
jac(4,4) = rs(i);
x = zeros(4,1);
x(2) = 1;
x(3) = -1/xd(i)-psip(i);
iter = 0;
inc = 1;
while max(abs(inc)) > 1e-8
if iter > 20
fm_disp(['Initialization of direct drive syn. gen. #', ...
num2str(i),' failed.'])
check = 0;
break
end
eqn(1,1) = x(1)*x(3)+x(2)*x(4)-Pc(i);
eqn(2,1) = x(2)*x(3)-x(1)*x(4);
eqn(3,1) = -x(1)-rs(i)*x(3)+xq(i)*x(4);
eqn(4,1) = x(2)+rs(i)*x(4)+xd(i)*x(3)-psip(i);
jac(1,1) = x(3);
jac(1,2) = x(4);
jac(1,3) = x(1);
jac(1,4) = x(2);
jac(2,1) = -x(4);
jac(2,2) = x(3);
jac(2,3) = x(2);
jac(2,4) = -x(1);
inc = -jac\eqn;
x = x + inc;
%disp(x')
iter = iter + 1;
end
vds = x(1);
vqs = x(2);
ids = x(3);
iqs = x(4);
omega_m = 1;
% theta_p
theta_p = Kp(i)*round(1000*(omega_m-1))/1000;
theta_p = max(theta_p,0);
% wind turbine state variables and constants
DAE.x(Ddsg.ids(i)) = ids; % ids
Ddsg.dat(i,2) = ids; % ids_ref
DAE.x(Ddsg.iqs(i)) = iqs; % iqs
Ddsg.dat(i,3) = iqs; % iqs_ref
DAE.x(Ddsg.omega_m(i)) = omega_m;
DAE.x(Ddsg.theta_p(i)) = theta_p;
% electrical torque
Tel = ((xq(i)-xd(i))*ids+psip(i))*iqs;
% Ps_ref
Ddsg.dat(i,6) = iqs*(psip(i)-xd(i)*ids);
% Qs_ref
Ddsg.dat(i,7) = xd(i)*((psip(i)/xd(i))^2-(ids-psip(i)/xd(i))^2);
% wind power [MW]
Pw = Tel*omega_m*Settings.mva*1e6;
if Pc(i) < 0
fm_disp([' * * Turbine power is negative at bus <',Varname.bus{Ddsg.bus(i)},'>.'])
fm_disp([' Wind speed <',num2str(Ddsg.wind(i)),'> cannot be initilized.'])
DAE.x(Wind.vw(Ddsg.wind(i))) = 1;
continue
end
% wind speed
iter = 0;
incvw = 1;
eqnvw = 1;
R = Ddsg.dat(i,4);
A = Ddsg.dat(i,5);
% initial guess for wind speed
vw = 0.9*Wind.con(Ddsg.wind(i),2);
while abs(eqnvw) > 1e-7
if iter > 50
wspeed = num2str(Ddsg.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,omega_m,theta_p,1)-Pw;
jacvw = windpower(rho(i),vw,A,R,omega_m,theta_p,2);
incvw = -eqnvw/jacvw(2);
vw = vw + incvw;
iter = iter + 1;
end
% average initial wind speed [p.u.]
DAE.x(Wind.vw(Ddsg.wind(i))) = vw/Wind.con(Ddsg.wind(i),2);
% find & delete static generators
if ~fm_rmgen(Ddsg.bus(i)), check = 0; end
end
if ~check
fm_disp(['Direct drive synchronous generator cannot be properly ' ...
'initialized.'])
else
fm_disp(['Direct drive synchronous generators initialized.'])
end
case 1
vds = -rs.*ids+omega_m.*xq.*iqs;
vqs = -rs.*iqs-omega_m.*(xd.*ids-psip);
ps = vds.*ids+vqs.*iqs;
DAE.gp = DAE.gp - sparse(Ddsg.bus,1,ps,Bus.n,1);
DAE.gq = DAE.gq - sparse(Ddsg.bus,1,V.*idc./cos(th)+tan(th).*ps,Bus.n,1);
case 2
c1 = cos(th);
s1 = sin(th);
t1 = s1./c1;
vds = -rs.*ids+omega_m.*xq.*iqs;
DAE.J22 = DAE.J22 - sparse(Ddsg.bus,Ddsg.bus,idc./c1,Bus.n,Bus.n);
DAE.J21 = DAE.J21 - sparse(Ddsg.bus,Ddsg.bus,t1.*V.*idc./c1+(1+t1.^2).*vds.*ids,Bus.n,Bus.n);
case 3
% wind speed in m/s
% -----------------
Vw = vw.*Wind.con(Ddsg.wind,2);
% mechanical torque
% -----------------
Pw = windpower(rho,Vw,A,R,omega_m,theta_p,1)/Settings.mva/1e6;
Tm = Pw./omega_m;
% motion equation
% ---------------
DAE.f(Ddsg.omega_m) = (0.5.*(Tm-(-xd.*ids+psip).*iqs-xq.*iqs.*ids))./Hm;
% pitch control equation
% ----------------------
% vary the pitch angle only by steps of 1% of the fn
phi = round(1000*(omega_m-1))/1000;
DAE.f(Ddsg.theta_p) = (Kp.*phi-theta_p)./Tp;
% non-windup limiter
idx = find(theta_p <= 0 & DAE.f(Ddsg.theta_p) < 0);
if idx, DAE.f(Ddsg.theta_p(idx)) = 0; end
DAE.x(Ddsg.theta_p) = max(DAE.x(Ddsg.theta_p),0);
% voltage control equation
% ------------------------
DAE.f(Ddsg.idc) = (-idc+Kv.*(Vref-V))./Tv;
idx = find(idc >= idc_max & DAE.f(Ddsg.idc) > 0);
if idx, DAE.f(Ddsg.idc(idx)) = 0; end
DAE.x(Ddsg.idc) = min(idc,idc_max);
idx = find(idc <= idc_min & DAE.f(Ddsg.idc) < 0);
if idx, DAE.f(Ddsg.idc(idx)) = 0; end
DAE.x(Ddsg.idc) = max(idc,idc_min);
% generator reactive power control equations
% ------------------------------------------
DAE.f(Ddsg.ids) = (psip./xd-sqrt(psip.*psip./xd./xd-Qref./omega_m./xd)-ids)./Teq;
idx = find(ids >= ids_max & DAE.f(Ddsg.ids) > 0);
if idx, DAE.f(Ddsg.ids(idx)) = 0; end
DAE.x(Ddsg.ids) = min(ids,ids_max);
idx = find(ids <= ids_min & DAE.f(Ddsg.ids) < 0);
if idx, DAE.f(Ddsg.ids(idx)) = 0; end
DAE.x(Ddsg.ids) = max(ids,ids_min);
% speed control equations
% -----------------------
Pm = Ddsg.con(:,3).*max(min(2*omega_m-1,1),0)/Settings.mva;
DAE.f(Ddsg.iqs) = (Pm./(psip-xd.*ids)./omega_m-iqs)./Tep;
idx = find(iqs >= iqs_max & DAE.f(Ddsg.iqs) > 0);
if idx, DAE.f(Ddsg.iqs(idx)) = 0; end
DAE.x(Ddsg.iqs) = min(iqs,iqs_max);
idx = find(iqs <= iqs_min & DAE.f(Ddsg.iqs) < 0);
if idx, DAE.f(Ddsg.iqs(idx)) = 0; end
DAE.x(Ddsg.iqs) = max(iqs,iqs_min);
case 4
w = omega_m;
Vw = vw.*Wind.con(Ddsg.wind,2);
dPwdx = windpower(rho,Vw,A,R,w,theta_p,2)./Settings.mva/1e6;
Pw = windpower(rho,Vw,A,R,w,theta_p,1)/Settings.mva/1e6;
Tm = Pw./w;
Tsp = Ddsg.con(:,3).*min(2*w-1,1)./w/Settings.mva;
Tsp = max(Tsp,0);
i2Hm = 0.5./Hm;
% d f / d y
% -----------
DAE.Fy = DAE.Fy + sparse(Ddsg.idc,Ddsg.bus+Bus.n,-Kv./Tv,DAE.n,2*Bus.n);
% d g / d x
% -----------
w = omega_m;
t1 = tan(th);
c1 = cos(th);
k1 = psip-xd.*ids;
k2 = xq.*iqs;
k3 = xd.*iqs;
k4 = xq.*ids;
k5 = rs.*ids;
k6 = rs.*iqs;
DAE.Gx = DAE.Gx - sparse(Ddsg.bus,Ddsg.omega_m,k2.*ids+k1.*iqs,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(Ddsg.bus+Bus.n,Ddsg.omega_m,t1.*k2.*ids+k1.*iqs,2*Bus.n,DAE.n);
% d f / d x
% -----------
% mechanical equation
% -------------------
DAE.Fx = DAE.Fx + sparse(Ddsg.omega_m,Ddsg.omega_m,(dPwdx(:,1)-Tm).*i2Hm./w,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Ddsg.omega_m,Wind.vw(Ddsg.wind),dPwdx(:,2).*i2Hm./w,DAE.n,DAE.n);
% voltage control
% ---------------
DAE.Fx = DAE.Fx + sparse(Ddsg.idc,Ddsg.idc,-1./Tv,DAE.n,DAE.n);
idx = find(ids >= ids_min & ids <= ids_max & DAE.f(Ddsg.ids) ~= 0);
if ~isempty(idx)
DAE.Gx = DAE.Gx - sparse(Ddsg.bus(idx)+Bus.n, ...
Ddsg.idc(idx),V(idx)./c1(idx),2*Bus.n,DAE.n);
end
% reactive power control
% ----------------------
DAE.Fx = DAE.Fx + sparse(Ddsg.ids,Ddsg.ids,-1./Teq,DAE.n,DAE.n);
idx = find(ids >= ids_min & ids <= ids_max & DAE.f(Ddsg.ids) ~= 0);
if ~isempty(idx)
i1 = Ddsg.ids(idx);
i2 = Ddsg.omega_m(idx);
i3 = Ddsg.bus(idx);
o1 = (-0.5./((psip(idx)./xd(idx)).^2-Qref(idx)./w(idx)./xd(idx)).^0.5 ...
.*Qref(idx)./w(idx).^2./xd(idx))./Teq(idx);
o2 = (xd(idx)-xq(idx)).*iqs(idx).*i2Hm(idx);
o3 = -2.*k5(idx)+w(idx).*k2(idx)-w(idx).*k3(idx);
o4 = -t1(idx).*k5(idx)+t1(idx).*(-k5(idx)+w(idx).*k2(idx))-w(idx).*k3(idx);
DAE.Fx = DAE.Fx + sparse(i1,i2,o1,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(i2,i1,o2,DAE.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(i3,i1,o3,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(i3+Bus.n,i1,o4,2*Bus.n,DAE.n);
end
% speed control
% -------------
DAE.Fx = DAE.Fx + sparse(Ddsg.iqs,Ddsg.iqs,-1./Tep,DAE.n,DAE.n);
tspo = 2*Ddsg.con(:,3)/Settings.mva;
idx = find(Tsp == 0 & (2*omega_m-1) >= 1);
if ~isempty(idx), tspo(idx) = 0; end
idx = find(iqs >= iqs_min & iqs <= iqs_max & DAE.f(Ddsg.iqs) ~= 0);
if ~isempty(idx)
i1 = Ddsg.iqs(idx);
i2 = Ddsg.omega_m(idx);
i3 = Ddsg.bus(idx);
o1 = (tspo(idx)-Tsp(idx))./w(idx)./ ...
(psip(idx)-xd(idx).*ids(idx))./Teq(idx);
o4 = Pw(idx).*xd(idx)./w(idx)./(psip(idx)-xd(idx).*ids(idx)).^2./Tep(idx);
o5 = ((xd(idx)-xq(idx)).*ids(idx)-psip(idx)).*i2Hm(idx);
o6 = w(idx).*k4(idx)-2.*k6(idx)+w(idx).*k1(idx);
o7 = t1(idx).*w(idx).*k4(idx)-2.*k6(idx)+w(idx).*k1(idx);
DAE.Fx = DAE.Fx + sparse(i1,i2,o1,DAE.n,DAE.n);
%DAE.Fx = DAE.Fx + sparse(i1,Ddsg.theta_p(idx),o2,DAE.n,DAE.n);
%DAE.Fx = DAE.Fx + sparse(i1,Wind.vw(Ddsg.wind(idx)),o3,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(i1,Ddsg.ids(idx),o4,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(i2,i1,o5,DAE.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(i3,i1,o6,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(i3+Bus.n,i1,o7,2*Bus.n,DAE.n);
end
% pitch angle control equation
% ----------------------------
DAE.Fx = DAE.Fx + sparse(Ddsg.theta_p,Ddsg.theta_p,-1./Tp,DAE.n,DAE.n);
idx = find(theta_p > 0 & DAE.f(Ddsg.theta_p) ~= 0);
if ~isempty(idx)
DAE.Fx = DAE.Fx + sparse(Ddsg.theta_p(idx),Ddsg.omega_m(idx), ...
-Kp(idx)./Tp(idx),DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Ddsg.omega_m(idx),Ddsg.theta_p(idx), ...
dPwdx(idx,3).*i2Hm(idx)./omega_m(idx),DAE.n,DAE.n);
end
case 5 % non-windup limiters
idx = find(theta_p <= 0 & DAE.f(Ddsg.theta_p) <= 0);
if ~isempty(idx)
k = Ddsg.theta_p(idx);
DAE.tn(k) = 0;
DAE.Ac(:,k) = 0;
DAE.Ac(k,:) = 0;
DAE.Ac(k,k) = speye(length(idx));
end
idx = find((idc >= idc_max | idc <= idc_min) & DAE.f(Ddsg.idc) == 0);
if ~isempty(idx)
k = Ddsg.idc(idx);
DAE.tn(k) = 0;
DAE.Ac(:,k) = 0;
DAE.Ac(k,:) = 0;
DAE.Ac(k,k) = speye(length(idx));
end
idx = find((ids >= ids_max | ids <= ids_min) & DAE.f(Ddsg.ids) == 0);
if ~isempty(idx)
k = Ddsg.ids(idx);
DAE.tn(k) = 0;
DAE.Ac(:,k) = 0;
DAE.Ac(k,:) = 0;
DAE.Ac(k,k) = speye(length(idx));
end
idx = find((iqs >= iqs_max | iqs <= iqs_min) & DAE.f(Ddsg.iqs) == 0);
if ~isempty(idx)
k = Ddsg.iqs(idx);
DAE.tn(k) = 0;
DAE.Ac(:,k) = 0;
DAE.Ac(k,:) = 0;
DAE.Ac(k,k) = speye(length(idx));
end
end
%--------------------------------------------------------
function output = windpower(rho,vw,Ar,R,omega,theta,type)
%--------------------------------------------------------
lambda = omega.*R./vw;
lambdai = 1./(1./(lambda+0.08*theta)-0.035./(theta.^3+1));
switch type
case 1 % Pw
cp = 0.22*(116./lambdai-0.4*theta-5).*exp(-12.5./lambdai);
output = 0.5*rho.*cp.*Ar.*vw.^3;
case 2 % d Pw / d x
output = zeros(length(omega),3);
a1 = exp(-12.5./lambdai);
a2 = (lambda+0.08*theta).^2;
a3 = 116./lambdai-0.4*theta-5;
a4 = -9.28./(lambda+0.08*theta).^2 + ...
12.180*theta.*theta./(theta.^3+1).^2-0.4;
a5 = 1.000./(lambda+0.08*theta).^2 - ...
1.3125*theta.*theta./(theta.^3+1).^2;
% d Pw / d omega_m
output(:,1) = R.*a1.*rho.*vw.*vw.*Ar.*(-12.760+1.3750*a3)./a2;
% d Pw / d vw
output(:,2) = (omega.*R.*(12.760-1.3750*a3)./a2 ...
+ 0.330*a3.*vw).*vw.*Ar.*rho.*a1;
% d Pw / d theta_p
output(:,3) = 0.110*rho.*(a4 + a3.*a5).*a1.*Ar.*vw.^3;
end
% -------------------------------------------------------------------
% function for creating warning messages
function ddsgwarn(idx, msg)
fm_disp(strcat('Warning: DDSG #',int2str(idx),msg))
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -