fm_sae2.m

一个较好的MATLAB潮流程序

function fm_sae2(flag)
% FM_SAE2 defines a subtransmission area equivalent
%        with two loads and two LTCs. This model
%        uses two state variables for the tap ratios
%        of the LTCs.
%
% FM_SAE2(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_SAE3
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global Bus DAE SAE2

for i = 1:SAE2.n

    if isempty(SAE2.m1)
        m1 = 1;
        m2 = 1;
    else
        m1 = DAE.x(SAE2.m1(i));
        m2 = DAE.x(SAE2.m2(i));
    end

    ha = Bus.int(round(SAE2.con(i,1)));
    hb = Bus.int(round(SAE2.con(i,2)));
    xt1 = SAE2.con(i,6);
    xt2 = SAE2.con(i,7);
    xa0 = SAE2.con(i,8);
    xb0 = SAE2.con(i,9);
    ab0 = xb0;
    a1 = SAE2.con(i,10);
    b1 = SAE2.con(i,11);
    a2 = SAE2.con(i,12);
    b2 = SAE2.con(i,13);
    h1 = SAE2.con(i,14);
    k1 = SAE2.con(i,15);
    vrif1 = SAE2.con(i,16);
    h2 = SAE2.con(i,17);
    k2 = SAE2.con(i,18);
    vrif2 = SAE2.con(i,19);
    mmax1 = SAE2.con(i,20);
    mmin1 = SAE2.con(i,21);
    mmax2 = SAE2.con(i,22);
    mmin2 = SAE2.con(i,23);

    va = DAE.V(ha);
    vb = DAE.V(hb);
    delta = DAE.a(ha);
    theta = DAE.a(hb);
    if isempty(SAE2.m1)
        hm1 = 1;
        hm2 = 1;
    else
        hm1 = SAE2.m1(i);
        hm2 = SAE2.m2(i);
    end

    switch flag
    case 0 % initialization

        if length(SAE2.con(1,:)) > 20
            xeq = zeros(4,1);
            A = [1 0 1 0; 1 0 0 1; 0 1 1 0; 0 1 0 1; 1 1 0 0; 0 0 1 1];
            x1 = SAE2.con(i,26) + SAE2.con(i,24);
            x2 = SAE2.con(i,26) + SAE2.con(i,25) + SAE2.con(i,28);
            x3 = SAE2.con(i,27) + SAE2.con(i,28) + SAE2.con(i,24);
            x4 = SAE2.con(i,27) + SAE2.con(i,25);
            x5 = SAE2.con(i,26) + SAE2.con(i,27) + SAE2.con(i,28);
            x6 = SAE2.con(i,24) + SAE2.con(i,25) + SAE2.con(i,28);
            xreal = [x1; x2; x3; x4; x5; x6];
            xeq = A\xreal;
            SAE2.con(i,8) = xeq(1);
            SAE2.con(i,9) = xeq(2);
            SAE2.con(i,6) = xeq(3);
            SAE2.con(i,7) = xeq(4);
        end

    case 1
        %  Pa :
        t1 = va*va;
        t2 = xa0*xa0;
        t4 = t1/t2;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t8 = t5/t6;
        t9 = va*vb;
        t10 = -delta+theta;
        t11 = cos(t10);
        t12 = 1/xa0;
        t13 = t11*t12;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13/ab0-t38);
        t42 = 1/xb0;
        t43 = xt1*xt1;
        t46 = xt2*xt2;
        Pa = va*t40/(t19+t29+t12+t42-1/t43*t24-1/t46*t34)*t12*(t37*(-va*t12-vb* ...
            t42*t11)-t40*vb*t42*sin(t10))/(t4+t8+2.0*t9*t13*t42);

        % C(Qa,optimized);
        t1 = va*va;
        t2 = 1/xa0;
        t4 = xa0*xa0;
        t6 = t1/t4;
        t7 = vb*vb;
        t8 = xb0*xb0;
        t10 = t7/t8;
        t11 = va*vb;
        t12 = delta-theta;
        t13 = cos(t12);
        t14 = t13*t2;
        t20 = 1/xt1;
        t21 = m1*m1;
        t25 = 1/(t20+a1/t21);
        t30 = 1/xt2;
        t31 = m2*m2;
        t35 = 1/(t30+a2/t31);
        t38 = -b1/m1*t20*t25-b2/m2*t30*t35;
        t39 = t38*t38;
        t41 = sqrt(t6+t10+2.0*t11*t14/ab0-t39);
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        Qa = t1*t2-va*t41/(t20+t30+t2+t43-1/t44*t25-1/t47*t35)*t2*(t38*vb*t43* ...
            sin(t12)+t41*(va*t2+vb*t43*t13))/(t6+t10+2.0*t11*t14*t43);

        % C(Pb,optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t4 = t1/t2;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t8 = t5/t6;
        t9 = va*vb;
        t10 = -delta+theta;
        t11 = cos(t10);
        t12 = 1/xa0;
        t13 = t11*t12;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13/ab0-t38);
        t42 = 1/xb0;
        t43 = xt1*xt1;
        t46 = xt2*xt2;
        Pb = vb*t40/(t19+t29+t12+t42-1/t43*t24-1/t46*t34)*t42*(t37*(-vb*t42-va* ...
            t12*t11)+t40*va*t12*sin(t10))/(t4+t8+2.0*t9*t13*t42);

        % C(Qb,optimized);
        t1 = vb*vb;
        t2 = 1/xb0;
        t4 = va*va;
        t5 = xa0*xa0;
        t7 = t4/t5;
        t8 = xb0*xb0;
        t10 = t1/t8;
        t11 = va*vb;
        t12 = delta-theta;
        t13 = cos(t12);
        t14 = 1/xa0;
        t15 = t13*t14;
        t21 = 1/xt1;
        t22 = m1*m1;
        t26 = 1/(t21+a1/t22);
        t31 = 1/xt2;
        t32 = m2*m2;
        t36 = 1/(t31+a2/t32);
        t39 = -b1/m1*t21*t26-b2/m2*t31*t36;
        t40 = t39*t39;
        t42 = sqrt(t7+t10+2.0*t11*t15/ab0-t40);
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        Qb = t1*t2-vb*t42/(t21+t31+t14+t2-1/t44*t26-1/t47*t36)*t2*(-t39*va*t14* ...
            sin(t12)+t42*(vb*t2+va*t14*t13))/(t7+t10+2.0*t11*t15*t2);

        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);

        % Calcolo dei termini di Jlfv

    case 2
        % C(diff(Pa,va),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t3 = 1/t2;
        t4 = t1*t3;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t8 = t5/t6;
        t9 = va*vb;
        t10 = -delta+theta;
        t11 = cos(t10);
        t12 = 1/xa0;
        t13 = t11*t12;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13*t14-t38);
        t41 = 1/xb0;
        t42 = xt1*xt1;
        t45 = xt2*xt2;
        t49 = 1/(t19+t29+t12+t41-1/t42*t24-1/t45*t34);
        t58 = t41*sin(t10);
        t61 = t12*(t37*(-va*t12-vb*t41*t11)-t40*vb*t58);
        t64 = t4+t8+2.0*t9*t13*t41;
        t65 = 1/t64;
        t68 = 1/t40;
        t71 = va*t3;
        t72 = vb*t11;
        t75 = 2.0*t71+2.0*t72*t12*t14;
        t80 = va*t40*t49;
        t89 = t64*t64;
        DAE.J12(ha,ha) = DAE.J12(ha,ha) +  t40*t49*t61*t65+va*t68*t49*t61*t65*t75/2+t80*t12*(-t37*t12-t68*vb* ...
            t58*t75/2)*t65-t80*t61/t89*(2.0*t71+2.0*t72*t12*t41);

        % C(diff(Pa,vb),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t4 = t1/t2;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t7 = 1/t6;
        t8 = t5*t7;
        t9 = va*vb;
        t10 = -delta+theta;
        t11 = cos(t10);
        t12 = 1/xa0;
        t13 = t11*t12;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13*t14-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        t51 = 1/(t19+t29+t12+t43-1/t44*t24-1/t47*t34);
        t59 = sin(t10);
        t60 = t43*t59;
        t63 = t12*(t37*(-va*t12-vb*t43*t11)-t40*vb*t60);
        t66 = t4+t8+2.0*t9*t13*t43;
        t67 = 1/t66;
        t68 = vb*t7;
        t69 = va*t11;
        t72 = 2.0*t68+2.0*t69*t12*t14;
        t77 = va*t40*t51;
        t89 = t66*t66;
        DAE.J12(ha,hb) = DAE.J12(ha,hb) + va*t41*t51*t63*t67*t72/2+t77*t12*(-t37*t43*t11-t41*vb*t60*t72/2-t40 ...
            *t43*t59)*t67-t77*t63/t89*(2.0*t68+2.0*t69*t12*t43);

        % C(diff(Pa,delta),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t3 = 1/t2;
        t4 = t1*t3;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t8 = t5/t6;
        t9 = va*vb;
        t10 = -delta+theta;
        t11 = cos(t10);
        t12 = 1/xa0;
        t13 = t11*t12;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13*t14-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        t51 = 1/(t19+t29+t12+t43-1/t44*t24-1/t47*t34);
        t52 = t51*t3;
        t59 = t40*vb;
        t60 = sin(t10);
        t61 = t43*t60;
        t63 = t37*(-va*t12-vb*t43*t11)-t59*t61;
        t66 = t4+t8+2.0*t9*t13*t43;
        t67 = 1/t66;
        t69 = vb*t60;
        t79 = t60*t60;
        t92 = t66*t66;
        DAE.J11(ha,ha) = DAE.J11(ha,ha) + t1*t41*t52*t63*t67*t69*t14+va*t40*t51*t12*(-t37*vb*t61-t41*t5*t43* ...
            t79*va*t12*t14+t59*t43*t11)*t67-2.0*t1*t40*t52*t63/t92*t69*t43;

        % C(diff(Pa,theta),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t3 = 1/t2;
        t4 = t1*t3;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t8 = t5/t6;
        t9 = va*vb;
        t10 = -delta+theta;
        t11 = cos(t10);
        t12 = 1/xa0;
        t13 = t11*t12;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13*t14-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        t51 = 1/(t19+t29+t12+t43-1/t44*t24-1/t47*t34);
        t52 = t51*t3;
        t59 = t40*vb;
        t60 = sin(t10);
        t61 = t43*t60;
        t63 = t37*(-va*t12-vb*t43*t11)-t59*t61;
        t66 = t4+t8+2.0*t9*t13*t43;
        t67 = 1/t66;
        t69 = vb*t60;
        t79 = t60*t60;
        t92 = t66*t66;
        DAE.J11(ha,hb) = DAE.J11(ha,hb) -t1*t41*t52*t63*t67*t69*t14+va*t40*t51*t12*(t37*vb*t61+t41*t5*t43* ...
            t79*va*t12*t14-t59*t43*t11)*t67+2.0*t1*t40*t52*t63/t92*t69*t43;

        % C(diff(Qa,va),optimized);
        t1 = 1/xa0;
        t2 = va*t1;
        t3 = va*va;
        t4 = xa0*xa0;
        t5 = 1/t4;
        t6 = t3*t5;
        t7 = vb*vb;
        t8 = xb0*xb0;
        t10 = t7/t8;
        t11 = va*vb;
        t12 = delta-theta;
        t13 = cos(t12);
        t14 = t13*t1;
        t15 = 1/ab0;
        t20 = 1/xt1;
        t21 = m1*m1;
        t25 = 1/(t20+a1/t21);
        t30 = 1/xt2;
        t31 = m2*m2;
        t35 = 1/(t30+a2/t31);
        t38 = -b1/m1*t20*t25-b2/m2*t30*t35;
        t39 = t38*t38;
        t41 = sqrt(t6+t10+2.0*t11*t14*t15-t39);
        t42 = 1/xb0;
        t43 = xt1*xt1;
        t46 = xt2*xt2;
        t50 = 1/(t20+t30+t1+t42-1/t43*t25-1/t46*t35);
        t58 = t2+vb*t42*t13;
        t61 = t1*(t38*vb*t42*sin(t12)+t41*t58);
        t64 = t6+t10+2.0*t11*t14*t42;
        t65 = 1/t64;
        t68 = 1/t41;
        t71 = va*t5;
        t72 = vb*t13;
        t75 = 2.0*t71+2.0*t72*t1*t15;
        t80 = va*t41*t50;
        t88 = t64*t64;
        DAE.J22(ha,ha) = DAE.J22(ha,ha) + 2.0*t2-t41*t50*t61*t65-va*t68*t50*t61*t65*t75/2-t80*t1*(t68*t58*t75 ...