【电力系统】基于Matlab实现风电光伏概率潮流计算

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法  神经网络预测 雷达通信  无线传感器

信号处理 图像处理 路径规划 元胞自动机 无人机  电力系统

⛄ 内容介绍

潮流计算是电力系统的一项重要分析功能,是进行故障计算,继电保护整定,安全分析的必要工具.本文提出了利用Matlab语言来进行电力系统潮流计算.通过算例,说明了该方法编程简便,运算效率高并符合人们的思维习惯,计算结果能满足工程计算需要,同时验证了该方法的有效性.

⛄ 部分代码

%相关原始数据格式说明如下：

%n——节点个数；n1——支路条数；isb——平衡节点号；H——PQ节点个数（为后面形成PVU存储PV节点初始电压用）；

%pr——误差精度。

%B1——支路参数矩阵，其中第一列和第二列是起始节点编号和终点节点编号,第三列、第四列、第五列、第六列分别为：

%支路电阻、电抗、变压器变比、电纳。（不考虑电导）

%B2——节点参数矩阵，其中第一列和第二列为节点编号和节点类型；第三列到第六列分别为：

%注入有功、注入无功、电压幅值、电压相位。

%节点类型分类如下：“0”为平衡节点,“1”为PQ，“2”为PV节点；“3”为PQ（V）节点，“4”为PI节点。

clc;

clear;

close all

Ploss_before=239;

V_before=[1.05000000000000;1.04539984601725;1.03826381242661;1.02656224417193;...

    1.01497413857938;0.986061941470649;0.980461558681093;0.972856430301411;...

    0.962880928647019;0.953612848625365;0.952264312724239;0.949911247256376;...

    0.940258149608386;0.936656329416138;0.934414308316597;0.932247854435629;...

    0.929012818991491;0.928050797399079;1.04459070499741;1.03910777072861;...

    1.03802728200737;1.03704915640893;1.03278914604802;1.02259101804636;...

    1.01750294074663;0.983041861628259;0.979024678756850;0.961078690739743;...

    0.948179211559331;0.942585582449513;0.936049346990836;0.934610825106120;0.934165000478165];

B2=[1 0 0 0 1.05 0;

    2 1 -0.01 -0.006 1 0;

    3 1 -0.009 -0.004 1 0;

    4 1 -0.012 -0.008 1 0;

    5 1 -0.006 -0.003 1 0;

    6 1 -0.006 -0.002 1 0;

    7 1 -0.02 -0.01 1 0;

    8 1 -0.02 -0.01 1 0;

    9 1 -0.006 -0.002 1 0;

    10 1 -0.006 -0.0035 1 0;

    11 1 -0.0045 -0.003 1 0;

    12 1 -0.006 -0.0035 1 0;

    13 1 -0.006 -0.0035 1 0;

    14 1 -0.012 -0.008 1 0;

    15 1 -0.006 -0.001 1 0;

    16 1 -0.006 -0.002 1 0;

    17 1 -0.006 -0.002 1 0;

    18 1 -0.009 -0.004 1 0;

    19 1 -0.009 -0.004 1 0;

    20 1 -0.009 -0.004 1 0;

    21 1 -0.009 -0.004 1 0;

    22 1 -0.009 -0.004 1 0;

    23 1 -0.009 -0.005 1 0;

    24 1 -0.042 -0.02 1 0;

    25 1 -0.042 -0.02 1 0;

    26 1 -0.006 -0.0025 1 0;

    27 1 -0.006 -0.0025 1 0;

    28 1 -0.006 -0.002 1 0;

    29 1 -0.012 -0.007 1 0;

    30 1 -0.02 -0.06 1 0;

    31 1 -0.015 -0.007 1 0;

    32 1 -0.021 -0.01 1 0;

    33 1 -0.006 -0.004 1 0];

%% 蒙特卡洛随机潮流计算样本规模

times=5e3;

%% 蒙特卡洛模拟采样随机变量样本

Pi=abs(B2(:,3));

Qi=abs(B2(:,4));

nPQi=length(Pi);

Nodenum=33;

% 负荷服从正态分布，变异系数ld_ero

ld_ero=0.05;

Pld_samp=zeros(Nodenum,times);

Qld_samp=zeros(Nodenum,times);

% 负荷有功、无功出力样本

for k=1:nPQi

    Pld_samp(k,:)=normrnd(Pi(k),Pi(k)*ld_ero,1,times);             

    Qld_samp(k,:)=normrnd(Qi(k),Qi(k)*ld_ero,1,times);

end

% 光伏有功服从Beta分布，单位功率因数

Ppv_samp=zeros(1,times);

% Beta分布的两个形状参数

