%% 21.8.22 RTS-79发输电系统可靠性评估
clear all;
clc;
tic
%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[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;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
num_sp=200000;
mpc=case24_ieee_rts;
num_gen=size(mpc.gen,1);
num_el=size(mpc.branch,1)+num_gen;
m_r=rand(num_el,num_sp);
gen_f0=load('RTS79_gen.dat'); % PMAX;outage_rate
branch_f=load('RTS79_branch.dat');
PD_N=sum(mpc.bus(:,PD)); %额定负荷
VAR=0.01; %失负荷概率方差系数
%% 负荷作为发出负功率的发电机,更改gen和gencost矩阵
idx_gen_load=find(mpc.bus(:,PD));
gen_load=zeros(size(idx_gen_load,1),size(mpc.gen,2));%gen_load
gen_load(:,[GEN_BUS,PG,VG,PMIN])=mpc.bus(idx_gen_load,[BUS_I,PD,VM,PD]);
gen_load(:,[PG,PMIN]) =-gen_load(:,[PG,PMIN]);
gen_load(:,GEN_STATUS)=1;
gen_load(:,MBASE)=100;
gen_load_cost=zeros(size(idx_gen_load,1),size(mpc.gencost,2));%gen_load_cost
gen_load_cost(:,MODEL)=2;
gen_load_cost(:,NCOST)=3;
gen_load_cost(:,COST+1)=10000;%成本一次项系数
mpc.gen=[mpc.gen;gen_load];
mpc.gencost=[mpc.gencost;gen_load_cost];
mpc.bus(:,[PD,QD])=0;
%% 形成发电机和线路的状态转移矩阵
for i=1:num_gen
for j=1:size(gen_f0,1)
if mpc.gen(i,PMAX)==gen_f0(j,1)
gen_f(i,1)=gen_f0(j,2);
end
end
end
s=[gen_f;branch_f];
%% 状态抽样
m_status=ones(num_el,num_sp);
for i=1:num_sp
idx_status=find(m_r(:,i)<s);
if ~isempty(idx_status)
m_status(idx_status,i)=0;
end
end
%% opf
idx_LOLP=0;%失负荷次数
idx_EENS=0;%失负荷量
mpopt = mpoption('verbose',0,'out.all', 0); %verbose——进程信息,默认(1);out.all——打印结果,默认(-1)
for i=1:num_sp
mpc.gen(1:num_gen,GEN_STATUS)=m_status(1:num_gen,i);
mpc.branch(:,BR_STATUS)=m_status(num_gen+1:end,i);
if ~isempty(find(m_status(:,i)==0)) %判断元件是否故障
[RESULTS, SUCCESS]=rundcopf(mpc,mpopt);
idx_PD=abs(sum(RESULTS.gen(num_gen+1:end,PG)));
if PD_N-idx_PD>0.1 %判断是否失负荷
idx_LOLP=idx_LOLP+1; %累计失负荷次数
idx_EENS=idx_EENS+PD_N-idx_PD; %累计失负荷量
end
end
I_LOLP=idx_LOLP/i; %失负荷概率(小时)
I_EENS=idx_EENS/i*8760; %电能不足期望值(MWh/年)
beta=sqrt(I_LOLP*(1-I_LOLP)/i)/I_LOLP; %LOLP 方差系数
if I_LOLP<1 && beta<=VAR %收敛性判断
break;
end
disp(['抽样次数 = ',num2str(i)]);
disp(['方差系数 = ',num2str(beta)]);
disp(['LOLP = ',num2str(I_LOLP)]);
disp(['EENS = ',num2str(I_EENS)]);
end
toc
评论0