clear;
%clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%數(shù)據(jù)加載
n=input（‘請(qǐng)輸入要計(jì)算的節(jié)點(diǎn)系統(tǒng)（5）:‘）
????load?Node5.txt;%節(jié)點(diǎn)數(shù)據(jù)
????load?Branch5.txt;%支路數(shù)據(jù)
????load?Generator5.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+nbr1）;%松弛變量
u=1.1*ones（2*ng+n+nbr1）;%松弛變量

w=-1.5*ones（2*ng+n+nbr1）;%拉格朗日乘子
z=ones（2*ng+n+nbr1）;%拉格朗日乘子
y=zeros（2*n1）;%拉格朗日乘子
y（1:2:2*n-1）=1e-3;
y（2:2:2*n）=-1e-3;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%計(jì)算不等式約束的上下限%%%%%%%%%%%%%%%%%%%%%%%%
%gmin
gmin=zeros（2*ng+n+nbr1）;
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+nbr1）;
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（nn）;
%%%%%%%%%%%%%%%%%%%%計(jì)算非對(duì)角元素%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%計(jì)算對(duì)角元素%%%%%%%%%%%%%%%%%%%%%%
for?ii=1:n%將支路導(dǎo)納加入到對(duì)角元素中
????for?jj=1:nbr
????????if?K（jj）==0&&（Nl（jj）==ii||Nr（jj）==ii）%非變壓器支路
????????????????Y（iiii）=Y（iiii）+1/Z（jj）;
????????else?if?K（jj）~=0&&（Nl（jj）==ii||Nr（jj）==ii）%變壓器支路
????????????????Y（iiii）=Y（iiii）+1/Z（jj）/K（jj）;
????????????end
????????end
????end
end

for?ii=1:nbr%將對(duì)地電納加入到對(duì)角元素中
????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

f