%*******************************************************
%            程序名称：电力系统无功优化程序
%            程序功能：主程序
%            程序算法：现代内点法，采用Sigmoid函数处理离散变量
%            目标函数：发电费用最小（1）、有功出力最小（2）、有功网损最小（3）、无功网损最小（4）
%            编者：    韦园清
%            编写时间：2010年10月
%*******************************************************
function []=OPF(filename)
clc;
%% 读取数据
data=dlmread(filename);
if data(1,1)<1780
    %读取1780节点系统格式类型的数据
    [nodeNum,lineNum,capacity,kmax,centrPara,accuracy,opfGoal,balNum,balNode,mtrLine,mtrGround,mtrTrans,runPara,Pg,Qg,Pd,Qd,pvNum,pvNode,pvV,...
        pvQmin,pvQmax,pgNumall,pgNodeall,pgNum,pgNode,a,b,c,Pmin,Pmax,capNum,capI,capGi,capBi,capK,capKmin,capKmax,Vmax,Vmin]=readData(filename);
    tic;
    %形成节点导纳矩阵
    [G,B,transY,transG,transB,transG1,transB1,transNum,lineNum,lineI,lineJ,lineB,transI,transJ,transK0,transKmin,transKmax]=formY(nodeNum,...
        lineNum,mtrLine,mtrGround,mtrTrans,capI,capBi,capK);
else
    %读取IEEE标准格式数据
    [nodeNum,lineNum,capacity,kmax,accuracy,balNum,balNode,mtrLine,mtrGround,mtrTrans,runPara,Pg,Qg,Pd,Qd,pvNum,pvNode,pvV,pvQmin,pvQmax,...
        mtrCapacitor,pgNum,pgNode,a,b,c,Pmin,Pmax,Vmin,Vmax,centrPara]=readData_RPP(data) ;
    tic;
     %形成节点导纳矩阵
    [G,B,lineI,lineJ,lineB,transI,transJ,transK0,transNum,transKmin,transKmax,capI,capGi,capBi,capK,capKmin,capKmax,capNum,transG,transB,...
        transG1,transB1,lineNum]=formY_RPP(nodeNum,lineNum,mtrLine,mtrGround,mtrTrans,mtrCapacitor);
end
%% 初始化
[e,f,PG,QR,Pg,Qg,T,pgvNum,xNum,m,r,l,u,z,w,x,y,para3,mu0,b0,bei,dmaxOut1,dmaxOut2,xz1,xz2,dmaxOut3,dmaxOut4,xz3,xz4,xz5,muOut,...
    opfGoal,times,secondBegin,ifconverge,k2,mismatch]=seting(nodeNum,lineNum,transNum,transK0,capNum,capK,pvQmin,pvQmax,pvNum,pvNode,pvV,...
    pgNum,pgNode,Pmin,Pmax,Pg,Qg);
