fcall.m

电力系统分析计算程序

function fcall(a)

global DAE Wind Settings

if ~a.n, return, end

omega_m = DAE.x(a.omega_m);
theta_p = DAE.x(a.theta_p);
idr = DAE.x(a.idr);
iqr = DAE.x(a.iqr);
vw = DAE.x(getidx(Wind,a.wind));
rho = getrho(Wind,a.wind);

V = DAE.y(a.vbus);
t = DAE.y(a.bus);
st = sin(t);
ct = cos(t);

rs = a.con(:,6);
xs = a.con(:,7);
rr = a.con(:,8);
xr = a.con(:,9);
xm = a.con(:,10);

i2Hm = a.u.*a.dat(:,3);
Kp = a.con(:,12);
Tp = a.con(:,13);
Kv = a.con(:,14);
Te = a.con(:,15);
R = a.dat(:,4);
A = a.dat(:,5);

as = rs.^2+a.dat(:,1).^2;
a13 = rs./as;
a23 = a.dat(:,1)./as;
a33 = a.dat(:,2);

vds = -V.*st;
vqs =  V.*ct;

ids =  -a13.*(vds-xm.*iqr)-a23.*(vqs+xm.*idr);
iqs =   a23.*(vds-xm.*iqr)-a13.*(vqs+xm.*idr);

% wind speed in m/s
iomega = 1./(omega_m+(~a.u));
Vw = vw.*getvw(Wind,a.wind);
Pw = windpower(a,rho,Vw,A,R,~a.u+omega_m,theta_p,1)/Settings.mva/1e6;
% mechanical torque
Tm = Pw.*iomega;

% motion equation
% ---------------
DAE.f(a.omega_m) = (Tm-xm.*(iqr.*ids-idr.*iqs)).*i2Hm;

% speed control equations
% -----------------------
Pm = a.con(:,3).*max(min(2*omega_m-1,1),0)/Settings.mva;
DAE.f(a.iqr) = a.u.*(-(xs+xm).*Pm./V./xm.*iomega-iqr-a.dat(:,7))./Te;
% non-windup limiter
%iqr_max = -a.con(:,21);
%iqr_min = -a.con(:,20);
iqr_max = a.dat(:,8);
iqr_min = a.dat(:,9);
idx = find(iqr <= iqr_min & DAE.f(a.iqr) < 0);
if idx, DAE.f(a.iqr(idx)) = 0; end
DAE.x(a.iqr) = max(DAE.x(a.iqr),iqr_min);
idx = find(iqr >= iqr_max & DAE.f(a.iqr) > 0);
if idx, DAE.f(a.iqr(idx)) = 0; end
DAE.x(a.iqr) = min(DAE.x(a.iqr),iqr_max);

% voltage control equation
% ------------------------
DAE.f(a.idr) = a.u.*(Kv.*(V-a.dat(:,6))-V./xm-idr);
% non-windup limiter
%idr_max = -a.con(:,23);
%idr_min = -a.con(:,22);
idr_max = a.dat(:,10);
idr_min = a.dat(:,11);
idx = find(idr <= idr_min & DAE.f(a.idr) < 0);
if idx, DAE.f(a.idr(idx)) = 0; end
DAE.x(a.idr) = max(DAE.x(a.idr),idr_min);
idx = find(idr >= idr_max & DAE.f(a.idr) > 0);
if idx, DAE.f(a.idr(idx)) = 0; end
DAE.x(a.idr) = min(DAE.x(a.idr),idr_max);

% pitch control equation
% ----------------------
% vary the pitch angle only by steps of 1% of the fn
phi = round(1000*(~a.u+omega_m-1))/1000;
DAE.f(a.theta_p) = a.u.*(Kp.*phi-theta_p)./Tp;
% non-windup limiter
idx = find(theta_p <= 0 & DAE.f(a.theta_p) < 0);
if idx, DAE.f(a.theta_p(idx)) = 0; end
DAE.x(a.theta_p) = max(DAE.x(a.theta_p),0);