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);
iqs = DAE.x(a.iqs);
idc = DAE.x(a.idc);
vw = DAE.x(getidx(Wind,a.wind));
rho = getrho(Wind,a.wind);

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

rs = a.con(:,6);
xd = a.con(:,7);
xq = a.con(:,8);
psip = a.con(:,9);

Hm = a.con(:,10);
Kp = a.con(:,11);
Tp = a.con(:,12);
Kv = a.con(:,13);
Tv = a.con(:,14);
Tep = a.con(:,15);
Teq = a.con(:,16);
R = a.dat(:,1);
A = a.dat(:,2);
Vref = a.dat(:,3);

idc_max = -a.con(:,24);
idc_min = -a.con(:,23);
iqs_max = a.con(:,21);
iqs_min = a.con(:,22);

iomega = 1./(~a.u+omega_m);

% wind speed in m/s
% -----------------
Vw = vw.*getvw(Wind,a.wind);

% mechanical torque
% -----------------
Pw = windpower(a,rho,Vw,A,R,~a.u+omega_m,theta_p,1)/Settings.mva/1e6;
Tm = Pw.*iomega;

% motion equation
% ---------------
DAE.f(a.omega_m) = 0.5*a.u.*(Tm-(psip-xd.*ids).*iqs-xq.*iqs.*ids)./Hm;

% 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) = (Kp.*phi-theta_p)./Tp;
% non-windup limiter
idx = find(theta_p <= 0 & DAE.f(a.theta_p) < 0 | (omega_m < 1 & ...
                                                  omega_m > 0));
if idx, DAE.f(a.theta_p(idx)) = 0; end
DAE.x(a.theta_p) = max(DAE.x(a.theta_p),0);

% voltage control equation
% ------------------------
DAE.f(a.idc) = a.u.*(-idc+Kv.*(Vref-V))./Tv;
idx = find(idc >= idc_max & DAE.f(a.idc) > 0);
if idx, DAE.f(a.idc(idx)) = 0; end
DAE.x(a.idc) = min(idc,idc_max);
idx = find(idc <= idc_min & DAE.f(a.idc) < 0);
if idx, DAE.f(a.idc(idx)) = 0; end
DAE.x(a.idc) = max(idc,idc_min);

% speed control equations
% -----------------------
Pm = a.con(:,3).*max(min(2*omega_m-1,1),0)/Settings.mva;
DAE.f(a.iqs) = a.u.*(Pm./(psip-xd.*ids).*iomega-iqs)./Tep;
idx = find(iqs >= iqs_max & DAE.f(a.iqs) > 0);
if idx, DAE.f(a.iqs(idx)) = 0; end
DAE.x(a.iqs) = min(iqs,iqs_max);
idx = find(iqs <= iqs_min & DAE.f(a.iqs) < 0);
if idx, DAE.f(a.iqs(idx)) = 0; end
DAE.x(a.iqs) = max(iqs,iqs_min);