📄 fm_synit.m
字号:
function fm_synit
% FM_SYNIT initialize Synchornous Machines state variables,
% field voltages and mechanical torques
%
% FM_SYNIT
%
%see also FM_SYN
%
%Author: Federico Milano
%Date: 11-Nov-2002
%Version: 1.0.0
%
%E-mail: fmilano@thunderbox.uwaterloo.ca
%Web-site: http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano
global DAE Syn Bus Varname Settings jay
if Syn.n == 0; return; end
% machine dynamic orders
ord = Syn.con(:,5);
% indexes of machine orders
is2 = find(ord == 2);
is3 = find(ord == 3);
is4 = find(ord == 4);
is51 = find(ord == 5.1);
is52 = find(ord == 5.2);
is53 = find(ord == 5.3);
is6 = find(ord == 6);
is8 = find(ord == 8);
% indexes of machine buses
bs2 = Syn.bus(is2);
bs3 = Syn.bus(is3);
bs4 = Syn.bus(is4);
bs51 = Syn.bus(is51);
bs52 = Syn.bus(is52);
bs53 = Syn.bus(is53);
bs6 = Syn.bus(is6);
bs8 = Syn.bus(is8);
% constants
Syn.c1 = zeros(Syn.n,1);
Syn.c2 = zeros(Syn.n,1);
Syn.c3 = zeros(Syn.n,1);
Syn.con([is2;is53;is8],20) = 0; % Kw = 0 for these machines
Syn.con([is2;is53;is8],21) = 0; % Kp = 0 for these machines
% check parameters
% -------------------------------------------------------------------
% check inertia M and x'd
idx = find(Syn.con(:,18) <= 0);
if ~isempty(idx)
Syn.con(idx,18) = 10;
if Settings.conv
Syn.con(idx,18) = Syn.con(idx,18).*Syn.con(idx,2)/Settings.mva;
end
synwarn(idx,' Inertia cannot be <= 0. M = 10 [kWs/kVa] will be used.')
end
idx = find(Syn.con(:,9) <= 0);
if ~isempty(idx)
Syn.con(idx,9) = 0.302;
if Settings.conv
Syn.con(idx,9)= (Syn.con(idx,9)./Syn.con(idx,2))*Settings.mva;
end
synwarn(idx,' x''d cannot be <= 0. x''d = 0.302 [p.u.] will be used.')
end
% check xd, T'd0 and T"d0
is = [is3; is4; is51; is52; is53; is6; is8];
idx = find(Syn.con(is,8) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,8) = 1.90;
if Settings.conv
Syn.con(idx,8)= (Syn.con(idx,8)./Syn.con(idx,2))*Settings.mva;
end
synwarn(idx,' xd cannot be <= 0. xd = 1.90 [p.u.] will be used.')
end
idx = find(Syn.con(is,11) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,11) = 8.00;
synwarn(idx,' T''d0 cannot be <= 0. T''d0 = 8.00 [s] will be used.')
end
% check T"d0, x"d and x"q
is = [is52; is6; is8];
idx = find(Syn.con(is,12) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,12) = 0.04;
synwarn(idx,' T"d0 cannot be <= 0. T"d0 = 0.04 [s] will be used.')
end
idx = find(Syn.con(is,10) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,10) = 0.204;
if Settings.conv
Syn.con(idx,10)= (Syn.con(idx,10)./Syn.con(idx,2))*Settings.mva;
end
synwarn(idx,' x"d cannot be <= 0. x"d = 0.204 [p.u.] will be used.')
end
idx = find(Syn.con(is,15) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,15) = 0.30;
if Settings.conv
Syn.con(idx,15)= (Syn.con(idx,15)./Syn.con(idx,2))*Settings.mva;
end
synwarn(idx,' x"q cannot be <= 0. x"q = 0.30 [p.u.] will be used.')
end
% check T'q0
is = [is4; is51; is6; is8];
idx = find(Syn.con(is,16) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,16) = 0.80;
synwarn(idx,' T''q0 cannot be <= 0. T''q0 = 0.80 [s] will be used.')
end
% check T"q0
is = [is51; is52; is6; is8];
idx = find(Syn.con(is,17) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,17) = 0.02;
synwarn(idx,' T"q0 cannot be <= 0. T"q0 = 0.02 [s] will be used.')
end
% check xq
is = [is3; is4; is51; is52; is53; is6; is8];
idx = find(Syn.con(is,13) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,13) = 1.70;
if Settings.conv
Syn.con(idx,13)= (Syn.con(idx,13)./Syn.con(idx,2))*Settings.mva;
end
synwarn(idx,' xq cannot be <= 0. xq = 1.70 [p.u.] will be used.')
end
% check x'q
is = [is4; is51; is6; is8];
idx = find(Syn.con(is,14) <= 0);
if ~isempty(idx)
idx = is(idx);
Syn.con(idx,14) = 0.50;
if Settings.conv
Syn.con(idx,14)= (Syn.con(idx,14)./Syn.con(idx,2))*Settings.mva;
end
synwarn(idx,' x''q cannot be <= 0. x''q = 0.50 [p.u.] will be used.')
end
% check Taa and saturation factors
% if saturation factors are defined, Taa = 0 is used.
idx = find(Syn.con(:,25) | Syn.con(:,26));
if ~isempty(idx)
Syn.con(idx,24) = 0;
end
% parameters
xl = -Syn.con(:,6);
ra = Syn.con(:,7);
xd = Syn.con(:,8);
xd1 = Syn.con(:,9);
xd2 = Syn.con(:,10);
xq = Syn.con(:,13);
xq1 = Syn.con(:,14);
xq2 = Syn.con(:,15);
Td10 = Syn.con(:,11);
Tq10 = Syn.con(:,16);
Td20 = Syn.con(:,12);
Tq20 = Syn.con(:,17);
Kp = Syn.con(:,21);
Taa = Syn.con(:,24);
% adjusting parameters for initialization
if ~isempty(is2)
xq(is2) = xd1(is2);
end
% rotor speeds
DAE.x(Syn.omega) = ones(Syn.n,1);
% check active and reactive power ratios
synbus = sort(Syn.bus);
n_old = -1;
for i = 1:Syn.n
n_new = synbus(i);
if n_new ~= n_old
idx = find(Syn.bus == n_new);
if length(idx) == 1
if Syn.con(idx,22) ~= 1
fm_disp(['Warning: Active power ratio of ', ...
'generator #', ...
num2str(idx),' must be 1'])
Syn.con(idx,22) = 1;
end
if Syn.con(idx,23) ~= 1
fm_disp(['Warning: Reactive power ratio of ', ...
'generator #', ...
num2str(idx),' must be 1'])
Syn.con(idx,23) = 1;
end
elseif length(idx) > 1
numsyn = length(idx);
ratiop = sum(Syn.con(idx,22));
ratioq = sum(Syn.con(idx,23));
if abs(ratiop-1) > 1e-5
fm_disp(['Warning: The sum of active power ', ...
'ratios of generators #', ...
num2str(idx'),' must be 1'])
Syn.con(idx,22) = 1/numsyn;
else
Syn.con(idx(1),22) = Syn.con(idx(1),22)-(ratiop-1);
end
if abs(ratioq-1) > 1e-5
fm_disp(['Warning: The sum of reactive power ', ...
'ratios of generators #', ...
num2str(idx'),' must be 1'])
Syn.con(idx,23) = 1/numsyn;
else
Syn.con(idx(1),23) = Syn.con(idx(1),23)-(ratioq-1);
end
end
n_old = n_new;
end
end
% powers and rotor angles
Syn.Pg = -Bus.Pg(Syn.bus).*Syn.con(:,22);
Syn.Qg = -Bus.Qg(Syn.bus).*Syn.con(:,23);
Syn.Pg0 = Syn.Pg;
Vg = DAE.V(Syn.bus);
ag = DAE.a(Syn.bus);
V = Vg.*exp(jay*ag);
S = -Syn.Pg + jay*Syn.Qg;
I = S./conj(V);
delta = angle(V + (ra + jay*(xq+xl)).*I);
DAE.x(Syn.delta) = delta;
% d and q-axis voltages and currents
Vdq = V.*exp(-jay*(delta-pi/2));
Idq = I.*exp(-jay*(delta-pi/2));
vd = real(Vdq);
vq = imag(Vdq);
Syn.Id = real(Idq);
Syn.Iq = imag(Idq);
% mechanical torques/powers
Syn.pm = (vq+ra.*Syn.Iq).*Syn.Iq+(vd+ra.*Syn.Id).*Syn.Id;
% remaining state variables and field voltages
if ~isempty(is2)
K = 1./(ra(is2).^2+xd1(is2).*xd1(is2)+xd1(is2).*xl(is2)+ ...
xl(is2).*xd1(is2)+xl(is2).^2);
Syn.c1(is2) = ra(is2).*K;
Syn.c2(is2) = (xd1(is2)+xl(is2)).*K;
Syn.c3(is2) = (xd1(is2)+xl(is2)).*K;
Syn.vf(is2) = vq(is2)+ra(is2).*Syn.Iq(is2)+ ...
(xd1(is2)+xl(is2)).*Syn.Id(is2);
end
if ~isempty(is3)
K = 1./(ra(is3).^2+xq(is3).*xd1(is3)+ ...
xq(is3).*xl(is3)+xl(is3).*xd1(is3)+xl(is3).^2);
Syn.c1(is3) = ra(is3).*K;
Syn.c2(is3) = (xd1(is3)+xl(is3)).*K;
Syn.c3(is3) = (xq(is3)+xl(is3)).*K;
DAE.x(Syn.e1q(is3)) = vq(is3)+ra(is3).*Syn.Iq(is3)+ ...
(xd1(is3)+xl(is3)).*Syn.Id(is3);
Syn.vf(is3) = Sat(1,DAE.x(Syn.e1q(is3)),is3)+ ...
(xd(is3)-xd1(is3)).*Syn.Id(is3);
end
if ~isempty(is4)
K = 1./(ra(is4).^2+xq1(is4).*xd1(is4)+ ...
xq1(is4).*xl(is4)+xl(is4).*xd1(is4)+xl(is4).^2);
Syn.c1(is4) = ra(is4).*K;
Syn.c2(is4) = (xd1(is4)+xl(is4)).*K;
Syn.c3(is4) = (xq1(is4)+xl(is4)).*K;
DAE.x(Syn.e1q(is4)) = vq(is4)+ra(is4).*Syn.Iq(is4)+ ...
(xd1(is4)+xl(is4)).*Syn.Id(is4);
DAE.x(Syn.e1d(is4)) = vd(is4)+ra(is4).*Syn.Id(is4)- ...
(xq1(is4)+xl(is4)).*Syn.Iq(is4);
Syn.vf(is4) = Sat(1,DAE.x(Syn.e1q(is4)),is4)+ ...
(xd(is4)-xd1(is4)).*Syn.Id(is4);
end
if ~isempty(is51)
K = 1./(ra(is51).^2+xd1(is51).*xd1(is51)+ ...
xd1(is51).*xl(is51)+xl(is51).*xd1(is51)+xl(is51).^2);
Syn.c1(is51) = ra(is51).*K;
Syn.c2(is51) = (xd1(is51)+xl(is51)).*K;
Syn.c3(is51) = (xd1(is51)+xl(is51)).*K;
DAE.x(Syn.e1q(is51)) = vq(is51)+ra(is51).*Syn.Iq(is51)+ ...
(xd1(is51)+xl(is51)).*Syn.Id(is51);
DAE.x(Syn.e2d(is51)) = vd(is51)+ra(is51).*Syn.Id(is51)- ...
(xd1(is51)+xl(is51)).*Syn.Iq(is51);
DAE.x(Syn.e1d(is51)) = (xq(is51)-xq1(is51)- ...
Tq20(is51).*xd1(is51).* ...
(xq(is51)-xq1(is51))./Tq10(is51)./ ...
xq1(is51)).*Syn.Iq(is51);
Syn.vf(is51) = Sat(1,DAE.x(Syn.e1q(is51)),is51)+ ...
(xd(is51)-xd1(is51)).*Syn.Id(is51);
end
if ~isempty(is52)
K = 1./(ra(is52).^2+(xq2(is52)+xl(is52)).*(xd2(is52)+xl(is52)));
Syn.c1(is52) = ra(is52).*K;
Syn.c2(is52) = (xd2(is52)+xl(is52)).*K;
Syn.c3(is52) = (xq2(is52)+xl(is52)).*K;
DAE.x(Syn.e2q(is52)) = vq(is52)+ra(is52).*Syn.Iq(is52)+ ...
(xd2(is52)+xl(is52)).*Syn.Id(is52);
DAE.x(Syn.e2d(is52)) = vd(is52)+ra(is52).*Syn.Id(is52)- ...
(xq2(is52)+xl(is52)).*Syn.Iq(is52);
k1 = xd(is52)-xd1(is52)-Td20(is52).*xd2(is52).* ...
(xd(is52)-xd1(is52))./Td10(is52)./xd1(is52);
k2 = xd1(is52)-xd2(is52)+Td20(is52).*xd2(is52).* ...
(xd(is52)-xd1(is52))./Td10(is52)./xd1(is52);
%Syn.vf(is52) = (k1+k2).*Syn.Id(is52)+DAE.x(Syn.e2q(is52));
DAE.x(Syn.e1q(is52)) = -k1.*Taa(is52)./Td10(is52).* ...
Syn.Id(is52)+(1-Taa(is52)./Td10(is52)).* ...
(DAE.x(Syn.e2q(is52)) + k2.*Syn.Id(is52));
Syn.vf(is52) = (k1.*Syn.Id(is52)+Sat(1,DAE.x(Syn.e1q(is52)),is52))./ ...
(1-Taa(is52)./Td10(is52));
end
if ~isempty(is53)
DAE.x(Syn.psid(is53)) = ra(is53).*Syn.Iq(is53) + vq(is53);
DAE.x(Syn.psiq(is53)) = -(ra(is53).*Syn.Id(is53) + vd(is53));
DAE.x(Syn.e1q(is53)) = DAE.x(Syn.psid(is53)) + ...
(xd(is53)+xl(is53)).*Syn.Id(is53);
Syn.vf(is53) = DAE.x(Syn.e1q(is53));
end
if ~isempty(is6)
K = 1./(ra(is6).^2+xq2(is6).*xd2(is6)+xq2(is6).* ...
xl(is6)+xl(is6).*xd2(is6)+xl(is6).^2);
Syn.c1(is6) = ra(is6).*K;
Syn.c2(is6) = (xd2(is6)+xl(is6)).*K;
Syn.c3(is6) = (xq2(is6)+xl(is6)).*K;
DAE.x(Syn.e2q(is6)) = vq(is6)+ra(is6).*Syn.Iq(is6)+ ...
(xd2(is6)+xl(is6)).*Syn.Id(is6);
DAE.x(Syn.e2d(is6)) = vd(is6)+ra(is6).*Syn.Id(is6)- ...
(xq2(is6)+xl(is6)).*Syn.Iq(is6);
DAE.x(Syn.e1d(is6)) = (xq(is6)-xq1(is6)-Tq20(is6).* ...
xq2(is6).*(xq(is6)-xq1(is6))./ ...
Tq10(is6)./xq1(is6)).*Syn.Iq(is6);
k1 = xd(is6)-xd1(is6)-Td20(is6).*xd2(is6).* ...
(xd(is6)-xd1(is6))./Td10(is6)./xd1(is6);
k2 = xd1(is6)-xd2(is6)+Td20(is6).*xd2(is6).* ...
(xd(is6)-xd1(is6))./Td10(is6)./xd1(is6);
%Syn.vf(is6) = (k1+k2).*Syn.Id(is6)+DAE.x(Syn.e2q(is6));
DAE.x(Syn.e1q(is6)) = DAE.x(Syn.e2q(is6))+k2.*Syn.Id(is6)- ...
Taa(is6)./Td10(is6).*((k1+k2).*Syn.Id(is6)+DAE.x(Syn.e2q(is6)));
Syn.vf(is6) = (k1.*Syn.Id(is6)+Sat(1,DAE.x(Syn.e1q(is6)),is6))./ ...
(1-Taa(is6)./Td10(is6));
end
if ~isempty(is8)
DAE.x(Syn.psid(is8)) = ra(is8).*Syn.Iq(is8) + vq(is8);
DAE.x(Syn.psiq(is8)) = -(ra(is8).*Syn.Id(is8) + vd(is8));
DAE.x(Syn.e2d(is8)) = -(DAE.x(Syn.psiq(is8)) + ...
(xq2(is8)+xl(is8)).*Syn.Iq(is8));
DAE.x(Syn.e2q(is8)) = DAE.x(Syn.psid(is8)) + ...
(xd2(is8)+xl(is8)).*Syn.Id(is8);
DAE.x(Syn.e1d(is8)) = (xq(is8)-xq1(is8)-Tq20(is8).* ...
xq2(is8).*(xq(is8)-xq1(is8))./ ...
Tq10(is8)./xq1(is8)).*Syn.Iq(is8);
k1 = xd(is8)-xd1(is8)-Td20(is8).*xd2(is8).* ...
(xd(is8)-xd1(is8))./Td10(is8)./xd1(is8);
k2 = xd1(is8)-xd2(is8)+Td20(is8).*xd2(is8).* ...
(xd(is8)-xd1(is8))./Td10(is8)./xd1(is8);
%Syn.vf(is8) = (k1+k2).*Syn.Id(is8)+DAE.x(Syn.e2q(is8));
DAE.x(Syn.e1q(is8)) = -k1.*Taa(is8)./Td10(is8).* ...
Syn.Id(is8)+(1-Taa(is8)./Td10(is8)).* ...
(DAE.x(Syn.e2q(is8)) + k2.*Syn.Id(is8));
Syn.vf(is8) = (k1.*Syn.Id(is8)+Sat(1,DAE.x(Syn.e1q(is8)),is8))./ ...
(1-Taa(is8)./Td10(is8));
end
% find & delete static generators
for j = 1:Syn.n
if ~fm_rmgen(Syn.bus(j)), check = 0; end
end
fm_disp('Initialization of Synchronous Machines completed.')
% -----------------------------------------------------------
% function for creating warning messages
% -----------------------------------------------------------
function synwarn(idx, msg)
global Syn
fm_disp(strcat('Warning: Synchronous Machine #', ...
int2str(idx),'(model ',num2str(Syn.con(idx,5)), ...
') at bus #',int2str(Syn.bus(idx)),msg))
% --------------------------------------------------------------
% Saturation function
% --------------------------------------------------------------
function output = Sat(flag,e1q,isx)
global Syn
b = [0.8*ones(length(isx),1), 1-Syn.con(isx,25), 1.2*(1-Syn.con(isx,26))];
c2 = b*[12.5; -25; 12.5];
c1 = b*[-27.5; 50; -22.5];
c0 = b*[15; -24; 10];
switch flag
case 1 % saturation function
output = e1q;
idx = find(output > 0.8);
if ~isempty(idx)
output(idx)=c2(idx).*e1q(idx).^2+c1(idx).*e1q(idx)+c0(idx);
end
case 2 % Jacobian
output = ones(length(isx),1);
idx = find(e1q(isx) > 0.8);
if ~isempty(idx)
output(idx)=2*c2(idx).*e1q(idx) + c1(idx);
end
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -