a.m

关于电力系统状态估计快速分解算法的MATLAB程序

clear;
p_time=0;
q_time=0;
p_chk=0;
q_chk=0;
check=0;
expt=input('请输入相角预期精度（小于1）');
v_expt=input('请输入电压幅值预期精度（小于1）');
display('开始导入数据...');
power=csvread('p_q.txt');
resist=csvread('b_g.txt');
p_b=csvread('p_b.txt');
q_b=csvread('q_b.txt');
a=csvread('angle.txt');
u=csvread('phase.txt');
v=csvread('v.txt');
R_a=csvread('R_a.txt');
R_r=csvread('R_r.txt');
display('导入成功！');
R_r=diag(R_r),
R_a=diag(R_a),
node=size(resist);
node=node(1),
g=real(resist),
b=imag(resist),
p=real(power),
q=imag(power),              %根据已知数据对电压幅值和相角进行初始化。
while check~=1
    info_p=(p_b)'*R_a*(p_b);
    info_q=(q_b)'*R_r*(q_b);               %生成信息矩阵
    for m=1:node
        for n=1:node
            pt(n)=u(m)*u(n)*(g(m,n)*cos(a(m)-a(n))+b(m,n)*sin(a(m)-a(n)));
        end
        ppt(m,m)=sum(pt);
    end                                   %节点注入有功功率
    for m=1:node
        for n=1:node
            if m<n
        ppt(m,n)=(u(m)^2-u(m)*u(n)*cos(a(m)-a(n)))*g(m,n)+b(m,n)*u(m)*u(n)*sin(a(n)-a(m));
     ppt(n,m)=(u(n)^2-u(m)*u(n)*cos(a(m)-a(n)))*g(m,n)+b(m,n)*u(m)*u(n)*sin(a(m)-a(n));
            end
        end
    end                                   %h(x)支路有功功率
    ppt=-ppt,
    count=1;
    for m=1:node
        for n=1:node
            if power(m,n)~=0
                rest_p(count,1)=p(m,n)-ppt(m,n);
                count=count+1;
            end
        end
    end                                   %生成相角残差数组
    rest_p,
    free_p=(p_b)'*R_a*rest_p,
    delt=(info_p\free_p)/(v^2),
    delt(size(delt)+1:node)=0;
    a=a+delt',
    p_time=p_time+1;
    m_delt=max(abs(delt));
    if m_delt<expt
        p_chk=1;
        if q_chk==1
            check=1;continue
        end
    end
    
    for m=1:node
        for n=1:node
            qt(n)=u(m)*u(n)*(g(m,n)*sin(a(m)-a(n))-b(m,n)*cos(a(m)-a(n)));
        end
        qqt(m,m)=sum(qt);
    end                                   %节点注入无功功率
    for m=1:node
        for n=1:node
            if m<n
                qqt(m,n)=-u(m)^2*b(m,n)+u(m)*u(n)*sin(a(n)-a(m))*g(m,n)+u(m)*u(n)*cos(a(m)-a(n))*b(m,n);
                qqt(n,m)=-u(n)^2*b(m,n)+u(m)*u(n)*cos(a(m)-a(n))*b(m,n)-u(m)*u(n)*sin(a(n)-a(m))*g(m,n);
            end
        end
    end                                    %h(x)支路无功功率
    qqt=-qqt,
    count=1;
    for m=1:node
        for n=1:node
            if power(m,n)~=0
                rest_q(count,1)=q(m,n)-qqt(m,n);
                count=count+1;
            end
        end
    end                                    %生成电压幅值残差数组
    rest_q,
    rest_v=v-u(3);
    rest_q_v=[rest_q;rest_v],
    free_q=(q_b)'*R_r*rest_q_v;
    delt_v=(info_q\free_q)/v;
    delt_v(size(delt_v)+1:node)=0;
    delt_v(1)=delt_v(1)/10;
    delt_v,
    u=u+delt_v'
    q_time=q_time+1;
    m_delt_v=max(abs(delt_v));
    if m_delt_v<v_expt
        q_chk=1;
        if p_chk==1
            check=1;continue
        end
    end
    if p_time>10*node^2
        display('fail')
        break
    end                       %为了防止不收敛情况下出现死循环
end
display('计算完毕!');
a,
u,
p_time,
q_time,                             %对计算结果进行输出