內(nèi)點(diǎn)法最優(yōu)潮流MATLAB算法

內(nèi)點(diǎn)法最優(yōu)潮流MATLAB算法clear;%clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%數(shù)據(jù)加載n=input('請(qǐng)輸入要計(jì)算的節(jié)點(diǎn)系統(tǒng)⑸：’)loadNode5.txt;%節(jié)點(diǎn)數(shù)據(jù)loadBranch5.txt;%支路數(shù)據(jù)loadGenerator5.txt;%發(fā)電機(jī)數(shù)據(jù)Node=Node5;Branch=Branch5;Generator=Generator5;%節(jié)點(diǎn)數(shù)據(jù)處理N=Node(:,1);%節(jié)點(diǎn)號(hào)Type=Node(:,2);%節(jié)點(diǎn)類型Uamp=Node(:,3);%節(jié)點(diǎn)電壓幅值Dlta=Node(:,4);%節(jié)點(diǎn)電壓相角Pd=Node(:,5);%節(jié)點(diǎn)負(fù)荷有功Qd=Node(:,6);%節(jié)點(diǎn)負(fù)荷無(wú)功Pg=Node(:,7);%節(jié)點(diǎn)出力有功Qg=Node(:,8);%節(jié)點(diǎn)出力無(wú)功Umax=Node(:,9);%節(jié)點(diǎn)電壓幅值上限Umin=Node(:,10);%節(jié)點(diǎn)電壓幅值下限Bc二Node(:,11);%節(jié)點(diǎn)補(bǔ)償電容電納值％支路數(shù)據(jù)處理Nbr=Branch(:,1);%支路號(hào)Nl=Branch(:,2);%支路首節(jié)點(diǎn)Nr二Branch(:,3);%支路末節(jié)點(diǎn)R=Branch(:,4);%支路電阻X=Branch(:,5);%支路電抗Z=R+1i*X;%支路阻抗二支路電阻十支路電抗Bn二Branch(:,6);%支路對(duì)地電納K=Branch(:,7);%支路變壓器變比，0表示無(wú)變壓器Ptmax=Branch(:,8);%線路傳輸功率上限%發(fā)電機(jī)數(shù)據(jù)處理Ng二Generator(:,1);%發(fā)電機(jī)序號(hào)Nbus=Generator(:,2);%所在母線號(hào)Pumax=Generator(:,3);%發(fā)電機(jī)有功出力上界Qumax=Generator(:,4);%發(fā)電機(jī)無(wú)功出力上界Pumin二Generator(:,5);%發(fā)電機(jī)有功出力下界Qumin=Generator(:,6);%發(fā)電機(jī)無(wú)功出力下界a2二Generator(:,7);%燃料耗費(fèi)曲線二次系數(shù)a1=Generator(:,8);%燃料耗費(fèi)曲線一次系數(shù)a0=Generator(:,9);%燃料耗費(fèi)曲線常數(shù)項(xiàng)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%n=length(N);%節(jié)點(diǎn)個(gè)數(shù)ng=length(Ng);%發(fā)電機(jī)臺(tái)數(shù)nbr=length(Nbr);%支路個(gè)數(shù)x=zeros(2*(ng+n),1);%控制變量+狀態(tài)變量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=T.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;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%計(jì)算不等式約束的上下限%%%%%%%%%%%%%%%%%%%%%%%%%gmingmin=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;%gmaxgmax=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;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成導(dǎo)納矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Y=zeros(n,n);%%%%%%%%%%%%%%%%%%%%計(jì)算非對(duì)角元素%%%%%%%%%%%%%%%%%%%%%forii=1:nbrifK(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));endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%計(jì)算對(duì)角元素%%%%%%%%%%%%%%%%%%%%%%forii=1:n%將支路導(dǎo)納加入到對(duì)角元素中forjj=1:nbrifK(jj)==0&&(Nl(jj)==ii||Nr(jj)==ii)%非變壓器支路Y(ii,ii)=Y(ii,ii)+1/Z(jj);elseifK(jj)~=0&&(Nl(jj)==ii||Nr(jj)==ii)%變壓器支路Y(ii,ii)=Y(ii,ii)+1/Z(jj)/K(jj);endendendendforii=1:nbr%將對(duì)地電納加入到對(duì)角元素中ifK(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);endendforii=1:nY(ii,ii)=Y(ii,ii)+i*Bc(ii);endG=real(Y);%電導(dǎo)B=imag(Y);%電納%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%k=0;%迭代次數(shù)Kmax=150;%最大迭代次數(shù)iteration=1e-4;%誤差精度delta=0.08;Gap=(l'*z-u'*w)*ones(Kmax,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%whilek<50%計(jì)算互補(bǔ)間隙GapGap(k+1)=l'*z-u'*w;ifGap>iterationmiu=delta*Gap(k+1)/(2*(2*ng+n+nbr));%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成系數(shù)矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%相角差計(jì)算%%%%%%%%%%%%%%%%%%%%%%theta=zeros(n,n);forii=1:nforjj=1:ntheta(ii,jj)=Dlta(ii)-Dlta(jj);endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1、等式約束雅克比矩陣%%%%%%%%%%%%%%%%pxh=zeros(2*(ng+n),2*n);%%%%%%%%%%%%%%%%%%%%%%%ah/aP%%%%%%%%%%%%%%%%%%%%%%%forii=1:ngpxh(Ng(ii),2*Nbus(ii)-1)=1;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ah/aQ%%%%%%%%%%%%%%%%%%%%%%%forii=1:ngpxh(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);forii=1:nforjj=1:nifii~=jj%i!二j時(shí)的情況%非對(duì)角元素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)));%對(duì)角元素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)));endendNN(ii,ii)=NN(ii,ii)-2*Uamp(ii)*G(ii,ii);LL(ii,ii)=LL(ii,ii)+2*Uamp(ii)*B(ii,ii);endpxh(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:電源無(wú)功出力上下限約束ag2aP=zeros(ng,ng);ag2aQ=eye(ng,ng);ag2ax=zeros(2*n,ng);%g3：節(jié)點(diǎn)電壓幅值上下限約束ag3aP=zeros(ng,n);ag3aQ=zeros(ng,n);ag3ax=zeros(2*n,n);forii=1:nag3ax(2*ii,ii)=1;end%g4:線路潮流上下限約束ag4aP=zeros(ng,nbr);ag4aQ=zeros(ng,nbr);ag4ax=zeros(2*n,nbr);forii=1:nforjj=1:nbrifNl(jj)==iiag4ax(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));endifNr(jj)==iiag4ax(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))));endendendpxg=[ag1aPag2aPag3aPag4aP;ag1aQag2aQag3aQag4aQ;aglaxag2axag3axag4ax];%此即為不等式約束的雅克比矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3、對(duì)角矩陣%%%%%%%%%%%%%%%%%%%%%%%%L_1Z=zeros(2*ng+n+nbr,2*ng+n+nbr);U_1W=zeros(2*ng+n+nbr,2*ng+n+nbr);forii=1:2*ng+n+nbrL_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);forii=1:nforjj=1:n%元素位置為:12ifii~=jj%34%對(duì)角線上與ii對(duì)應(yīng)的元素%ApApb(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(i1門(mén)))))；％對(duì)角線處1號(hào)元素Apb(2*ii-1,2*ii,ii)=Apb(2*ii-1,2*ii,ii)+Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-8(打門(mén)))*。0$(m0士@(打門(mén)))))；％對(duì)角線處2號(hào)元素%%3號(hào)元素與之相等%AqAqb(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))-8(打門(mén)))*。0$(m0土2(打門(mén)))))；％對(duì)角線處1號(hào)元素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)));%對(duì)角線處2號(hào)元素%%3號(hào)元素與之相等%對(duì)角線上與jj對(duì)應(yīng)的元素%ApApb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(i1門(mén)))))；％對(duì)角線處1號(hào)元素Apb(2*jj-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%對(duì)角線處2號(hào)元素Apb(2*jj,2*jj-1,ii)=Apb(2*jj-1,2*jj,ii);%3號(hào)元素與2號(hào)元素相等%AqAqb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-8(打門(mén)))*。。$(m0士2(打門(mén)))))；％對(duì)角線處1號(hào)元素Aqb(2*jj-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%對(duì)角線處2號(hào)元素Aqb(2*jj,2*jj-1,ii)=Aqb(2*jj-1,2*jj,ii);%3號(hào)元素與2號(hào)元素相等%4號(hào)元素為0%非對(duì)角線行元素%ApApb(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)));%非對(duì)角線行處1號(hào)元素Apb(2*ii-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非對(duì)角線行處2號(hào)元素Apb(2*ii,2*jj-1,ii)=-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非對(duì)角線行處3號(hào)元素Apb(2*ii,2*jj,ii)=-(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非對(duì)角線行處4號(hào)元素%AqAqb(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)));%非對(duì)角線行處1號(hào)元素Aqb(2*ii-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非對(duì)角線行處2號(hào)元素Aqb(2*ii,2*jj-1,ii)=Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非對(duì)角線行處3號(hào)元素Aqb(2*ii,2*jj,ii)=-(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非對(duì)角線行處4號(hào)元素%非對(duì)角線列元素%ApApb(2*jj-1,2*ii-1,ii)=Apb(2*ii-1,2*jj-1,ii);%非對(duì)角線列處1號(hào)元素Apb(2*jj-1,2*ii,ii)=Apb(2*ii,2*jj-1,ii);%非對(duì)角線列處2號(hào)元素Apb(2*jj,2*ii-1,ii)=Apb(2*ii-1,2*jj,ii);%非對(duì)角線列處3號(hào)元素Apb(2*jj,2*ii,ii)=Apb(2*ii,2*jj,ii);%%非對(duì)角線列處4號(hào)元素%AqAqb(2*jj-1,2*ii-1,ii)=Aqb(2*ii-1,2*jj-1,ii);%非對(duì)角線列處1號(hào)元素Aqb(2*jj-1,2*ii,ii)=Aqb(2*ii,2*jj-1,ii);%非對(duì)角線列處2號(hào)元素Aqb(2*jj,2*ii-1,ii)=Aqb(2*ii-1,2*jj,ii);%非對(duì)角線列處3號(hào)元素Aqb(2*jj,2*ii,ii)=Aqb(2*ii,2*jj,ii);%%非對(duì)角線列處4號(hào)元素endend%對(duì)角線上與ii對(duì)應(yīng)的元素%ApApb(2*ii,2*ii-1,ii)二人口6(2*打-1,2*打，打)；％對(duì)角線處3號(hào)元素與2號(hào)元素相等人口6(2*打，2*打，㈢)=-2*6(打，打)；％對(duì)角線處4號(hào)元素%AqAqb(2*ii,2*ii-1,ii)=Aqb(2*ii-1,2*ii,ii);%對(duì)角線處3號(hào)元素與2號(hào)元素相等Aqb(2*ii,2*ii,ii)=2*B(ii,ii);%對(duì)角線處4號(hào)元素endforii=1:nA=A+Apb(:,:,ii)*y(2*ii-1)+Aqb(:,:,ii)*y(2*ii);endH2(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);forii=1:nbr%對(duì)角線上iiApc(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));%對(duì)角線上jjApc(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;%非對(duì)角線ijApc(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)));%非對(duì)角線jiApc(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);endH3(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;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成常數(shù)項(xiàng)%%%%%%%%%%%%%%%%%%%%%%%%%%Lyh=zeros(2*n,1);forii=1:nh(2*ii-1)=Pg(ii)-Pd(ii);h(2*ii)=Qg(ii)-Qd(ii);forjj=1:nh(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)));endendLy=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));forii=1:nbrgx(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));endLz=gx-l-gmin;%LwLw=gx+u-gmax;%Lle=ones(2*ng+n+nbr,1);LZ=zeros(2*ng+n+nbr,2*ng+n+nbr);forii=1:2*ng+n+nbr;LZ(ii,ii)=l(ii)*z(ii);endLl=LZ*e-miu*e;%LuUW=zeros(2*ng+n+nbr,2*ng+n+nbr);forii=1:2*ng+n+nbrUW(ii,ii)=u(ii)*w(ii);endLu=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);forii=1:2*ng+n+nbrLx41(ii)=(Ll(ii)+z(ii)*Lz(ii))/l(ii);Lx42(ii)=(Lu(ii)-w(ii)*Lw(ii))/u(ii);endLx4=pxg*(Lx41+Lx42);Lx=Lx1-Lx2-Lx3;Lxx=Lx+Lx4;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%