www.pudn.com > IPOPF.rar > opf2.m, change:2011-02-25,size:20064b


 
clear; 
%clc; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%数据加载 
n=input('请输入要计算的节点系统(5,14):') 
if n==5 
    load Node5.txt;%节点数据 
    load Branch5.txt;%支路数据 
    load Generator5.txt;%发电机数据 
    Node=Node5; 
    Branch=Branch5; 
    Generator=Generator5; 
else if n==14 
    load Node14.txt;%节点数据 
    load Branch14.txt;%支路数据 
    load Generator14.txt;%发电机数据 
    Node=Node14; 
    Branch=Branch14; 
    Generator=Generator14; 
    end 
end 
 
%节点数据处理 
N=Node(:,1);%节点号 
Type=Node(:,2);%节点类型 
Uamp=Node(:,3);%节点电压幅值 
Dlta=Node(:,4);%节点电压相角 
Pd=Node(:,5);%节点负荷有功 
Qd=Node(:,6);%节点负荷无功 
Pg=Node(:,7);%节点出力有功 
Qg=Node(:,8);%节点出力无功 
Umax=Node(:,9);%节点电压幅值上限 
Umin=Node(:,10);%节点电压幅值下限 
Bc=Node(:,11);%节点补偿电容电纳值 
%支路数据处理 
Nbr=Branch(:,1);%支路号 
Nl=Branch(:,2);%支路首节点 
Nr=Branch(:,3);%支路末节点 
R=Branch(:,4);%支路电阻 
X=Branch(:,5);%支路电抗 
Z=R+1i*X;%支路阻抗=支路电阻+支路电抗 
Bn=Branch(:,6);%支路对地电纳 
K=Branch(:,7);%支路变压器变比,0表示无变压器 
Ptmax=Branch(:,8);%线路传输功率上限 
 
%发电机数据处理 
Ng=Generator(:,1);%发电机序号 
Nbus=Generator(:,2);%所在母线号 
Pumax=Generator(:,3);%发电机有功出力上界 
Qumax=Generator(:,4);%发电机无功出力上界 
Pumin=Generator(:,5);%发电机有功出力下界 
Qumin=Generator(:,6);%发电机无功出力下界 
a2=Generator(:,7);%燃料耗费曲线二次系数 
a1=Generator(:,8);%燃料耗费曲线一次系数 
a0=Generator(:,9);%燃料耗费曲线常数项 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
n=length(N);%节点个数 
ng=length(Ng);%发电机台数 
nbr=length(Nbr);%支路个数 
 
x=zeros(2*(ng+n),1);%控制变量+状态变量 
x(1:ng)=Pg(Nbus); 
x(ng+1:2*ng)=Qg(Nbus); 
x((2*ng+2):2:2*(ng+n))=Uamp; 
x((2*ng+1):2:2*(ng+n)-1)=Dlta; 
l=0.8*ones(2*ng+n+nbr,1);%松弛变量 
u=1.1*ones(2*ng+n+nbr,1);%松弛变量 
 
w=-1.5*ones(2*ng+n+nbr,1);%拉格朗日乘子 
z=ones(2*ng+n+nbr,1);%拉格朗日乘子 
y=zeros(2*n,1);%拉格朗日乘子 
y(1:2:2*n-1)=1e-3; 
y(2:2:2*n)=-1e-3; 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算不等式约束的上下限%%%%%%%%%%%%%%%%%%%%%%%% 
%gmin 
gmin=zeros(2*ng+n+nbr,1); 
gmin(1:ng)=Pumin; 
gmin(ng+1:2*ng)=Qumin; 
gmin(2*ng+1:2*ng+n)=Umin; 
gmin(2*ng+n+1:2*ng+n+nbr)=-Ptmax;   
%gmax 
gmax=zeros(2*ng+n+nbr,1); 
gmax(1:ng)=Pumax; 
gmax(ng+1:2*ng)=Qumax; 
gmax(2*ng+1:2*ng+n)=Umax; 
gmax(2*ng+n+1:2*ng+n+nbr)=Ptmax;         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成导纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
Y=zeros(n,n); 
%%%%%%%%%%%%%%%%%%%%计算非对角元素%%%%%%%%%%%%%%%%%%%%% 
for ii=1:nbr 
    if K(ii)==0%非变压器支路 
        Y(Nl(ii),Nr(ii))=-1/Z(ii); 
        Y(Nr(ii),Nl(ii))=Y(Nl(ii),Nr(ii)); 
    else%变压器支路 
        Y(Nl(ii),Nr(ii))=-1/Z(ii)/K(ii); 
        Y(Nr(ii),Nl(ii))= Y(Nl(ii),Nr(ii)); 
    end 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