disp (['迭代次数    互补间隙     最优潮流目标'])
%% 进入循环体求解
for k=0:kmax
    [H1,N1,diagE,diagF,hx,Gap,max_hx,fx,Pg,Qg,xz1,dmaxOut1,xz2,dmaxOut2,xz3,dmaxOut3,xz4,dmaxOut4,PGR]=midResult(nodeNum,pgNode,opfGoal,...
        z,l,u,w,Pg,Qg,e,f,G,B,a,b,c,PG,Pd,Qd,k,xz1,dmaxOut1,xz2,dmaxOut2,xz3,dmaxOut3,xz4,dmaxOut4,times,accuracy,mismatch);
    if Gap<accuracy&&max(abs(hx))<mismatch     %判断互补间隙和最大不平衡量是否达到精度
        if times==1
            times=2;
        end
        if times==2
            if secondBegin==1
                [x,secondBegin,mut,oneTrans,Zt,tij,Tn,kt0,r,l,u,z,w,y,oneCap,capKn,Bci,Zb,mub,q0,Tn1,capKn1,Kn,capK,firstT,firstcapK,...
                    firstfx,firstV,Gap]=snd_setting(x,secondBegin,transNum,T,transKmin,transKmax,r,m,capNum,capK,pgvNum,mu0,capBi,fx,e,f);
            end
            if times==2
                tij_Bci=[tij;Bci];    
                index_0= tij_Bci<0.5;    %获取靠近0的值下标
                index_1=find(tij_Bci>=0.5); %获取靠近1的值下标
                tij_0=max(tij_Bci(index_0));%获取靠近0的值中最大的量
                tij_1=max(ones(size(index_1))-tij_Bci(index_1));%获取靠近1的值中最小的量
                if isempty(tij_0)
                    tmax=tij_1;
                else if isempty(tij_1)
                        tmax=tij_0;
                    else
                        tmax=max(tij_0,tij_1);%获取所有离散变量远离0或1最大的量
                    end
                end
                if (tmax<accuracy)    %判断离散变量是否达到精度
                    ifconverge=1;
                    break
                end
            end
        end
    end
    
    mu=(centrPara*Gap)/(2*r);    %计算扰动因子
    if times==1    %预计算过程
        %形成修正方程系数矩阵
        [HSB,dhx,dgx,dgxT,GIJ,BIJ,L1,U1,LZ,UW,transG,transB]=formHSB(nodeNum,balNode,opfGoal,lineNum,transNum,lineI,lineJ,transI,transJ,T,...
            transG,transB,pgNode,pvNode,G,B,e,f,pgNum,pvNum,pgvNum,xNum,capNum,capI,a,m,r,l,u,z,w,y,N1,H1,diagE,diagF,balNum);
       %形成扰动KKT方程不平衡量矩阵
        [const,Lz,Lw,L1,U1]=Constant(nodeNum,opfGoal,lineNum,transNum,transI,transJ,transG,transB,T,transKmin,transKmax,lineI,lineJ,balNode,...
            mu,G,B,GIJ,BIJ,e,f,Pmin,Pmax,a,b,PG,QR,pvQmin,pvQmax,pgNum,pgvNum,xNum,capNum,capK,capKmin,capKmax,pgNode,pvNode,dhx,dgx,L1,U1,LZ,...
            UW,l,u,y,r,hx,Vmax,Vmin,capBi);
        %求解方程并更新变量
        [z,l,u,w,x,y,PG,QR,Pg,Qg,e,f,T,capK,G,B]=solve(HSB,const,mu,r,nodeNum,pgNum,pvNum,xNum,transNum,capNum,capI,capK,pgNode,pvNode,Pg,...
            Qg,m,x,y,z,l,u,w,Lz,Lw,LZ,UW,dgxT,L1,U1,G,B,transI,transJ,T,transG1,transB1);
    else    %连续化处理过程
        k2=k2+1;
        xz5=[xz5 k2];    %保留每次迭代光滑参数的值，用于作图
        muOut=[muOut mub(2)];
        %形成修正方程系数矩阵
        [HSB,dhx,dgx,dgxT,GIJ,BIJ,L1,U1,LZ,UW,transG,transB,dT_Zt,mut,mub]=snd_formHSB(nodeNum,balNode,opfGoal,lineNum,transNum,lineI,...
            lineJ,transI,transJ,T,transG,transB,pgNode,pvNode,G,B,e,f,pgNum,pvNum,pgvNum,xNum,a,m,r,l,u,z,w,y,kt0,mut,Zt,oneTrans,oneCap,Zb,...
            mub,q0,capBi,capNum,capI,tij,Bci,N1,H1,diagE,diagF,b0,balNum);
        %形成扰动KKT方程不平衡量矩阵
        [const,Lz,Lw,L1,U1]=snd_Constant(nodeNum,opfGoal,lineNum,transNum,transI,transJ,transG,transB,lineI,lineJ,balNode,mu,G,B,GIJ,BIJ,e,f,...
            Pmin,Pmax,a,b,PG,QR,pvQmin,pvQmax,pgNum,pgvNum,xNum,pgNode,pvNode,dhx,dgx,L1,U1,LZ,UW,l,u,y,r,dT_Zt,capNum,hx,Vmax,Vmin,Zt,Zb,bei,...
            mu0,oneTrans,oneCap);
        %求解方程并更新变量
        [z,l,u,w,x,y,PG,QR,Pg,Qg,e,f,T,G,B,Zt,tij,Zb,Bci,capK,mut,mub]=snd_solve(HSB,const,mu,r,nodeNum,pgNum,pvNum,pgvNum,xNum,transNum,pgNode,...
            pvNode,Pg,Qg,m,x,y,z,l,u,w,Lz,Lw,LZ,UW,dgxT,L1,U1,G,B,transI,transJ,T,transG1,transB1,kt0,Tn,mut,oneTrans,q0,capKn,mub,oneCap,capK,...
            capNum,capI,capBi);
    end
end
if k>=kmax
    disp '最优潮流计算不收敛！'
    disp(k)
else
    disp '最优潮流计算收敛！'
    disp(['迭代',num2str(k),'次！'])
end
toc;
%% 输出计算结果
OutResult(e,f,T,k,Gap,fx,transI,transJ,filename,nodeNum,PGR,pgNode,QR,pvNode,capI,capK,firstfx,firstT,firstcapK,Tn,Tn1,capKn,capKn1,ifconverge,...
    max_hx,G,B,balNode,toc,firstV,capBi);

%% 作互补间隙和最大不平衡量收敛曲线图
% plot(xz1,dmaxOut1,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% hold on;
% plot(xz2,dmaxOut2,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% hold on;
% xlabel('迭代次数')
% ylabel('互补间隙')
% plot(xz3,dmaxOut3,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% hold on;
% plot(xz4,dmaxOut4,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% xlabel('迭代次数')
% ylabel('最大不平衡量')

%% 作mut在迭代过程中的变化图
% plot(xz5,muOut,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% xlabel('迭代次数')
% ylabel('μ')
% hold on;

%% 做二次优化的互补间隙图
% plot(xz5,dmaxOut2,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% hold on;
% xlabel('迭代次数')
% ylabel('互补间隙')