For the four bus power system in power world shown below whi

For the four bus power system in power world shown below, which can be downloaded from Blackboard: Reset it to Flat Start, then using Newton -Raphson power flow method to solve the power flow (determine the voltage magnitude and angle at each bus); Using Fast Decoupled Newton -Raphson method to solve the same problem; Using DC power flow to solve the same problem. Please use table to summarize the result for each question (such as voltage magnitude, voltage angle, real power, reactive power). Please add the screenshot from Power world for each question. Using DC power flow method and the data from the Power world model to solve the same problem manually. Compare the results with Power world solution.

Solution

newton\'s

clc
clear;
n=3;
v=[1.04 1.0 1.04];
y=[5.88228-j*23.50514 -2.9427+j*11.7676 -2.9427+j*11.7676
    -2.9427+j*11.7676 5.88228-j*23.50514 -2.9427+j*11.7676
    -2.9427+j*11.7676 -2.9427+j*11.7676 5.88228-j*23.50514];
bus=ones(n,1);
qlmx=zeros(n,1);
qlmn=zeros(n,1);
vmg=zeros(n,1);
qlmx(3)=1.5;
qlmn(3)=0;
vmg(2)=1.04;
diff=10;
nf=1;
ps=[inf 0.5 -1.5];
qs=[inf 1 0];
vp=v;
while (diff>0.0001 || nf==1),
eq=1;
for i=2:n,
    bus(3)=2;
    sc(i)=0;
    sv=0;
    for k=1:n,
        sv=sv+y(i,k)*v(k);
    end
    sc(i)=v(i)*conj(sv);
    p(i)=real(sc(i));
    q(i)=imag(sc(i));
    if bus(i)==2,
        if (q(i)>qlmx(i) || q(i)<qlmn(i)),
               if(q(i)<qlmn(i)),
                   q(i)=qlmn(i);
               else
                   q(i)=qlmx(i);
               end
               bus(i)=1;
           else
               bus(i)=2;
        end
    end
    if bus(i)==1,
        er(eq)=\'p\';
        es(eq)=i;
        mm(eq)=ps(i)-p(i);
        er(eq+1)=\'q\';
        es(eq+1)=i;
        mm(eq+1)=qs(i)-q(i);
        cr(eq)=\'d\';
        cs(eq)=i;
        cr(eq+1)=\'v\';
        cs(eq+1)=i;
        eq=eq+2;
    else
        er(eq)=\'p\';
        es(eq)=i;
        cr(eq)=\'d\';
        cs(eq)=i;
        mm(eq)=ps(i)-p(i);
        eq=eq+1;
    end
end
mm
eq=eq-1;
nfq=eq;
up=zeros(eq,1);
abs(v)
abs(vp)
vp=v;
pause
for ceq=1:eq,
    for cc=1:eq,
        am=real(y(es(ceq),cs(cc))*v(cs(cc)));
        bm=imag(y(es(ceq),cs(cc))*v(cs(cc)));
        ei=real(v(es(ceq)));
        fi=imag(v(es(ceq)));
        if er(ceq)==\'p\' && cr(cc)==\'d\',
            if es(ceq)~=cs(cc),
                H=am*fi-bm*ei;
            else
                H=-q(es(ceq))-imag(y(es(ceq),cs(ceq)))*abs(v(es(ceq))^2);
            end
            jacob(ceq,cc)=H;
        end
            if er(ceq)==\'p\' && cr(cc)==\'v\',
                if es(ceq)~=cs(cc),
                    N=am*ei+bm*fi;
                else
                    N=p(es(ceq))+real(y(es(ceq),cs(ceq)))*abs(v(es(ceq))^2);
                end
                jacob(ceq,cc)=N;
            end
            if er(ceq)==\'q\' && cr(cc)==\'d\',
                if es(ceq)~=cs(cc),
                    J=-(am*ei+bm*fi);
                else
                    J=p(es(ceq))-real(y(es(ceq),cs(ceq)))*abs(v(es(ceq))^2);
                end
                jacob(ceq,cc)=J;
            end
            if er(ceq)==\'q\' && cr(cc)==\'v\',
                if es(ceq)~=cs(cc),
                    L=am*fi-bm*ei;
                else
                    L=q(es(ceq))-imag(y(es(ceq),cs(ceq)))*abs(v(es(ceq))^2);
                end
                jacob(ceq,cc)=L;
            end
    end
