fm_sae2.m

一个较好的MATLAB潮流程序

        t30 = m2*m2;
        t34 = 1/(t29+a2/t30);
        t37 = -b1/m1*t19*t24-b2/m2*t29*t34;
        t38 = t37*t37;
        t40 = sqrt(t5+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);
        t53 = t37*va;
        t54 = sin(t10);
        t55 = t12*t54;
        t60 = vb*t43+va*t12*t11;
        t62 = -t53*t55+t40*t60;
        t67 = t5+t8+2.0*t9*t13*t43;
        t68 = 1/t67;
        t91 = t67*t67;
        DAE.J21(hb,hb) = DAE.J21(hb,hb) -t1*t41*t51*t43*t62*t68*va*t55*t14-vb*t40*t51*t43*(t53*t13+t41*t60* ...
            va*vb*t54*t12*t14+t40*va*t55)*t68+2.0*t1*t40*t51*t7*t62/t91*va*t54*t12;

    case 3

        % calcolo delle equazioni differenziali DAE.f

        % C(dm1,optimized);
        t2 = va*va;
        t3 = xa0*xa0;
        t6 = vb*vb;
        t7 = xb0*xb0;
        t13 = 1/xa0;
        t18 = 1/m1;
        t20 = 1/xt1;
        t21 = m1*m1;
        t22 = 1/t21;
        t25 = 1/(t20+a1*t22);
        t30 = 1/xt2;
        t31 = m2*m2;
        t35 = 1/(t30+a2/t31);
        t39 = pow(-b1*t18*t20*t25-b2/m2*t30*t35,2.0);
        t42 = xt1*xt1;
        t43 = 1/t42;
        t45 = xt2*xt2;
        t49 = pow(t20+t30+t13+1/xb0-t43*t25-1/t45*t35,2.0);
        t53 = b1*b1;
        t56 = sqrt((t2/t3+t6/t7+2.0*va*vb*cos(delta-theta)*t13/ab0-t39)/t49*t43-t53*t22);
        DAE.f(hm1) = -h1*m1+k1*(t56*t18*t25-vrif1);

        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

        % C(dm2,optimized);
        t2 = va*va;
        t3 = xa0*xa0;
        t6 = vb*vb;
        t7 = xb0*xb0;
        t13 = 1/xa0;
        t20 = 1/xt1;
        t21 = m1*m1;
        t25 = 1/(t20+a1/t21);
        t28 = 1/m2;
        t30 = 1/xt2;
        t31 = m2*m2;
        t32 = 1/t31;
        t35 = 1/(t30+a2*t32);
        t39 = pow(-b1/m1*t20*t25-b2*t28*t30*t35,2.0);
        t42 = xt1*xt1;
        t45 = xt2*xt2;
        t46 = 1/t45;
        t49 = pow(t20+t30+t13+1/xb0-1/t42*t25-t46*t35,2.0);
        t53 = b2*b2;
        t56 = sqrt((t2/t3+t6/t7+2.0*va*vb*cos(delta-theta)*t13/ab0-t39)/t49*t46-t53*t32);
        DAE.f(hm2) = -h2*m2+k2*(t56*t28*t35-vrif2);

        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


    case 4

        %Calcolo dei termini di DAE.Fx

        % C(diff(dm1,m1),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t12 = 1/xa0;
        t17 = 1/m1;
        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*t17*t25-b2/m2*t29*t34;
        t38 = t37*t37;
        t39 = t1/t2+t5/t6+2.0*va*vb*cos(delta-theta)*t12/ab0-t38;
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t47 = t19+t29+t12+1/xb0-t42*t24-1/t44*t34;
        t48 = t47*t47;
        t49 = 1/t48;
        t52 = b1*b1;
        t55 = sqrt(t39*t49*t42-t52*t21);
        t60 = t20*t20;
        t61 = 1/t60;
        t63 = t23*t23;
        t64 = 1/t63;
        t75 = t41*t41;
        t78 = t64*a1;
        t80 = 1/t20/m1;
        DAE.Fx(hm1,hm1) = -h1+k1*(1/t55*t17*t24*(-2.0*t37*(b1*t21*t25-2.0*b1*t61*t19*t64*a1)* ...
            t49*t42+4.0*t39/t48/t47/t75*t78*t80+2.0*t52*t80)/2-t55*t21*t24+2.0*t55*t61*t78);

        % C(diff(dm1,m2),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t12 = 1/xa0;
        t17 = 1/m1;
        t19 = 1/xt1;
        t20 = m1*m1;
        t21 = 1/t20;
        t24 = 1/(t19+a1*t21);
        t29 = 1/xt2;
        t30 = m2*m2;
        t31 = 1/t30;
        t33 = t29+a2*t31;
        t34 = 1/t33;
        t35 = t29*t34;
        t37 = -b1*t17*t19*t24-b2/m2*t35;
        t38 = t37*t37;
        t39 = t1/t2+t5/t6+2.0*va*vb*cos(delta-theta)*t12/ab0-t38;
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t47 = t19+t29+t12+1/xb0-t42*t24-t45*t34;
        t48 = t47*t47;
        t49 = 1/t48;
        t52 = b1*b1;
        t55 = sqrt(t39*t49*t42-t52*t21);
        t61 = t30*t30;
        t64 = t33*t33;
        t65 = 1/t64;
        DAE.Fx(hm1,hm2) = k1/t55*t17*t24*(-2.0*t37*(b2*t31*t35-2.0*b2/t61*t29*t65*a2)*t49*t42 ...
            +4.0*t39/t48/t47*t42*t45*t65*a2/t30/m2)/2;

        % C(diff(dm2,m1),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t12 = 1/xa0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t21 = 1/t20;
        t23 = t19+a1*t21;
        t24 = 1/t23;
        t25 = t19*t24;
        t27 = 1/m2;
        t29 = 1/xt2;
        t30 = m2*m2;
        t31 = 1/t30;
        t34 = 1/(t29+a2*t31);
        t37 = -b1/m1*t25-b2*t27*t29*t34;
        t38 = t37*t37;
        t39 = t1/t2+t5/t6+2.0*va*vb*cos(delta-theta)*t12/ab0-t38;
        t41 = xt1*xt1;
        t42 = 1/t41;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t47 = t19+t29+t12+1/xb0-t42*t24-t45*t34;
        t48 = t47*t47;
        t49 = 1/t48;
        t52 = b2*b2;
        t55 = sqrt(t39*t49*t45-t52*t31);
        t61 = t20*t20;
        t64 = t23*t23;
        t65 = 1/t64;
        DAE.Fx(hm2,hm1) = k2/t55*t27*t34*(-2.0*t37*(b1*t21*t25-2.0*b1/t61*t19*t65*a1)*t49*t45 ...
            +4.0*t39/t48/t47*t45*t42*t65*a1/t20/m1)/2;

        % C(diff(dm2,m2),optimized);
        t1 = va*va;
        t2 = xa0*xa0;
        t5 = vb*vb;
        t6 = xb0*xb0;
        t12 = 1/xa0;
        t19 = 1/xt1;
        t20 = m1*m1;
        t24 = 1/(t19+a1/t20);
        t27 = 1/m2;
        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*t27*t35;
        t38 = t37*t37;
        t39 = t1/t2+t5/t6+2.0*va*vb*cos(delta-theta)*t12/ab0-t38;
        t41 = xt1*xt1;
        t44 = xt2*xt2;
        t45 = 1/t44;
        t47 = t19+t29+t12+1/xb0-1/t41*t24-t45*t34;
        t48 = t47*t47;
        t49 = 1/t48;
        t52 = b2*b2;
        t55 = sqrt(t39*t49*t45-t52*t31);
        t60 = t30*t30;
        t61 = 1/t60;
        t63 = t33*t33;
        t64 = 1/t63;
        t75 = t44*t44;
        t78 = t64*a2;
        t80 = 1/t30/m2;
        DAE.Fx(hm2,hm2) = -h2+k2*(1/t55*t27*t34*(-2.0*t37*(b2*t31*t35-2.0*b2*t61*t29*t64*a2)* ...
            t49*t45+4.0*t39/t48/t47/t75*t78*t80+2.0*t52*t80)/2-t55*t31*t34+2.0*t55*t61*t78);

        % calcolo dei termini di DAE.Gx

        % C(diff(Pa,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;
        t57 = -va*t12-vb*t43*t11;
        t60 = sin(t10);
        t63 = t37*t57-t40*vb*t43*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 = va*t40;
        t84 = t50*t50;
        DAE.Gx(ha,hm1) = -va*t41*t51*t12*t63*t67*t37*t79+2.0*t83/t84*t12*t63*t67*t45*t75*a1 ...
            /t20/m1+t83*t51*t12*(t79*t57+t41*vb*t43*t60*t37*t79)*t67;

        % C(diff(Pa,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 = -va*t12-vb*t43*t11;
        t60 = sin(t10);
        t63 = t37*t57-t40*vb*t43*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 = va*t40;
        t84 = t50*t50;
        DAE.Gx(ha,hm2) = -va*t41*t51*t12*t63*t67*t37*t79+2.0*t83/t84*t12*t63*t67*t48*t75*a2 ...
            /t30/m2+t83*t51*t12*(t79*t57+t41*vb*t43*t60*t37*t79)*t67;

        % C(diff(Qa,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 = t43*sin(t10);
        t61 = va*t12+vb*t43*t11;
        t63 = t37*vb*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 = va*t40;
        t84 = t50*t50;
        DAE.Gx(Bus.n+ha,hm1) = va*t41*t51*t12*t63*t67*t80-2.0*t83/t84*t12*t63*t67*t45*t75*a1/t20/ ...
            m1-t83*t51*t12*(t79*vb*t56-t41*t61*t80)*t67;

        % C(diff(Qa,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 = t43*sin(t10);
        t61 = va*t12+vb*t43*t11;
        t63 = t37*vb*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 = va*t40;
        t84 = t50*t50;
        DAE.Gx(Bus.n+ha,hm2) = va*t41*t51*t12*t63*t67*t80-2.0*t83/t84*t12*t63*t67*t48*t75*a2/t30/ ...
            m2-t83*t51*t12*(t79*vb*t56-t41*t61*t80)*t67;

        % C(diff(Pb,m1),optimized);
        t1 = va*va;
        t2 = xa0*xa0;