matpower切负荷

function [load_shed EXIT] = load_shedding（mpc）
%该函数根据输入的系统数据mpc，计算得到各节点切负荷情况
%mpc：系统数据，其中包括线路、发电机状态
%load_shed：两列的矩阵，第一列为节点编号，第二列为对应切负荷量，单位MW
%EXIT：优化算法收敛标志，如果为1则表明找到最优解

[PQ PV REF NONE BUS_I ~ PD QD GS BS BUS_AREA VM ...
    VA base_KV ZONE VMAX VMIN LAM_P LAM_Q MU_VMAX MU_VMIN] = idx_bus;
[F_BUS T_BUS BR_R BR_X BR_B RATE_A RATE_B ...
    RATE_C TAP SHIFT BR_STATUS PF QF PT QT MU_SF MU_ST ...
    ANGMIN ANGMAX MU_ANGMIN MU_ANGMAX] = idx_brch;
[GEN_BUS PG QG QMAX QMIN VG Mbase GEN_STATUS PMAX PMIN ...
    MU_PMAX MU_PMIN MU_QMAX MU_QMIN ~ PC2 QC1MIN QC1MAX ...
    QC2MIN QC2MAX RAMP_AGC RAMP_10 RAMP_30 RAMP_Q APF] = idx_gen;



%% 形成节点、线路、发电机参数，定义基准功率
bus = mpc.bus;
gen = mpc.gen;
branch = mpc.branch;
baseMVA = mpc.baseMVA;
%修改线路，去掉停运线路
index = find（branch（: BR_STATUS） == 0）;
branch（index :） = [];
%节点数和当前支路数
Nb = size（bus 1）;
Nl = size（branch 1）;

%% 形成B0导纳矩阵
stat = branch（: BR_STATUS）;                    %% ones at in-service branches支路状态
b = stat ./ branch（: BR_X）;                    %% series susceptance电纳
tap = ones（Nl 1）;                              %% default tap ratio = 1
i = find（branch（: TAP））;                       %% indices of non-zero tap ratios （TAP=RATIO）ratio~=0的序号
tap（i） = branch（i TAP）;                        %% assign non-zero tap ratios
b = b ./ tap;
%构建支路节点关联矩阵
f = branch（: F_BUS）;                           %% list of “from“ buses
t = branch（: T_BUS）;