fm_sae2.m

一个较好的MATLAB潮流程序

        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;
        t21 = 1/t20;
        t23 = t19+a1*t21;
        t24 = 1/t23;
        t25 = t19*t24;
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t25-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13/ab0-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t45 = 1/t44;
        t47 = xt2*xt2;
        t50 = t19+t29+t12+t43-t45*t24-1/t47*t34;
        t51 = 1/t50;
        t57 = -vb*t43-va*t12*t11;
        t60 = sin(t10);
        t63 = t37*t57+t40*va*t12*t60;
        t67 = 1/(t4+t8+2.0*t9*t13*t43);
        t71 = t20*t20;
        t74 = t23*t23;
        t75 = 1/t74;
        t79 = b1*t21*t25-2.0*b1/t71*t19*t75*a1;
        t83 = vb*t40;
        t84 = t50*t50;
        DAE.Gx(hb,hm1) = -vb*t41*t51*t43*t63*t67*t37*t79+2.0*t83/t84*t43*t63*t67*t45*t75*a1 ...
            /t20/m1+t83*t51*t43*(t79*t57-t41*va*t12*t60*t37*t79)*t67;

        % C(diff(Pb,m2),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;
        t31 = 1/t30;
        t33 = t29+a2*t31;
        t34 = 1/t33;
        t35 = t29*t34;
        t37 = -b1/m1*t19*t24-b2/m2*t35;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13/ab0-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        t48 = 1/t47;
        t50 = t19+t29+t12+t43-1/t44*t24-t48*t34;
        t51 = 1/t50;
        t57 = -vb*t43-va*t12*t11;
        t60 = sin(t10);
        t63 = t37*t57+t40*va*t12*t60;
        t67 = 1/(t4+t8+2.0*t9*t13*t43);
        t71 = t30*t30;
        t74 = t33*t33;
        t75 = 1/t74;
        t79 = b2*t31*t35-2.0*b2/t71*t29*t75*a2;
        t83 = vb*t40;
        t84 = t50*t50;
        DAE.Gx(hb,hm2) = -vb*t41*t51*t43*t63*t67*t37*t79+2.0*t83/t84*t43*t63*t67*t48*t75*a2 ...
            /t30/m2+t83*t51*t43*(t79*t57-t41*va*t12*t60*t37*t79)*t67;

        % C(diff(Qb,m1),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;
        t21 = 1/t20;
        t23 = t19+a1*t21;
        t24 = 1/t23;
        t25 = t19*t24;
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t25-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13/ab0-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t45 = 1/t44;
        t47 = xt2*xt2;
        t50 = t19+t29+t12+t43-t45*t24-1/t47*t34;
        t51 = 1/t50;
        t56 = t12*sin(t10);
        t61 = vb*t43+va*t12*t11;
        t63 = -t37*va*t56+t40*t61;
        t67 = 1/(t4+t8+2.0*t9*t13*t43);
        t71 = t20*t20;
        t74 = t23*t23;
        t75 = 1/t74;
        t79 = b1*t21*t25-2.0*b1/t71*t19*t75*a1;
        t80 = t37*t79;
        t83 = vb*t40;
        t84 = t50*t50;
        DAE.Gx(Bus.n+hb,hm1) = vb*t41*t51*t43*t63*t67*t80-2.0*t83/t84*t43*t63*t67*t45*t75*a1/t20/ ...
            m1-t83*t51*t43*(-t79*va*t56-t41*t61*t80)*t67;

        % C(diff(Qb,m2),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;
        t31 = 1/t30;
        t33 = t29+a2*t31;
        t34 = 1/t33;
        t35 = t29*t34;
        t37 = -b1/m1*t19*t24-b2/m2*t35;
        t38 = t37*t37;
        t40 = sqrt(t4+t8+2.0*t9*t13/ab0-t38);
        t41 = 1/t40;
        t43 = 1/xb0;
        t44 = xt1*xt1;
        t47 = xt2*xt2;
        t48 = 1/t47;
        t50 = t19+t29+t12+t43-1/t44*t24-t48*t34;
        t51 = 1/t50;
        t56 = t12*sin(t10);
        t61 = vb*t43+va*t12*t11;
        t63 = -t37*va*t56+t40*t61;
        t67 = 1/(t4+t8+2.0*t9*t13*t43);
        t71 = t30*t30;
        t74 = t33*t33;
        t75 = 1/t74;
        t79 = b2*t31*t35-2.0*b2/t71*t29*t75*a2;
        t80 = t37*t79;
        t83 = vb*t40;
        t84 = t50*t50;
        DAE.Gx(Bus.n+hb,hm2) = vb*t41*t51*t43*t63*t67*t80-2.0*t83/t84*t43*t63*t67*t48*t75*a2/t30/ ...
            m2-t83*t51*t43*(-t79*va*t56-t41*t61*t80)*t67;


        % calcolo dei termini di DAE.Fy

        % C(diff(dm1,va),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t3 = 1/t2;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t11 = cos(delta-theta);
        t12 = 1/xa0;
        t14 = 1/ab0;
        t17 = 1/m1;
        t19 = 1/xt1;
        t20 = m1*m1;
        t21 = 1/t20;
        t24 = 1/(t19+a1*t21);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t38 = pow(-b1*t17*t19*t24-b2/m2*t29*t34,2.0);
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t48 = pow(t19+t29+t12+1/xb0-t42*t24-1/t44*t34,2.0);
        t49 = 1/t48;
        t52 = b1*b1;
        t55 = sqrt((t1*t3+t5/t6+2.0*va*vb*t11*t12*t14-t38)*t49*t42-t52*t21);
        DAE.Fy(hm1,Bus.n+ha) = k1/t55*t17*t24*(2.0*va*t3+2.0*vb*t11*t12*t14)*t49*t42/2;

        % C(diff(dm1,vb),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t7 = 1/t6;
        t11 = cos(delta-theta);
        t12 = 1/xa0;
        t14 = 1/ab0;
        t17 = 1/m1;
        t19 = 1/xt1;
        t20 = m1*m1;
        t21 = 1/t20;
        t24 = 1/(t19+a1*t21);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t38 = pow(-b1*t17*t19*t24-b2/m2*t29*t34,2.0);
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t48 = pow(t19+t29+t12+1/xb0-t42*t24-1/t44*t34,2.0);
        t49 = 1/t48;
        t52 = b1*b1;
        t55 = sqrt((t1/t2+t5*t7+2.0*va*vb*t11*t12*t14-t38)*t49*t42-t52*t21);
        DAE.Fy(hm1,Bus.n+hb) = k1/t55*t17*t24*(2.0*vb*t7+2.0*va*t11*t12*t14)*t49*t42/2;

        % C(diff(dm1,delta),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t10 = delta-theta;
        t12 = 1/xa0;
        t14 = 1/ab0;
        t17 = 1/m1;
        t19 = 1/xt1;
        t20 = m1*m1;
        t21 = 1/t20;
        t24 = 1/(t19+a1*t21);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t38 = pow(-b1*t17*t19*t24-b2/m2*t29*t34,2.0);
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t48 = pow(t19+t29+t12+1/xb0-t42*t24-1/t44*t34,2.0);
        t49 = 1/t48;
        t52 = b1*b1;
        t55 = sqrt((t1/t2+t5/t6+2.0*va*vb*cos(t10)*t12*t14-t38)*t49*t42-t52*t21);
        DAE.Fy(hm1,ha) = -k1/t55*t17*t24*va*t12*vb*sin(t10)*t14*t49*t42;

        % C(diff(dm1,theta),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t10 = delta-theta;
        t12 = 1/xa0;
        t14 = 1/ab0;
        t17 = 1/m1;
        t19 = 1/xt1;
        t20 = m1*m1;
        t21 = 1/t20;
        t24 = 1/(t19+a1*t21);
        t29 = 1/xt2;
        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t38 = pow(-b1*t17*t19*t24-b2/m2*t29*t34,2.0);
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t48 = pow(t19+t29+t12+1/xb0-t42*t24-1/t44*t34,2.0);
        t49 = 1/t48;
        t52 = b1*b1;
        t55 = sqrt((t1/t2+t5/t6+2.0*va*vb*cos(t10)*t12*t14-t38)*t49*t42-t52*t21);
        DAE.Fy(hm1,hb) = k1/t55*t17*t24*va*t12*vb*sin(t10)*t14*t49*t42;

        % C(diff(dm2,va),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t3 = 1/t2;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t11 = cos(delta-theta);
        t12 = 1/xa0;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t27 = 1/m2;
        t29 = 1/xt2;
        t30 = m2*m2;
        t31 = 1/t30;
        t34 = 1/(t29+a2*t31);
        t38 = pow(-b1/m1*t19*t24-b2*t27*t29*t34,2.0);
        t41 = xt1*xt1;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t48 = pow(t19+t29+t12+1/xb0-1/t41*t24-t45*t34,2.0);
        t49 = 1/t48;
        t52 = b2*b2;
        t55 = sqrt((t1*t3+t5/t6+2.0*va*vb*t11*t12*t14-t38)*t49*t45-t52*t31);
        DAE.Fy(hm2,Bus.n+ha) = k2/t55*t27*t34*(2.0*va*t3+2.0*vb*t11*t12*t14)*t49*t45/2;

        % C(diff(dm2,vb),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t7 = 1/t6;
        t11 = cos(delta-theta);
        t12 = 1/xa0;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t27 = 1/m2;
        t29 = 1/xt2;
        t30 = m2*m2;
        t31 = 1/t30;
        t34 = 1/(t29+a2*t31);
        t38 = pow(-b1/m1*t19*t24-b2*t27*t29*t34,2.0);
        t41 = xt1*xt1;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t48 = pow(t19+t29+t12+1/xb0-1/t41*t24-t45*t34,2.0);
        t49 = 1/t48;
        t52 = b2*b2;
        t55 = sqrt((t1/t2+t5*t7+2.0*va*vb*t11*t12*t14-t38)*t49*t45-t52*t31);
        DAE.Fy(hm2,Bus.n+hb) = k2/t55*t27*t34*(2.0*vb*t7+2.0*va*t11*t12*t14)*t49*t45/2;

        % C(diff(dm2,delta),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t10 = delta-theta;
        t12 = 1/xa0;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t27 = 1/m2;
        t29 = 1/xt2;
        t30 = m2*m2;
        t31 = 1/t30;
        t34 = 1/(t29+a2*t31);
        t38 = pow(-b1/m1*t19*t24-b2*t27*t29*t34,2.0);
        t41 = xt1*xt1;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t48 = pow(t19+t29+t12+1/xb0-1/t41*t24-t45*t34,2.0);
        t49 = 1/t48;
        t52 = b2*b2;
        t55 = sqrt((t1/t2+t5/t6+2.0*va*vb*cos(t10)*t12*t14-t38)*t49*t45-t52*t31);
        DAE.Fy(hm2,ha) = -k2/t55*t27*t34*va*t12*vb*sin(t10)*t14*t49*t45;

        % C(diff(dm2,theta),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t10 = delta-theta;
        t12 = 1/xa0;
        t14 = 1/ab0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t27 = 1/m2;
        t29 = 1/xt2;
        t30 = m2*m2;
        t31 = 1/t30;
        t34 = 1/(t29+a2*t31);
        t38 = pow(-b1/m1*t19*t24-b2*t27*t29*t34,2.0);
        t41 = xt1*xt1;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t48 = pow(t19+t29+t12+1/xb0-1/t41*t24-t45*t34,2.0);
        t49 = 1/t48;
        t52 = b2*b2;
        t55 = sqrt((t1/t2+t5/t6+2.0*va*vb*cos(t10)*t12*t14-t38)*t49*t45-t52*t31);
        DAE.Fy(hm2,hb) = k2/t55*t27*t34*va*t12*vb*sin(t10)*t14*t49*t45;

        if ((m1 >= mmax1 | m1 <= mmin1) & DAE.f(hm1) == 0)
            DAE.Fx(hm1,:) = zeros(1,DAE.n);
            DAE.Fx(hm1,hm1) = -1;
            %                  DAE.Fx(hm2,hm1) = 0;
            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.Fx(hm1,hm2) = 0;
            DAE.Fy(hm2,:) = 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
    end


end

%===============================================================================
function  output = pow(a, b)

output = a.^b;