%%%%%%%%%%%%%%%%%%%%计算对角元素%%%%%%%%%%%%%%%%%%%%%% 
for ii=1:n%将支路导纳加入到对角元素中 
    for jj=1:nbr 
        if K(jj)==0&&(Nl(jj)==ii||Nr(jj)==ii)%非变压器支路 
                Y(ii,ii)=Y(ii,ii)+1/Z(jj); 
        else if K(jj)~=0&&(Nl(jj)==ii||Nr(jj)==ii)%变压器支路 
                Y(ii,ii)=Y(ii,ii)+1/Z(jj)/K(jj); 
            end 
        end 
    end 
end 
 
for ii=1:nbr%将对地电纳加入到对角元素中 
    if K(ii)==0%非变压器支路 
        Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+1i*Bn(ii); 
        Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+1i*Bn(ii); 
    else%变压器支路         
        Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+(K(ii)-1)/K(ii)/Z(ii); 
        Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+(1-K(ii))/K(ii)/K(ii)/Z(ii); 
    end 
end 
 
for ii=1:n 
    Y(ii,ii)=Y(ii,ii)+i*Bc(ii); 
end 
 
G=real(Y);%电导 
B=imag(Y);%电纳 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 
 
k=0;%迭代次数 
Kmax=150;%最大迭代次数 
iteration=1e-4;%误差精度 
delta=0.08; 
Gap=(l'*z-u'*w)*ones(Kmax,1); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
while k<50 
     
    %计算互补间隙Gap 
    Gap(k+1)=l'*z-u'*w; 
    if Gap>iteration 
        miu=delta*Gap(k+1)/(2*(2*ng+n+nbr)); 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%形成系数矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%相角差计算%%%%%%%%%%%%%%%%%%%%%% 
        theta=zeros(n,n); 
        for ii=1:n 
            for jj=1:n 
                theta(ii,jj)=Dlta(ii)-Dlta(jj); 
            end 
        end 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%1、等式约束雅克比矩阵%%%%%%%%%%%%%%%% 
        pxh=zeros(2*(ng+n),2*n); 
        %%%%%%%%%%%%%%%%%%%%%%%ah/aP%%%%%%%%%%%%%%%%%%%%%%% 
        for ii=1:ng 
            pxh(Ng(ii),2*Nbus(ii)-1)=1; 
        end 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
        %%%%%%%%%%%%%%%%%%%%%%%ah/aQ%%%%%%%%%%%%%%%%%%%%%%% 
        for ii=1:ng 
            pxh(Ng(ii)+ng,2*Nbus(ii))=1; 
        end 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%%%%%ah/ax%%%%%%%%%%%%%%%%%%%%%%% 
        HH=zeros(n,n); 
        JJ=zeros(n,n); 
        NN=zeros(n,n);         
        LL=zeros(n,n); 
 
        for ii=1:n 
            for jj=1:n 
                if ii~=jj%i!=j时的情况 
                     
                    %非对角元素 
                    HH(ii,jj)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); 
                    JJ(ii,jj)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); 
                    NN(ii,jj)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));                     
                    LL(ii,jj)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); 
                     
                    %对角元素 
                    HH(ii,ii)=HH(ii,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); 
                    JJ(ii,ii)=JJ(ii,ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); 
                    NN(ii,ii)=NN(ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); 
                    LL(ii,ii)=LL(ii,ii)-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); 
                end 
            end 
            NN(ii,ii)=NN(ii,ii)-2*Uamp(ii)*G(ii,ii); 
            LL(ii,ii)=LL(ii,ii)+2*Uamp(ii)*B(ii,ii); 
        end 
         
        pxh(1+2*ng:2:2*(n+ng)-1,1:2:2*n-1)=HH'; 
        pxh(1+2*ng:2:2*(n+ng)-1,2:2:2*n)=JJ'; 
        pxh(2+2*ng:2:2*(n+ng),1:2:2*n-1)=NN'; 
        pxh(2+2*ng:2:2*(n+ng),2:2:2*n)=LL'; 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%%%%2、不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%% 
        %g1:电源有功出力上下限约束 
        ag1aP=eye(ng,ng); 
        ag1aQ=zeros(ng,ng); 
        ag1ax=zeros(2*n,ng); 
        %g2:电源无功出力上下限约束 
        ag2aP=zeros(ng,ng); 
        ag2aQ=eye(ng,ng); 
        ag2ax=zeros(2*n,ng); 
        %g3:节点电压幅值上下限约束 
        ag3aP=zeros(ng,n); 
        ag3aQ=zeros(ng,n); 
        ag3ax=zeros(2*n,n); 
        for ii=1:n 
            ag3ax(2*ii,ii)=1; 
        end 
        %g4:线路潮流上下限约束 
        ag4aP=zeros(ng,nbr); 
        ag4aQ=zeros(ng,nbr); 
        ag4ax=zeros(2*n,nbr); 
        for ii=1:n 
            for jj=1:nbr 
                if Nl(jj)==ii 
                    ag4ax(2*ii-1,jj)=-Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))); 
                    ag4ax(2*ii,jj)=Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj))))-2*Uamp(Nl(jj))*G(Nl(jj),Nr(jj)); 
                end 
                if Nr(jj)==ii 
                    ag4ax(2*ii-1,jj)=Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))); 
                    ag4ax(2*ii,jj)=Uamp(Nl(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj)))); 
                end 
            end 
        end 
        pxg=[ag1aP ag2aP ag3aP ag4aP; 
             ag1aQ ag2aQ ag3aQ ag4aQ; 
             ag1ax ag2ax ag3ax ag4ax];%此即为不等式约束的雅克比矩阵 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%3、对角矩阵%%%%%%%%%%%%%%%%%%%%%%%% 
        L_1Z=zeros(2*ng+n+nbr,2*ng+n+nbr); 
        U_1W=zeros(2*ng+n+nbr,2*ng+n+nbr); 
        for ii=1:2*ng+n+nbr 
            L_1Z(ii,ii)=z(ii)/l(ii); 
            U_1W(ii,ii)=w(ii)/u(ii); 
        end 
               
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%海森伯矩阵%%%%%%%%%%%%%%%%%%%%%%%%%% 
        %将海森伯矩阵分为4块:H1,H2,H3,H4 
        %%%%%%%%%%%%%%%%%%%%%H1%%%%%%%%%%%%%%%%%%%%%% 
        A2=diag(a2); 
        H1=zeros(2*(ng+n),2*(ng+n)); 
        H1(1:ng,1:ng)=2*A2; 
        %%%%%%%%%%%%%%%%%%%%H2%%%%%%%%%%%%%%%%%%%%%% 
        H2=zeros(2*(ng+n),2*(ng+n)); 
        A=zeros(2*n,2*n); 
        Apb=zeros(2*n,2*n,n); 
        Aqb=zeros(2*n,2*n,n); 
                
        for ii=1:n 
            for jj=1:n                              %元素位置为:1  2 
                if ii~=jj                           %           3  4 
                %对角线上与ii对应的元素             
                    %Ap 
                    Apb(2*ii-1,2*ii-1,ii)=Apb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%对角线处1号元素 
                    Apb(2*ii-1,2*ii,ii)=Apb(2*ii-1,2*ii,ii)+Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处2号元素%%3号元素与之相等 
                    %Aq 
                    Aqb(2*ii-1,2*ii-1,ii)=Aqb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素 
                    Aqb(2*ii-1,2*ii,ii)=Aqb(2*ii-1,2*ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%对角线处2号元素%%3号元素与之相等 
 
                %对角线上与jj对应的元素 
                    %Ap 
                    Apb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%对角线处1号元素 
                    Apb(2*jj-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));          %对角线处2号元素 
                    Apb(2*jj,2*jj-1,ii)=Apb(2*jj-1,2*jj,ii);                                                        %3号元素与2号元素相等 
                    %Aq 
                    Aqb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素 
                    Aqb(2*jj-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));           %对角线处2号元素 
                    Aqb(2*jj,2*jj-1,ii)=Aqb(2*jj-1,2*jj,ii);                                                        %3号元素与2号元素相等 
                                                                                                                    %4号元素为0 
                %非对角线行元素 
                    %Ap 
                    Apb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处1号元素 
                    Apb(2*ii-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处2号元素 
                    Apb(2*ii,2*jj-1,ii)=-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处3号元素 
                    Apb(2*ii,2*jj,ii)=-(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处4号元素 
                    %Aq 
                    Aqb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处1号元素 
                    Aqb(2*ii-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处2号元素 
                    Aqb(2*ii,2*jj-1,ii)=Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处3号元素 
                    Aqb(2*ii,2*jj,ii)=-(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处4号元素 
 
                %非对角线列元素 
                    %Ap 
                    Apb(2*jj-1,2*ii-1,ii)=Apb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素 
                    Apb(2*jj-1,2*ii,ii)=Apb(2*ii,2*jj-1,ii);%非对角线列处2号元素 
                    Apb(2*jj,2*ii-1,ii)=Apb(2*ii-1,2*jj,ii);%非对角线列处3号元素 
                    Apb(2*jj,2*ii,ii)=Apb(2*ii,2*jj,ii);%%非对角线列处4号元素 
                    %Aq 
                    Aqb(2*jj-1,2*ii-1,ii)=Aqb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素 
                    Aqb(2*jj-1,2*ii,ii)=Aqb(2*ii,2*jj-1,ii);%非对角线列处2号元素 
                    Aqb(2*jj,2*ii-1,ii)=Aqb(2*ii-1,2*jj,ii);%非对角线列处3号元素 
                    Aqb(2*jj,2*ii,ii)=Aqb(2*ii,2*jj,ii);%%非对角线列处4号元素 
                end 
            end 
            %对角线上与ii对应的元素 
                %Ap 
                Apb(2*ii,2*ii-1,ii)=Apb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等 
                Apb(2*ii,2*ii,ii)=-2*G(ii,ii);%对角线处4号元素 
                %Aq 
                Aqb(2*ii,2*ii-1,ii)=Aqb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等 
                Aqb(2*ii,2*ii,ii)=2*B(ii,ii);%对角线处4号元素 
        end 
         
        for ii=1:n 
            A=A+Apb(:,:,ii)*y(2*ii-1)+Aqb(:,:,ii)*y(2*ii); 
        end 
         
        H2(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n))=A; 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%%%%%H3%%%%%%%%%%%%%%%%%%%%%% 
         
        H3=zeros(2*(ng+n),2*(ng+n)); 
        A3=zeros(2*n,2*n); 
        Apc=zeros(2*n,2*n,nbr); 
        for ii=1:nbr 
            %对角线上ii 
            Apc(2*Nl(ii)-1,2*Nl(ii)-1,ii)=-Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nl(ii)-1,2*Nl(ii),ii)=-Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nl(ii),2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nl(ii),ii); 
            Apc(2*Nl(ii),2*Nl(ii),ii)=-2*G(Nl(ii),Nr(ii)); 
             
            %对角线上jj 
            Apc(2*Nr(ii)-1,2*Nr(ii)-1,ii)=-Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nr(ii)-1,2*Nr(ii),ii)=Uamp(Nl(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nr(ii),2*Nr(ii)-1,ii)=Apc(2*Nr(ii)-1,2*Nr(ii),ii); 
            Apc(2*Nr(ii),2*Nr(ii),ii)=0; 
             
            %非对角线ij 
            Apc(2*Nl(ii)-1,2*Nr(ii)-1,ii)=Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nl(ii)-1,2*Nr(ii),ii)=-Uamp(Nl(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nl(ii),2*Nr(ii)-1,ii)=Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))); 
            Apc(2*Nl(ii),2*Nr(ii),ii)=G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))); 
             
            %非对角线ji 
            Apc(2*Nr(ii)-1,2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nr(ii)-1,ii); 
            Apc(2*Nr(ii)-1,2*Nl(ii),ii)=Apc(2*Nl(ii),2*Nr(ii)-1,ii); 
            Apc(2*Nr(ii),2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nr(ii),ii); 
            Apc(2*Nr(ii),2*Nl(ii),ii)=Apc(2*Nl(ii),2*Nr(ii),ii); 
             
            %求和 
            c=z+w; 
            A3=A3+Apc(:,:,ii)*c(2*ng+n+ii); 
        end 
         
        H3(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n))=A3; 
        %%%%%%%%%%%%%%%%%%%%%%%H4%%%%%%%%%%%%%%%%%%%%%%%%%   
        H4=pxg*(L_1Z-U_1W)*pxg'; 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
        H=-H1+H2+H3-H4; 
 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成常数项%%%%%%%%%%%%%%%%%%%%%%%%% 
         
    %Ly 
        h=zeros(2*n,1); 
        for ii=1:n 
            h(2*ii-1)=Pg(ii)-Pd(ii); 
            h(2*ii)=Qg(ii)-Qd(ii);             
            for jj=1:n 
                h(2*ii-1)=h(2*ii-1)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); 
                h(2*ii)=h(2*ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); 
            end 
        end 
        Ly=h; 
         
    %Lz 
        %g(x) 
        gx=zeros(2*ng+n+nbr,1); 
        gx(1:ng)=x(1:ng); 
        gx(ng+1:2*ng)=x(ng+1:2*ng); 
        gx(2*ng+1:2*ng+n)=x(2*ng+2:2:2*(ng+n));         
        for ii=1:nbr 
            gx(2*ng+n+ii)=Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))))-Uamp(Nl(ii))*Uamp(Nl(ii))*G(Nl(ii),Nr(ii)); 
        end   
         
        Lz=gx-l-gmin; 
         
    %Lw 
        Lw=gx+u-gmax; 
         
    %Ll 
        e=ones(2*ng+n+nbr,1); 
        LZ=zeros(2*ng+n+nbr,2*ng+n+nbr); 
        for ii=1:2*ng+n+nbr; 
            LZ(ii,ii)=l(ii)*z(ii); 
        end 
        Ll=LZ*e-miu*e;                 
    %Lu 
        UW=zeros(2*ng+n+nbr,2*ng+n+nbr); 
        for ii=1:2*ng+n+nbr 
            UW(ii,ii)=u(ii)*w(ii); 
        end 
        Lu=UW*e+miu*e;                   
    %Lx' 
        Lx1=zeros(2*(ng+n),1); 
        Lx1(1:ng)=2*a2.*x(1:ng)+a1; 
        Lx2=pxh*y; 
        Lx3=pxg*c;       
        Lx41=zeros(2*(ng+n),1); 
        Lx42=zeros(2*(ng+n),1); 
        for ii=1:2*ng+n+nbr 
            Lx41(ii)=(Ll(ii)+z(ii)*Lz(ii))/l(ii); 
            Lx42(ii)=(Lu(ii)-w(ii)*Lw(ii))/u(ii); 
        end 
        Lx4=pxg*(Lx41+Lx42); 
        Lx=Lx1-Lx2-Lx3; 
        Lxx=Lx+Lx4; 
         
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%求出修正量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
    %dx,dy 
        Hxy=[H pxh;pxh' zeros(2*n,2*n)]; 
        LxLy=[Lxx;-Ly]; 
        I=eye(2*(ng+n)+2*n); 
        dxdy=I/Hxy*LxLy; 
        dx=dxdy(1:2*(ng+n)); 
        dy=dxdy(2*(ng+n)+1:2*(ng+n)+2*n); 
    %dl 
        dl=pxg'*dx+Lz; 
    %du 
        du=-pxg'*dx-Lw; 
    %dz 
        dz=zeros(2*ng+n+nbr,1); 
        for ii=1:2*ng+n+nbr 
            dz(ii)=(-Ll(ii)-z(ii)*dl(ii))/l(ii); 
        end 
         
    %dw 
        dw=zeros(2*ng+n+nbr,1); 
        for ii=1:2*ng+n+nbr 
            dw(ii)=(-Lu(ii)-w(ii)*du(ii))/u(ii); 
        end 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%计算alfap和alfad%%%%%%%%%%%%%%%%%%%%%%%%                
       alfap=1; 
        alfad=1; 
        for ii=1:2*ng+n+nbr 
            if dl(ii)<0&&-l(ii)/dl(ii)<alfap 
                alfap=-l(ii)/dl(ii); 
            end 
            if du(ii)<0&&-u(ii)/du(ii)<alfap 
                alfap=-u(ii)/du(ii); 
            end 
            if dz(ii)<0&&-z(ii)/dz(ii)<alfad 
                alfad=-z(ii)/dz(ii); 
            end 
            if dw(ii)>0&&-w(ii)/dw(ii)<alfad 
                alfad=-w(ii)/dw(ii); 
            end 
        end 
        alfap=0.9995*alfap; 
        alfad=0.9995*alfad;       
         
        x=x+alfap*dx; 
        l=l+alfap*dl; 
        u=u+alfap*du; 
         
        y=y+alfad*dy; 
        z=z+alfad*dz;                 
        w=w+alfad*dw; 
         
        %迭代功率、电压幅值和相角 
        for ii=1:ng 
            Pg(Nbus(ii))=x(ii); 
            Qg(Nbus(ii))=x(ng+ii); 
        end 
        for ii=1:n 
            Uamp(ii)=x(2*(ng+ii)); 
            Dlta(ii)=x(2*(ng+ii)-1); 
        end 
         
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        k=k+1; 
    else 
        break; 
    end 
     
end 
 
fcost=0; 
for ii=1:ng 
    fcost=fcost+a2(ii)*Pg(Nbus(ii))*Pg(Nbus(ii))+a1(ii)*Pg(Nbus(ii))+a0(ii); 
end 
fcost 
k 
plot(0:k,Gap(1:k+1),':*'); 
Pg 
Qg 
Uamp 
 
for ii=1:n 
    if Type(ii)==3 
        Dlta=Dlta-Dlta(ii)*ones(n,1); 
    end 
end 
Dlta