📄 fm_sae3.m
字号:
function fm_sae3(flag)
% FM_SAE3 define a subtransmission area equivalent
% with three loads and three LTCs. This model
% uses three state variables for the tap ratios
% of the LTCs.
%
% FM_SAE3(FLAG)
% FLAG = 0 initialization
% FLAG = 1 algebraic equations
% FLAG = 2 algebraic Jacobians
% FLAG = 3 differential equations
% FLAG = 4 state matrix
% FLAG = 5 non-windup limits
%
%see also FM_SAE1 and FM_SAE2
%
%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 Bus SAE3
for i = 1:SAE3.n
if isempty(SAE3.m1)
m1 = 1;
m2 = 1;
m3 = 1;
else
m1 = DAE.x(SAE3.m1(i));
m2 = DAE.x(SAE3.m2(i));
m3 = DAE.x(SAE3.m3(i));
end
ha = Bus.int(round(SAE3.con(i,1)));
hb = Bus.int(round(SAE3.con(i,2)));
xt1 = SAE3.con(i,6);
xt2 = SAE3.con(i,7);
xt3 = SAE3.con(i,8);
xa0 = SAE3.con(i,9);
xb0 = SAE3.con(i,10);
if xt1 == 0; xt1 = 0.2; end
if xt2 == 0; xt2 = 0.2; end
if xt3 == 0; xt3 = 0.2; end
if xa0 == 0; xa0 = 0.2; end
if xb0 == 0; xb0 = 0.2; end
a1 = SAE3.con(i,11);
b1 = SAE3.con(i,12);
a2 = SAE3.con(i,13);
b2 = SAE3.con(i,14);
a3 = SAE3.con(i,15);
b3 = SAE3.con(i,16);
h1 = SAE3.con(i,17);
k1 = SAE3.con(i,18);
vrif1 = SAE3.con(i,19);
h2 = SAE3.con(i,20);
k2 = SAE3.con(i,21);
vrif2 = SAE3.con(i,22);
h3 = SAE3.con(i,23);
k3 = SAE3.con(i,24);
vrif3 = SAE3.con(i,25);
mmax1 = SAE3.con(i,26);
mmin1 = SAE3.con(i,27);
mmax2 = SAE3.con(i,28);
mmin2 = SAE3.con(i,29);
mmax3 = SAE3.con(i,30);
mmin3 = SAE3.con(i,31);
va = DAE.V(ha);
vb = DAE.V(hb);
delta = DAE.a(ha);
theta = DAE.a(hb);
if isempty(SAE3.m1)
hm1 = 1;
hm2 = 1;
hm3 = 1;
else
hm1 = SAE3.m1(i);
hm2 = SAE3.m2(i);
hm3 = SAE3.m3(i);
end
t1 = 1/m1;
t3 = 1/xt1;
t4 = m1*m1;
t5 = 1/t4;
t7 = t3+a1*t5;
t8 = 1/t7;
t9 = t3*t8;
t11 = 1/m2;
t13 = 1/xt2;
t14 = m2*m2;
t15 = 1/t14;
t17 = t13+a2*t15;
t18 = 1/t17;
t19 = t13*t18;
t21 = 1/m3;
t23 = 1/xt3;
t24 = m3*m3;
t25 = 1/t24;
t27 = t23+a3*t25;
t28 = 1/t27;
t29 = t23*t28;
t31 = b1*t1*t9+b2*t11*t19+b3*t21*t29;
t32 = 1/xa0;
t33 = 1/xb0;
t34 = xt1*xt1;
t35 = 1/t34;
t37 = xt2*xt2;
t38 = 1/t37;
t40 = xt3*xt3;
t41 = 1/t40;
t43 = t3+t23+t13+t32+t33-t35*t8-t38*t18-t41*t28;
t44 = 1/t43;
t45 = t31*t44;
t46 = va*t32;
t48 = va*vb;
t49 = delta-theta;
t50 = sin(t49);
t52 = t44*t32;
t53 = t52*t33;
t54 = t48*t50*t53;
t56 = vb*t33;
t59 = xa0*xa0;
t60 = 1/t59;
t62 = t32-t60*t44;
t63 = va*va;
t65 = cos(t49);
t67 = t53*t48*t65;
t69 = xb0*xb0;
t70 = 1/t69;
t72 = t33-t70*t44;
t73 = vb*vb;
t78 = vb*t50*t53;
t81 = va*t50*t53;
t84 = t52*t56*t65;
t86 = t33*va;
t88 = t52*t86*t65;
t94 = t46+t56;
t95 = t94*t44;
t97 = t1*t3*t8;
t104 = t11*t13*t18;
t111 = t21*t23*t28;
t116 = t43*t43;
t117 = 1/t116;
t118 = t94*t117;
t119 = t4*t4;
t120 = 1/t119;
t124 = t7*t7;
t134 = 1/t124;
t136 = t3*t134*a1;
t144 = k1*t94*t117*t1*t3;
t146 = t17*t17;
t147 = 1/t146;
t148 = t147*a2;
t150 = 1/t14/m2;
t151 = t148*t150;
t155 = t27*t27;
t156 = 1/t155;
t157 = t156*a3;
t159 = 1/t24/m3;
t160 = t157*t159;
t166 = k2*t94*t117*t11*t13;
t168 = t134*a1;
t170 = 1/t4/m1;
t171 = t168*t170;
t174 = t14*t14;
t175 = 1/t174;
t189 = t13*t147*a2;
t200 = k3*t94*t117*t21*t23;
t207 = t24*t24;
t208 = 1/t207;
t222 = t23*t156*a3;
t232 = (-b1*t5*t9+2.0*b1*t120*t136)*t44;
t234 = t31*t117;
t235 = t234*t46;
t238 = t35*t134*a1*t170;
t242 = t48*t50*t117*t32;
t245 = t242*t33*t35*t171;
t252 = (-b2*t15*t19+2.0*b2*t175*t189)*t44;
t256 = t38*t147*a2*t150;
t260 = t242*t33*t38*t151;
t267 = (-b3*t25*t29+2.0*b3*t208*t222)*t44;
t271 = t41*t156*a3*t159;
t275 = t242*t33*t41*t160;
t277 = t60*t117;
t284 = t32*t117*t86*vb;
t287 = t284*t65*t35*t171;
t295 = t284*t65*t38*t151;
t303 = t284*t65*t41*t160;
t306 = t234*t56;
t315 = t70*t117;
switch flag
case 0
if length(SAE3.con(1,:)) > 28
xeq = zeros(5,1);
A = [1 0 1 0 0; 1 0 0 1 0; 1 0 0 0 1; 0 1 1 0 0; 0 1 0 1 0; ...
0 1 0 0 1; 0 0 1 1 0; 0 0 0 1 1; 0 0 1 0 1; 1 1 0 0 0];
x1 = SAE3.con(i,35) + SAE3.con(i,32);
x2 = SAE3.con(i,35) + SAE3.con(i,33) + SAE3.con(i,37);
x3 = SAE3.con(i,35) + SAE3.con(i,34) + SAE3.con(i,37) + SAE3.con(i,38);
x4 = SAE3.con(i,36) + SAE3.con(i,32) + SAE3.con(i,38) + SAE3.con(i,37);
x5 = SAE3.con(i,36) + SAE3.con(i,33) + SAE3.con(i,38);
x6 = SAE3.con(i,36) + SAE3.con(i,34);
x7 = SAE3.con(i,32) + SAE3.con(i,37) + SAE3.con(i,33);
x8 = SAE3.con(i,34) + SAE3.con(i,38) + SAE3.con(i,33);
x9 = SAE3.con(i,32) + SAE3.con(i,37) + SAE3.con(i,38) + SAE3.con(i,34);
x10= SAE3.con(i,35) + SAE3.con(i,37) + SAE3.con(i,38) + SAE3.con(i,36);
xreal = [x1; x2; x3; x4; x5; x6; x7; x8; x9; x10];
xeq = A\xreal;
SAE3.con(i,9) = xeq(1);
SAE3.con(i,10) = xeq(2);
SAE3.con(i,6) = xeq(3);
SAE3.con(i,7) = xeq(4);
SAE3.con(i,8) = xeq(5);
end
case 1
Pa = t45*t46+t54;
Pb = t45*t56-t54;
Qa = t62*t63-t67;
Qb = t72*t73-t67;
DAE.gp(ha) = Pa + DAE.gp(ha);
DAE.gq(ha) = Qa + DAE.gq(ha);
DAE.gp(hb) = Pb + DAE.gp(hb);
DAE.gq(hb) = Qb + DAE.gq(hb);
case 2
DAE.J12(ha,ha) = DAE.J12(ha,ha) + t45*t32+t78;
DAE.J12(ha,hb) = DAE.J12(ha,hb) + t81;
DAE.J11(ha,ha) = DAE.J11(ha,ha) + t67;
DAE.J11(ha,hb) = DAE.J11(ha,hb) - t67;
DAE.J22(ha,ha) = DAE.J22(ha,ha) + 2.0*t62*va-t84;
DAE.J22(ha,hb) = DAE.J22(ha,hb) - t88;
DAE.J21(ha,ha) = DAE.J21(ha,ha) + t54;
DAE.J21(ha,hb) = DAE.J21(ha,hb) - t54;
DAE.J12(hb,ha) = DAE.J12(hb,ha) - t78;
DAE.J12(hb,hb) = DAE.J12(hb,hb) + t45*t33-t81;
DAE.J11(hb,ha) = DAE.J11(hb,ha) - t67;
DAE.J11(hb,hb) = DAE.J11(hb,hb) + t67;
DAE.J22(hb,ha) = DAE.J22(hb,ha) - t84;
DAE.J22(hb,hb) = DAE.J22(hb,hb) + 2.0*t72*vb-t88;
DAE.J21(hb,ha) = DAE.J21(hb,ha) + t54;
DAE.J21(hb,hb) = DAE.J21(hb,hb) - t54;
case 3
DAE.f(hm1) = -h1*m1+k1*(t95*t97-vrif1);
DAE.f(hm2) = -h2*m2+k2*(t95*t104-vrif2);
DAE.f(hm3) = -h3*m3+k3*(t95*t111-vrif3);
if (m1 >= mmax1)
m1 = mmax1;
DAE.x(hm1) = m1;
if(DAE.f(hm1) > 0); DAE.f(hm1) = 0;end
end
if (m1 <= mmin1)
m1 = mmin1;
DAE.x(hm1) = m1;
if(DAE.f(hm1) < 0); DAE.f(hm1) = 0;end
end
if (m2 >= mmax2)
m2 = mmax2;
DAE.x(hm2) = m2;
if(DAE.f(hm2) > 0); DAE.f(hm2) = 0;end
end
if (m2 <= mmin2)
m2 = mmin2;
DAE.x(hm2) = m2;
if(DAE.f(hm2) < 0); DAE.f(hm2) = 0;end
end
if (m3 >= mmax3)
m3 = mmax3;
DAE.x(hm3) = m3;
if(DAE.f(hm3) > 0); DAE.f(hm3) = 0;end
end
if (m3 <= mmin3)
m3 = mmin3;
DAE.x(hm3) = m3;
if(DAE.f(hm3) < 0); DAE.f(hm3) = 0;end
end
case 4
DAE.Fx(hm1,hm1) = -h1+k1*(2.0*t118*t120/t34/xt1/t124/t7*a1-t95*t5*t3*t8+2.0*t95*t120*t136);
DAE.Fx(hm1,hm2) = 2.0*t144*t8*t38*t151;
DAE.Fx(hm1,hm3) = 2.0*t144*t8*t41*t160;
DAE.Fx(hm2,hm1) = 2.0*t166*t18*t35*t171;
DAE.Fx(hm2,hm2) = -h2+k2*(2.0*t118*t175/t37/xt2/t146/t17*a2-t95*t15*t13*t18+2.0*t95*t175*t189);
DAE.Fx(hm2,hm3) = 2.0*t166*t18*t41*t160;
DAE.Fx(hm3,hm1) = 2.0*t200*t28*t35*t171;
DAE.Fx(hm3,hm2) = 2.0*t200*t28*t38*t151;
DAE.Fx(hm3,hm3) = -h3+k3*(2.0*t118*t208/t40/xt3/t155/t27*a3-t95*t25*t23*t28+2.0*t95*t208*t222);
DAE.Gx(ha,hm1) = t232*t46+2.0*t235*t238+2.0*t245;
DAE.Gx(ha,hm2) = t252*t46+2.0*t235*t256+2.0*t260;
DAE.Gx(ha,hm3) = t267*t46+2.0*t235*t271+2.0*t275;
DAE.Gx(ha+Bus.n,hm1) = -2.0*t277*t35*t168*t170*t63-2.0*t287;
DAE.Gx(ha+Bus.n,hm2) = -2.0*t277*t38*t148*t150*t63-2.0*t295;
DAE.Gx(ha+Bus.n,hm3) = -2.0*t277*t41*t157*t159*t63-2.0*t303;
DAE.Gx(hb,hm1) = t232*t56+2.0*t306*t238-2.0*t245;
DAE.Gx(hb,hm2) = t252*t56+2.0*t306*t256-2.0*t260;
DAE.Gx(hb,hm3) = t267*t56+2.0*t306*t271-2.0*t275;
DAE.Gx(hb+Bus.n,hm1) = -2.0*t315*t35*t168*t170*t73-2.0*t287;
DAE.Gx(hb+Bus.n,hm2) = -2.0*t315*t38*t148*t150*t73-2.0*t295;
DAE.Gx(hb+Bus.n,hm3) = -2.0*t315*t41*t157*t159*t73-2.0*t303;
DAE.Fy(hm1,Bus.n+ha) = k1*t32*t44*t97;
DAE.Fy(hm1,Bus.n+hb) = k1*t33*t44*t97;
DAE.Fy(hm2,Bus.n+ha) = k2*t32*t44*t104;
DAE.Fy(hm2,Bus.n+hb) = k2*t33*t44*t104;
DAE.Fy(hm3,Bus.n+ha) = k3*t32*t44*t111;
DAE.Fy(hm3,Bus.n+hb) = k3*t33*t44*t111;
if ((m1 >= mmax1 | m1 <= mmin1) & DAE.f(hm1) == 0)
DAE.Fx(hm1,:) = zeros(1,DAE.n);
DAE.Fx(hm1,hm1) = -1;
DAE.Fy(hm1,:) = zeros(1,2*Bus.n);
end
if ((m2 >= mmax2 | m2 <= mmin2) & DAE.f(hm2) == 0)
DAE.Fx(hm2,:) = zeros(1,DAE.n);
DAE.Fx(hm2,hm2) = -1;
DAE.Fy(hm2,:) = zeros(1,2*Bus.n);
end
if ((m3 >= mmax3 | m3 <= mmin3) & DAE.f(hm3) == 0)
DAE.Fx(hm3,:) = zeros(1,DAE.n);
DAE.Fx(hm3,hm3) = -1;
DAE.Fy(hm3,:) = zeros(1,2*Bus.n);
end
case 5
if ((m1 >= mmax1 | m1 <= mmin1) & DAE.f(hm1) == 0)
DAE.tn(hm1) = 0;
DAE.Ac(hm1,:) = zeros(1,DAE.n+2*Bus.n);
DAE.Ac(hm1,hm1) = DAE.Fx(hm1,hm1);
end
if ((m2 >= mmax2 | m2 <= mmin2) & DAE.f(hm2) == 0)
DAE.tn(hm2) = 0;
DAE.Ac(hm2,:) = zeros(1,DAE.n+2*Bus.n);
DAE.Ac(hm2,hm2) = DAE.Fx(hm2,hm2);
end
if ((m3 >= mmax3 | m3 <= mmin3) & DAE.f(hm3) == 0)
DAE.tn(hm3) = 0;
DAE.Ac(hm3,:) = zeros(1,DAE.n+2*Bus.n);
DAE.Ac(hm3,hm3) = DAE.Fx(hm3,hm3);
end
end
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -