✅作者简介：热爱科研的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)<vci
Pwt_samp(i)=0;
end
if wt_samp(i)>vci&&wt_samp(i)<vN
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.