a_pv=2.06;  %0.45

b_pv=2.5;  %9.18

% 光伏发电相关参数：组件总面积S_pv、光电转换率prey_pv、最大光强rmax(MW/m2)

S_pv=10000;

prey_pv=0.14;

rmax=0.8e-3;

% 光伏有功出力样本

pv_samp(1,:)=betarnd(a_pv,b_pv,1,times);

 Ppv_samp(1,:)=pv_samp(1,:)*rmax*S_pv*prey_pv*1000;

 % weibull分布的两个形状参数

k_wt=2.25;

c_wt=4.5;

wt_samp =wblrnd(c_wt,k_wt,1,times);

PN_wt=500;

    vci=3 ;

    vN=10.5;

    vco=22 ;

    for i=1:times

        if wt_samp(i)

            Pwt_samp(i)=0;

        end

         if wt_samp(i)>vci&&wt_samp(i)

            Pwt_samp(i)=(wt_samp(i)-vci)/(vN-wt_samp(i))*PN_wt;

            if   Pwt_samp(i)>PN_wt

                 Pwt_samp(i)=PN_wt;

            end

         end

         if wt_samp(i)>vN&&wt_samp(i)< vco

            Pwt_samp(i)=PN_wt;

         end

         if wt_samp(i)> vco

            Pwt_samp(i)=0;

         end

    end

    

    

   j=1; 

    for i=1:times

plot(centers, counts / sum(counts))%画直方图，x是各个区间中间值，y是对应的概率（区间数/总数量）

xlabel('光照强度');

ylabel('概率');

title('光照强度概率分布图')

[counts,centers] = hist(wt_samp , 100);%返回指定数据，特定间隔（如7）的各个分割区间的数量和中间值% 画频率分布直方图

figure(2)

plot(centers, counts / sum(counts))%画直方图，x是各个区间中间值，y是对应的概率（区间数/总数量）

 xlabel('风速');

ylabel('概率');

title('风速概率分布图')

  [counts,centers] = hist(Ploss, 100);%返回指定数据，特定间隔（如7）的各个分割区间的数量和中间值% 画频率分布直方图

figure(3)

plot(centers, counts / sum(counts))%画直方图，x是各个区间中间值，y是对应的概率（区间数/总数量）

 xlabel('网损');

ylabel('概率');

title('网损概率分布图')

flage=5;  %选取不同节点得到节点电压概率分布图

 V_flage=V_amp(:,flage);

  [counts,centers] = hist(V_flage, 100);%返回指定数据，特定间隔（如7）的各个分割区间的数量和中间值% 画频率分布直方图

figure(4)

plot(centers, counts / sum(counts))%画直方图，x是各个区间中间值，y是对应的概率（区间数/总数量）

 xlabel('电压幅值');

ylabel('概率');

title('电压幅值概率分布图')

% ymin=min( V_flage); 

% ymax=max( V_flage); 

% x=linspace(ymin,ymax,100); %将最大最小区间分成20个等分点(19等分),然后分别计算各个区间的个数

% yy=hist( V_flage,x); %计算各个区间的个数

% yy=yy/length( V_flage); %计算各个区间样本数所占总样本数的比例（概率）

% bar(x,yy) %画出概率密度分布图

% 

% s=0; 

% for i=2:length(x) 

% s=[s,trapz(x([1:i]),yy([1:i]))]; % 使用trapz模拟积分 CDF

% end

% figure;

% plot(x,s,'*-')

flage=25;  %选取不同节点得到节点电压概率分布图

 V_flage=V_amp(:,flage);

  [counts,centers] = hist(V_flage, 100);%返回指定数据，特定间隔（如7）的各个分割区间的数量和中间值% 画频率分布直方图

figure(7)

plot(centers, counts / sum(counts))%画直方图，x是各个区间中间值，y是对应的概率（区间数/总数量）

 xlabel('电压幅值');

ylabel('概率');

title('电压幅值概率分布图')

 V_aver=mean(V_amp,1);  %均值

  Ploss_aver=mean(Ploss,2);

 d_v=std(V_amp,0,1);  %方差

  figure(5)

plot(V_before);

hold on

plot( V_aver);

legend('接入前','接入后')

xlabel('节点序号');

ylabel('电压幅值');

⛄ 运行结果

⛄ 参考文献

[1]徐劲松, 宁玉琳, 杨永锋. 基于Matlab的电力系统PQ分解法潮流计算研究[J]. 电气传动自动化, 2011, 33(2):9.

[2]郗忠梅, 李有安, 赵法起,等. 基于Matlab的电力系统潮流计算[J]. 山东农业大学学报：自然科学版, 2010(2):4.

⛄ Matlab代码关注

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料