ncp_sigmoid/Sigmoid函数法/OPF.m

%*******************************************************
%            <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵͳ<EFBFBD>޹<EFBFBD><EFBFBD>Ż<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%            <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%            <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>㷨<EFBFBD><EFBFBD><EFBFBD>ִ<EFBFBD><EFBFBD>ڵ㷨<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Sigmoid<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɢ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%            Ŀ<EFBFBD>꺯<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD>1<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>й<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD>2<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>й<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD>3<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>޹<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD>4<EFBFBD><EFBFBD>
%            <EFBFBD><EFBFBD><EFBFBD>ߣ<EFBFBD>    Τ԰<EFBFBD><EFBFBD>
%            <EFBFBD><EFBFBD>дʱ<EFBFBD>䣺2010<EFBFBD><EFBFBD>10<EFBFBD><EFBFBD>
%*******************************************************
function []=OPF(filename)
clc;
%% <EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
data=dlmread(filename);
if data(1,1)<1780
    %<EFBFBD><EFBFBD>ȡ1780<EFBFBD>ڵ<EFBFBD>ϵͳ<EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD>͵<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    [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;
    %<EFBFBD>γɽڵ㵼<EFBFBD>ɾ<EFBFBD><EFBFBD><EFBFBD>
    [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
    %<EFBFBD><EFBFBD>ȡIEEE<EFBFBD><EFBFBD>׼<EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    [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;
     %<EFBFBD>γɽڵ㵼<EFBFBD>ɾ<EFBFBD><EFBFBD><EFBFBD>
    [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
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD>
[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 (['<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>    <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>϶     <EFBFBD><EFBFBD><EFBFBD>ų<EFBFBD><EFBFBD><EFBFBD>Ŀ<EFBFBD><EFBFBD>'])
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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     %<EFBFBD>жϻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>϶<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD>ﵽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        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;    %<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>ֵ<EFBFBD>±<EFBFBD>
                index_1=find(tij_Bci>=0.5); %<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>1<EFBFBD><EFBFBD>ֵ<EFBFBD>±<EFBFBD>
                tij_0=max(tij_Bci(index_0));%<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
                tij_1=max(ones(size(index_1))-tij_Bci(index_1));%<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>1<EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
                if isempty(tij_0)
                    tmax=tij_1;
                else if isempty(tij_1)
                        tmax=tij_0;
                    else
                        tmax=max(tij_0,tij_1);%<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɢ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Զ<EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
                    end
                end
                if (tmax<accuracy)    %<EFBFBD>ж<EFBFBD><EFBFBD><EFBFBD>ɢ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD>ﵽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
                    ifconverge=1;
                    break
                end
            end
        end
    end
    
    mu=(centrPara*Gap)/(2*r);    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ŷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    if times==1    %Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        %<EFBFBD>γ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        [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);
       %<EFBFBD>γ<EFBFBD><EFBFBD>Ŷ<EFBFBD>KKT<EFBFBD><EFBFBD><EFBFBD>̲<EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        [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);
        %<EFBFBD><EFBFBD><EFBFBD>ⷽ<EFBFBD>̲<EFBFBD><EFBFBD><EFBFBD><EFBFBD>±<EFBFBD><EFBFBD><EFBFBD>
        [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    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        k2=k2+1;
        xz5=[xz5 k2];    %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ÿ<EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD>⻬<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ
        muOut=[muOut mub(2)];
        %<EFBFBD>γ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        [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);
        %<EFBFBD>γ<EFBFBD><EFBFBD>Ŷ<EFBFBD>KKT<EFBFBD><EFBFBD><EFBFBD>̲<EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
        [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);
        %<EFBFBD><EFBFBD><EFBFBD>ⷽ<EFBFBD>̲<EFBFBD><EFBFBD><EFBFBD><EFBFBD>±<EFBFBD><EFBFBD><EFBFBD>
        [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 '<EFBFBD><EFBFBD><EFBFBD>ų<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>㲻<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>'
    disp(k)
else
    disp '<EFBFBD><EFBFBD><EFBFBD>ų<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>'
    disp(['<EFBFBD><EFBFBD><EFBFBD><EFBFBD>',num2str(k),'<EFBFBD>Σ<EFBFBD>'])
end
toc;
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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);

%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>϶<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ
% 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('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>')
% ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>϶')
% 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('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>')
% ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>')

%% <EFBFBD><EFBFBD>mut<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>еı仯ͼ
% plot(xz5,muOut,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% xlabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>')
% ylabel('<EFBFBD><EFBFBD>')
% hold on;

%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ż<EFBFBD><EFBFBD>Ļ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>϶ͼ
% plot(xz5,dmaxOut2,'-ks','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',3)
% hold on;
% xlabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>')
% ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>϶')