end
jacob
pause
up=inv(jacob)*mm\';
nfq=1;
for i=2:n,
    if bus(i)==1,
        deldif=up(nfq);
        newdelta=angle(v(i))+deldif;
        vdif=up(nfq+1)
        newvtg=abs(v(i))+vdif;
        v(i)=newvtg*(cos(newdelta)+j*sin(newdelta));
        nfq=nfq+2;
    else
        deldif=up(nfq);
        newdelta=angle(v(i))+deldif;
        v(i)=abs(v(i))*(cos(newdelta)+j*sin(newdelta));
        nfq=nfq+1;
    end
end
diff=max(abs(abs(v(2:n))-abs(vp(2:n))));
nf=nf+1;
diff
newdelta
q
nf
end

fast decoupled

clc
clear;
n=3;
v=[1.04 1.0 1.04];
y=[5.88228-j*23.50514 -2.9427+j*11.7676 -2.9427+j*11.7676
    -2.9427+j*11.7676 5.88228-j*23.50514 -2.9427+j*11.7676
    -2.9427+j*11.7676 -2.9427+j*11.7676 5.88228-j*23.50514];
bus=ones(n,1);
qlmx=zeros(n,1);
qlmn=zeros(n,1);
vmg=zeros(n,1);
qlmx(3)=1.5;
qlmn(3)=0;
vmg(2)=1.04;
diff=10;nf=1;
ps=[inf 0.5 -1.5];
qs=[inf 1 0];
vp=v;
while (diff>0.0001 || nf==1),
eq=1;
abs(v)
abs(vp)
vp=v;
for i=2:n,
    bus(3)=2;
    sc(i)=0;sv=0;
    for k=1:n,
        sv=sv+y(i,k)*v(k);
    end
    sc(i)=v(i)*conj(sv);
    p(i)=real(sc(i));
    q(i)=imag(sc(i));
    if bus(i)==2,
        if (q(i)>qlmx(i) || q(i)<qlmn(i)),
               if(q(i)<qlmn(i)),
                   q(i)=qlmn(i);
               else
                   q(i)=qlmx(i);
               end
               bus(i)=1;
           else
               bus(i)=2;
        end
    end
    if bus(i)==1,
        er(eq)=\'p\';
        es(eq)=i;
        mm(eq)=(ps(i)-p(i))/abs(v(i));
        er(eq+1)=\'q\';
        es(eq+1)=i;
        mm(eq+1)=(qs(i)-q(i))/abs(v(i));
        cr(eq)=\'d\';
        cs(eq)=i;
        cr(eq+1)=\'v\';
        cs(eq+1)=i;
        eq=eq+2;
    else
        er(eq)=\'p\';
        es(eq)=i;
        cr(eq)=\'d\';
        cs(eq)=i;
        mm(eq)=(ps(i)-p(i))/abs(v(i));
        eq=eq+1;
    end
end
mm
eq=eq-1;
nfq=eq;
up=zeros(eq,1);
for ceq=1:eq,
    for cc=1:eq,
        bm=imag(y(es(ceq),cs(cc)));
       if er(ceq)==\'p\' && cr(cc)==\'d\',
            if es(ceq)~=cs(cc),
                H=-bm;
            else
                H=-imag(y(es(ceq),cs(ceq)));
            end
            jacob(ceq,cc)=H;
        end
          
            if er(ceq)==\'q\' && cr(cc)==\'v\',
                if es(ceq)~=cs(cc),
                    L=-bm;
                else
                    L=-imag(y(es(ceq),cs(ceq)));
                end
                jacob(ceq,cc)=L;
            end
    end
end
jacob
pause
up=inv(jacob)*mm\';
nfq=1;
for i=2:n,
    if bus(i)==1,
        deldif=up(nfq);
        newdelta=angle(v(i))+deldif;
        vdif=up(nfq+1)
        newvtg=abs(v(i))+vdif;
        v(i)=newvtg*(cos(newdelta)+j*sin(newdelta));
        nfq=nfq+2;
    else
        deldif=up(nfq);
        newdelta=angle(v(i))+deldif;
        v(i)=abs(v(i))*(cos(newdelta)+j*sin(newdelta));
        nfq=nfq+1;
    end
end
diff=max(abs(abs(v(2:n))-abs(vp(2:n))));
nf=nf+1;
newdelta
q
nf
end