内点法无功优化

%附录 优化程序源代码
clear;
%clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%数据加载
    load Node14.txt;%节点数据
    load Branch14.txt;%支路数据
    load Generator14.txt;%发电机数据
    Node=Node14;
    Branch=Branch14;
    Generator=Generator14;
%节点数据处理
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+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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算不等式约束的上下限%%%%%%%%%%%%%%%%%%%%%%%%
%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;        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成导纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=zeros（nn）;
%%%%%%%%%%%%%%%%%%%%计算非对角元素%%%%%%%%%%%%%%%%%%%%%
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（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%将对地电纳加入到对角元素中
    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（iiii）=Y