📄 fm_dfig.m
字号:
pv = -ids.*st + iqs.*ct;
pt = -ids.*vqs + iqs.*vds;
%qv = ids.*ct + iqs.*st;
%qt = ids.*vds + iqs.*vqs;
%qv = -Dfig.con(:,10).*idr./Dfig.dat(:,1)-2*V./Dfig.dat(:,1);
idsv = a13.*st - a23.*ct;
idst = a13.*vqs - a23.*vds;
iqsv = -a23.*st - a13.*ct;
iqst = -a23.*vqs - a13.*vds;
k = (1-omega_m).*xm;
vdrv = k.*iqsv;
vdrt = k.*iqst;
vqrv = -k.*idsv;
vqrt = -k.*idst;
j11 = vds.*idst + vqs.*iqst + vdrt.*idr + vqrt.*iqr + pt;
j12 = vds.*idsv + vqs.*iqsv + vdrv.*idr + vqrv.*iqr + pv;
%j21 = vqs.*idst - vds.*iqst + qt;
%j22 = vqs.*idsv - vds.*iqsv + qv;
j22 = -(Dfig.con(:,10).*idr+2*V)./Dfig.dat(:,1);
DAE.J11 = DAE.J11 - sparse(Dfig.bus,Dfig.bus,j11,Bus.n,Bus.n);
DAE.J12 = DAE.J12 - sparse(Dfig.bus,Dfig.bus,j12,Bus.n,Bus.n);
%DAE.J21 = DAE.J21 - sparse(Dfig.bus,Dfig.bus,j21,Bus.n,Bus.n);
DAE.J22 = DAE.J22 - sparse(Dfig.bus,Dfig.bus,j22,Bus.n,Bus.n);
case 3
% wind speed in m/s
Vw = vw.*Wind.con(Dfig.wind,2);
Pw = windpower(rho,Vw,A,R,omega_m,theta_p,1)/Settings.mva/1e6;
% mechanical torque
Tm = Pw./omega_m;
% motion equation
% ---------------
DAE.f(Dfig.omega_m) = (Tm-xm.*(iqr.*ids-idr.*iqs)).*i2Hm;
% speed control equations
% -----------------------
Pm = Dfig.con(:,3).*max(min(2*omega_m-1,1),0)/Settings.mva;
DAE.f(Dfig.iqr) = (-(xs+xm).*Pm./V./xm./omega_m-iqr-Dfig.dat(:,7))./Te;
% non-windup limiter
%iqr_max = -Dfig.con(:,21);
%iqr_min = -Dfig.con(:,20);
iqr_max = Dfig.dat(:,8);
iqr_min = Dfig.dat(:,9);
idx = find(iqr <= iqr_min & DAE.f(Dfig.iqr) < 0);
if idx, DAE.f(Dfig.iqr(idx)) = 0; end
DAE.x(Dfig.iqr) = max(DAE.x(Dfig.iqr),iqr_min);
idx = find(iqr >= iqr_max & DAE.f(Dfig.iqr) > 0);
if idx, DAE.f(Dfig.iqr(idx)) = 0; end
DAE.x(Dfig.iqr) = min(DAE.x(Dfig.iqr),iqr_max);
% voltage control equation
% ------------------------
DAE.f(Dfig.idr) = Kv.*(V-Dfig.dat(:,6))-V./xm-idr;
% non-windup limiter
%idr_max = -Dfig.con(:,23);
%idr_min = -Dfig.con(:,22);
idr_max = Dfig.dat(:,10);
idr_min = Dfig.dat(:,11);
idx = find(idr <= idr_min & DAE.f(Dfig.idr) < 0);
if idx, DAE.f(Dfig.idr(idx)) = 0; end
DAE.x(Dfig.idr) = max(DAE.x(Dfig.idr),idr_min);
idx = find(idr >= idr_max & DAE.f(Dfig.idr) > 0);
if idx, DAE.f(Dfig.idr(idx)) = 0; end
DAE.x(Dfig.idr) = min(DAE.x(Dfig.idr),idr_max);
% pitch control equation
% ----------------------
% vary the pitch angle only by steps of 1% of the fn
phi = round(1000*(omega_m-1))/1000;
DAE.f(Dfig.theta_p) = (Kp.*phi-theta_p)./Tp;
% non-windup limiter
idx = find(theta_p <= 0 & DAE.f(Dfig.theta_p) < 0);
if idx, DAE.f(Dfig.theta_p(idx)) = 0; end
DAE.x(Dfig.theta_p) = max(DAE.x(Dfig.theta_p),0);
case 4
Vw = vw.*Wind.con(Dfig.wind,2);
dPwdx = windpower(rho,Vw,A,R,omega_m,theta_p,2)./Settings.mva/1e6;
Pw = windpower(rho,Vw,A,R,omega_m,theta_p,1)/Settings.mva/1e6;
Tm = Pw./omega_m;
Tsp = Dfig.con(:,3).*min(2*omega_m-1,1)./omega_m/Settings.mva;
Tsp = max(Tsp,0);
slip = 1-omega_m;
iqr_min = -Dfig.con(:,3)/Settings.mva;
% d f / d y
% -----------
idsv = a13.*st - a23.*ct;
idst = a13.*vqs - a23.*vds;
iqsv = -a23.*st - a13.*ct;
iqst = -a23.*vqs - a13.*vds;
ot = xm.*(idr.*iqst-iqr.*idst).*i2Hm;
ov = xm.*(idr.*iqsv-iqr.*idsv).*i2Hm;
iqrv = (xs+xm).*Tsp./Te./V./V./xm;
DAE.Fy = DAE.Fy + sparse(Dfig.omega_m,Dfig.bus,ot,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Dfig.omega_m,Dfig.bus+Bus.n,ov,DAE.n,2*Bus.n);
% d g / d x
% -----------
idsidr = -a23.*xm;
idsiqr = a13.*xm;
iqsidr = -a13.*xm;
iqsiqr = -a23.*xm;
vdridr = -rr + slip.*xm.*iqsidr;
vdriqr = slip.*(a33 + xm.*iqsiqr);
vqriqr = -rr - slip.*xm.*idsiqr;
vqridr = -slip.*(a33 + xm.*idsidr);
vdrom = -(a33.*iqr+xm.*iqs);
vqrom = a33.*idr+xm.*ids;
pidr = vds.*idsidr + vqs.*iqsidr + idr.*vdridr + vdr + iqr.*vqridr;
piqr = vds.*idsiqr + vqs.*iqsiqr + idr.*vdriqr + vqr + iqr.*vqriqr;
pom = idr.*vdrom + iqr.*vqrom;
%qidr = vqs.*idsidr - vds.*iqsidr;
%qiqr = vqs.*idsiqr - vds.*iqsiqr;
qidr = -xm.*V./Dfig.dat(:,1);
DAE.Gx = DAE.Gx - sparse(Dfig.bus,Dfig.omega_m,pom,2*Bus.n,DAE.n);
% d f / d x
% -----------
oidr = xm.*(idr.*iqsidr+iqs-iqr.*idsidr).*i2Hm;
oiqr = xm.*(idr.*iqsiqr-ids-iqr.*idsiqr).*i2Hm;
% mechanical equation
% -------------------
DAE.Fx = DAE.Fx + sparse(Dfig.omega_m,Dfig.omega_m, ...
(dPwdx(:,1)-Tm).*i2Hm./omega_m,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Dfig.omega_m,Wind.vw(Dfig.wind), ...
dPwdx(:,2).*i2Hm./omega_m,DAE.n,DAE.n);
% pitch angle control equation
% ----------------------------
DAE.Fx = DAE.Fx - sparse(Dfig.theta_p,Dfig.theta_p,1./Tp,DAE.n,DAE.n);
idx = find(theta_p > 0 & DAE.f(Dfig.theta_p) ~= 0);
if ~isempty(idx)
DAE.Fx = DAE.Fx + sparse(Dfig.theta_p(idx),Dfig.omega_m(idx), ...
Kp(idx)./Tp(idx),DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(Dfig.omega_m(idx),Dfig.theta_p(idx), ...
dPwdx(idx,3).*i2Hm(idx)./omega_m(idx),DAE.n,DAE.n);
end
% speed control equation
% ----------------------
kiqr = -(xs+xm)./V./xm./Te;
tspo = 2*Dfig.con(:,3)/Settings.mva;
idx = find(Tsp == 0 & (2*omega_m-1) >= 1);
if ~isempty(idx), tspo(idx) = 0; end
%iqr_max = -Dfig.con(:,21);
%iqr_min = -Dfig.con(:,20);
iqr_max = Dfig.dat(:,8);
iqr_min = Dfig.dat(:,9);
idx = find(iqr >= iqr_min & iqr <= iqr_max & DAE.f(Dfig.iqr) ~= 0);
if ~isempty(idx)
ki = Dfig.iqr(idx);
ko = Dfig.omega_m(idx);
kb = Dfig.bus(idx);
DAE.Fx = DAE.Fx + sparse(ki,ko,kiqr(idx).* ...
(tspo(idx)-Tsp(idx))./omega_m(idx),DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(ko,ki,oiqr(idx),DAE.n,DAE.n);
DAE.Fy = DAE.Fy + sparse(ki,kb+Bus.n,iqrv(idx),DAE.n,2*Bus.n);
%DAE.Gx = DAE.Gx - sparse(kb+Bus.n,ki,qiqr(idx),2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(kb,ki,piqr(idx),2*Bus.n,DAE.n);
end
DAE.Fx = DAE.Fx + sparse(Dfig.iqr,Dfig.iqr,-1./Te,DAE.n,DAE.n);
% voltage control equation
% ------------------------
%idr_max = -Dfig.con(:,23);
%idr_min = -Dfig.con(:,22);
idr_max = Dfig.dat(:,10);
idr_min = Dfig.dat(:,11);
idx = find(idr >= idr_min & idr <= idr_max & DAE.f(Dfig.idr) ~= 0);
if ~isempty(idx)
ki = Dfig.idr(idx);
ko = Dfig.omega_m(idx);
kb = Dfig.bus(idx);
DAE.Fx = DAE.Fx - sparse(ko,ki,oidr(idx),DAE.n,DAE.n);
DAE.Fy = DAE.Fy + sparse(ki,kb+Bus.n,Kv(idx)-1./xm(idx),DAE.n,2*Bus.n);
DAE.Gx = DAE.Gx - sparse(kb,ki,pidr(idx),2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx - sparse(kb+Bus.n,ki,qidr(idx),2*Bus.n,DAE.n);
end
DAE.Fx = DAE.Fx + sparse(Dfig.idr,Dfig.idr,-1,DAE.n,DAE.n);
case 5 % non-windup limiters
idx = find(theta_p <= 0 & DAE.f(Dfig.theta_p) <= 0);
if ~isempty(idx)
k = Dfig.theta_p(idx);
DAE.tn(k) = 0;
DAE.Ac(:,k) = 0;
DAE.Ac(k,:) = 0;
DAE.Ac(k,k) = speye(length(idx));
end
%iqr_max = -Dfig.con(:,21);
%iqr_min = -Dfig.con(:,20);
iqr_max = Dfig.dat(:,8);
iqr_min = Dfig.dat(:,9);
idx = find((iqr <= iqr_min | iqr >= iqr_max) & DAE.f(Dfig.iqr) == 0);
if ~isempty(idx)
k = Dfig.iqr(idx);
DAE.tn(k) = 0;
DAE.Ac(:,k) = 0;
DAE.Ac(k,:) = 0;
DAE.Ac(k,k) = speye(length(idx));
end
%idr_max = -Dfig.con(:,23);
%idr_min = -Dfig.con(:,22);
idr_max = Dfig.dat(:,10);
idr_min = Dfig.dat(:,11);
idx = find((idr <= idr_min | idr >= idr_max) & DAE.f(Dfig.idr) == 0);
if ~isempty(idx)
k = Dfig.idr(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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -