function [windspeed_avr active_power_dev windspeed_fin_avr active_power_fin_avr CP]= Powercurve_cal_v03(filename1,filename2,filename3,filename4)
%读取excel的xlsx文件
% B列为功率 C列为风速
data1=xlsread(filename1,1,'B2:C52000');
data2=xlsread(filename2,1,'B2:C52000');
data3=xlsread(filename3,1,'B2:C52000');
data4=xlsread(filename4,1,'B2:C52000');
%将四个Excel表合成一个表格
data=[data1;data2;data3;data4];
% 对数据进行筛选 风速<3m/s 功率<10kw的数据不要
n=length(data(:,2));
for j=1:1:n
if data(j,2)>=3 && data(j,1)>=10
data_opt(j,1)=data(j,2);
data_opt(j,2)=data(j,1);
end
end
%对风速和功率进行排序并返回索引
[windspeed,ind]=sort(data_opt(:,1));
% 求0.5m/s的时间间隔的平均风速和相对应的平均功率,标准差
m=length(ind);
k=0;
u=0;
h=0;
windspeed_sum=0;
active_power_sum=0;
windspeed_avr=0;
active_power_avr=0;
active_power_error_sum=0;
windspeed_fin_sum=0;
active_power_fin_sum=0;
P=1.062;
A=pi*(88/2)^2;
for windspeed_num=3:0.5:19.5
for i=1:1:m
%初步计算平均风速和平均功率
active_power(i)=data_opt(ind(i),2);
if windspeed(i)>=windspeed_num && windspeed(i)
k=k+1;
windspeed_sum= windspeed_sum+windspeed(i);
active_power_sum=active_power_sum+active_power(i);
end
q=(windspeed_num-3)/0.5+1;
windspeed_avr(q)=windspeed_sum/k;
active_power_avr(q)=active_power_sum/k;
%计算标准差
if windspeed(i)>=windspeed_num && windspeed(i)
h=h+1;
active_power_error(i)=(active_power(i)-active_power_avr(q))^2;
active_power_error_sum=active_power_error_sum+active_power_error(i);
end
active_power_dev(q)=sqrt(active_power_error_sum/h);
%剔除散点和平均功率差的绝对值与标准差比值大于4的点,并计算最终的平均风速/功率曲线/发电机功率系数
if windspeed(i)>=windspeed_num && windspeed(i)
e=abs(active_power(i)-active_power_avr(q))/active_power_dev(q);
if e<=4
u=u+1;
windspeed_fin_sum=windspeed_fin_sum+windspeed(i);
active_power_fin_sum=active_power_fin_sum+active_power(i);
end
end
windspeed_fin_avr(q)=windspeed_fin_sum/u;
active_power_fin_avr(q)=active_power_fin_sum/u;
%计算发电机功率系数
CP(q)=1000*active_power_fin_avr(q)/(0.5*P*A*(windspeed_fin_avr(q)^3));
end
windspeed_sum=0;
active_power_sum=0;
active_power_error_sum=0;
windspeed_fin_sum=0;
active_power_fin_sum=0;
k=0;
h=0;
u=0;
end
%画出初步平均风速和平均功率的曲线
figure
plot(windspeed_avr, active_power_avr,'r .-');
set(gca,'xtick',0:1:20)
set(gca,'ytick',0:100:1600)
grid on
xlabel('windspeed_avr m/s')
ylabel('active_power_avr (kW)')
title('power curve')
%画出标准差的曲线
figure
plot(windspeed_avr,active_power_dev,'g .-');
%bar((2*windspeed_num+0.5)/2,active_power_dev);
set(gca,'xtick',0:1:20)
set(gca,'ytick',0:20:200)
grid on
xlabel('windspeed_avr m/s')
ylabel('active_power_dev')
title('power standard deviation')
%画出最终的平均风速和平均功率曲线
figure
plot(windspeed_fin_avr, active_power_fin_avr,'k .-');
set(gca,'xtick',0:1:20)
set(gca,'ytick',0:100:1600)
grid on
xlabel('windspeed_fin_avr m/s')
ylabel('active_power_fin_avr (kW)')
title('final power curve')
%画出最终的发电机功率系数
figure
plot(windspeed_fin_avr,CP,'m .-');
set(gca,'xtick',0:1:20)
set(gca,'ytick',0:0.1:0.6)
grid on
xlabel('windspeed_fin_avr m/s')
ylabel('active_power_dev')
title('CP系数')
%画出采集的散点图
figure
plot(data(:,2),data(:,1),'.');
set(gca,'xtick',0:1:22)
set(gca,'ytick',0:100:1600)
grid on
xlabel('windspeed m/s')
ylabel('active_power')
title('power point ')
![de65960cc196f655ea66fdec7a8b3139.png](https://img-blog.csdnimg.cn/img_convert/de65960cc196f655ea66fdec7a8b3139.png)