WIND1.zip_wind farm_安全分析_风 潮流_风电 潮流_风电潮流计算

clc; clear; clear all; branchiformation=csvread('bpalinedat14.csv');%读取支路信息 branchnum=size(branchiformation,1);%读取支路数 nodey1=linefy(branchiformation);%生成节点导纳矩阵 nodeinformation=csvread('bpabusdat.csv');%读取节点信息 nodenum1=size(nodeinformation,1);%读取节点数 nodeinf1=csvread('information11.csv');%读取风电节点所服从的分布类型 lineofnodeinf1=size(nodeinf1,1);%显示矩阵行数 lineofnodeinf2=size(nodeinf1,2);%显示矩阵列数 randompower=zeros(lineofnodeinf1,2);%生成行数为风电节点数，列数为2列的零矩阵 randompower(:,1)=nodeinf1(:,1);%将上述矩阵的第一列全部置为要修改潮流的风电节点 %===========================大循环========================================= % cutoutoflimitno=[]; % cutoutoflimitno1=[]; % cutoutoflimitnoa=[]; cutoutoflimityesD=[1,1,1]; cutoutoflimityesB=[1,1,1]; cutoutoflimitCC=[]; countout=[]; countout000=[]; countout1=[]; count222=0; % mod=[]; % mod13=[]; % mod4=[]; % cutmod=[]; % cutmod13=[]; % cutmod4=[]; totalnum=500; for totalcycle1=1:totalnum count111=0; [randompower]=randompowerRP(lineofnodeinf1,nodeinf1,randompower); %============================一下调用modify4.exe修改随机注入功率下的dat文件================ nodeData1=randompower%nodeData1为传递函数传递节点数 initialBpadatfile='bpa14.dat'; modifiedBpaDatfile=strcat('mod',initialBpadatfile); modifyNodePower(nodeData1,initialBpadatfile,modifiedBpaDatfile);%调用midify4修改节点注入功率 pause(1) %%%%修改！ %==========================以下调用mdify15.exe计算一下本次注入功率状态下的潮流============== modifiedBpaDatfile1='modbpa14.dat'; [calcheck11,voltageresult1]=flowcal1(modifiedBpaDatfile1,nodenum1); if calcheck11==300 end; modvoltagephaseangle=load('readmodbpa14.txt');%返回电压相角结果 % mod=[mod;modvoltagephaseangle]; % mod13=[mod13;modvoltagephaseangle(13,1)]; % mod4=[mod4;modvoltagephaseangle(13,1)]; %[mods2]=PowerFlow1(modvoltagephaseangle,branchnum,branchiformation)%返回潮流结果 %=============以下循环进行支路开断模拟===================================================== modifiedBpaDatfile2='modbpa14.dat'; cutoutoflimit1=[]; cutoutoflimit1B=[]; cutoutoflimit1C=[]; for counttemp1=1:branchnum cutdat1=linecut(modifiedBpaDatfile2,counttemp1); [calcheck11,voltageresult1,countAAA]=flowcal2(cutdat1,nodenum1,countout,counttemp1); if calcheck11==300 cutmodvoltagephaseangle=load('readcutmodbpa14.txt');%将生成的开断支路潮流结果读入cutmodvoltagephaseangle中 cutmodvoltage=cutmodvoltagephaseangle(:,1);%的开断支路潮流电压模值结果 cutmodphaseangle=cutmodvoltagephaseangle(:,2); % cutmod=[cutmod;cutmodvoltagephaseangle];%的开断支路潮流相角结果 % cutmod13=[cutmod13;cutmodvoltagephaseangle(13,1)]; % cutmod4=[cutmod4;cutmodvoltagephaseangle(4,1)]; modvoltage=modvoltagephaseangle(:,1); %[cutp2,cutq2,cutoutoflimit2,cutoutoflimit2B,cutoutoflimit2C,count]=PowerFlow2(cutmodvoltagephaseangle,branchnum,mods2,branchiformation,modvoltage,cutmodvoltage,cutmodphaseangle,nodenum1); [cutoutoflimit2C,count]=PowerFlow22(branchnum,modvoltage,cutmodvoltage,nodenum1); % cutoutoflimit1=[cutoutoflimit1;cutoutoflimit2];%视载功率越限的支路号 % cutoutoflimit1B=[cutoutoflimit1B;cutoutoflimit2B];%电压相位差越限 cutoutoflimit1C=[cutoutoflimit1C;cutoutoflimit2C];%电压模值越限 generatorlimited=csvread('generatorlimited.csv'); generatornum=size(generatorlimited,1); %Inequalityconstraints3(generatornum,generatorlimited,cutp2,cutq2) countout=[countout;countAAA]; count111=count111+count; end end % disp('系统中支路潮流超过视在功率最大允许限值：'); % [cutoutoflimitD]=Limitindicators(cutoutoflimit1) % disp('相邻节点间的电压相位差超过最大允许偏差：'); % [cutoutoflimitB]=Limitindicators(cutoutoflimit1B) disp('节点电压模值超过了所允许的上下限：'); [cutoutoflimitC]=Limitindicators1(cutoutoflimit1C) % [cutoutoflimitYD]=Limitindicators2(cutoutoflimitD); % if isempty(cutoutoflimitYD) % cutoutoflimitYD=[1,1,1]; % end % [cutoutoflimitYB]=Limitindicators2(cutoutoflimitB); % if isempty(cutoutoflimitYB) % cutoutoflimitYB=[1,1,1]; % end % cutoutoflimityesD=[cutoutoflimityesD;cutoutoflimitYD]; % cutoutoflimityesB=[cutoutoflimityesB;cutoutoflimitYB]; % cutoutoflimitO=[cutoutoflimitD;cutoutoflimitB]; % cutoutoflimityes=cutoutoflimitO%[cutoutoflimityes]=Limitindicators2(cutoutoflimitO); % cutyessize=size(cutoutoflimityes,1); % cutnosize=size(cutoutoflimitno,1); % if (isempty(cutoutoflimityes))&&(isempty(cutoutoflimitno)) % else % if cutyessize>=cutnosize; % cutoutoflimitno(cutyessize,3)=0; % else % if cutyessize<cutnosize; % cutoutoflimityes(cutnosize,3)=0; % end % end % end % cutoutoflimitno1=[cutoutoflimitno;cutoutoflimityes]; % cutoutoflimitno1(all(cutoutoflimitno1==0,2),:)=[] ; % [cutoutoflimitnoa]=Limitindicators2(cutoutoflimitno1); cutoutoflimitCC=[cutoutoflimitCC;cutoutoflimitC]; countout000=[countout000;countout]; % disp('支路开端后拥有的最小奇异值为min，此支路为：') % min=min(ssss) % [number1,number2]=find(ssss==min); % [branchiformation(number2,1),branchiformation(number2,2)] count222=count222+count111; end % disp('支路开断后进行潮流计算结果不收敛，此支路为：') % [cutoutoflimit]=Limitindicators3(countout000,totalnum,branchnum) % disp('发生支路越限的支路及越限次数以及概率为：') % [cutoutoflimitbranch]=Limitindicators3(cutoutoflimitnoa,totalnum,branchnum,count222) % disp('发生越限的节点及越限次数以及概率为：') % [cutoutoflimitnum]=Limitindicators4(cutoutoflimitCC,totalnum,branchnum) % disp('发生支路越限的支路及越限次数以及概率为：') % [cutoutoflimitbranch]=Limitindicators3(cutoutoflimityesD,totalnum,branchnum,count222) % disp('发生支路越限的支路及相角越限次数以及概率为：') % [cutoutoflimitbranch]=Limitindicators3(cutoutoflimityesB,totalnum,branchnum,count222) disp('发生越限的节点及越限次数以及概率为：') [cutoutoflimitnum]=Limitindicators4(cutoutoflimitCC,totalnum,branchnum) %AAAA(mod,cutmod);