【物理应用】基于matlab非序贯蒙特卡洛法评估风电系统【含matlab源码 766期】

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【物理应用】基于matlab非序贯蒙特卡洛法评估风电系统【含matlab源码 766期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、部分源代码

clc
clear
%% 1.计算风速weibull分布
% 数据处理
tic
load data;
mu=mean(speed);%原始数据的统计参数
sigma=sqrt(var(speed));% 计算威布尔分布参数
parmhat=wblfit(speed);
k=parmhat(2);
c=parmhat(1);
% k=(sigma/mu)^-1.086;
% c=mu/gamma(1+1/k);% 威布尔分布拟合
[y,x]=hist(speed,ceil(max(speed)/0.5));%x是区间中心数,组距-1.5
prob1=y/8760/0.5;%计算原始数据概率密度 ，频数除以数据种数，除以组距
prob2=(k/c)*(x/c).^(k-1).*exp(-(x/c).^k);%威布尔分布figure(1)
title('Weibull分布拟合图');bar(x,prob1,1)
hold on
plot(x,prob2,'r')
legend('历史数据','Weibull拟合结果')
% legend('Weibull拟合结果')
hold off
save('result_weibull.mat')
toc
% c=cumsum(prob2)
% plot(x,c)%% 2.ARMA模型预测风速
%
% clc
% clear
% load data
% mu=mean(speed);%原始数据的统计参数
% sigma=sqrt(var(speed));
%
% % 2.1数据标准化
% Stdspeed=(speed-mu)/sigma;
% % Stdspeed=iddata(Stdspeed);
% AIC=zeros(10,1);
% % for n=2:10
% for n=2:7
%     sys=armax(Stdspeed,[n,n-1]);
%     NoiseStd=sqrt(sys.NoiseVariance);
%     e=normrnd(0,NoiseStd,8760,1);%产生白噪声序列[时序标号 对应时序的值]
%     y=zeros(8760,1);%预测后的标准风速[时序标号 对应时序的值]
%
%
% %   2.2利用arma模型产生预测值,1-8760项
%     y(1)=e(1);%第一项等于第一个噪声
%     for i=2:n %2-n项通项公式一致
%         y(i)=-sys.A(2:i)*y(i-1:-1:1)+sys.C(1:i)*e(i:-1:1);
%     end
%
%     for i=n+1:8760% n+1至8760项通项公式一致
%         y(i)=-sys.A(2:n+1)*y(i-1:-1:i-n)+sys.C(1:n)*e(i:-1:i-n+1);
%     end
%
% %    2.3计算AIC值
% % 计算残差
%     s=0;
%     for i=1:8760
%     s=s+(Stdspeed(i)-y(i))^2;
%     end
%
%     AIC(n,1)=8760*log(s)+2*n;
% end
%     [minAIC,optn]=min(AIC(2:7));%找到AIC最小的阶数
%
% %   2.3得到最优阶数，代入arma模型预测风速
% sys=armax(Stdspeed,[optn,optn-1]);
% Noise.Std=sqrt(sys.NoiseVariance);
% e=normrnd(0,Noise.Std,8760,1);%产生白噪声序列[时序标号 对应时序的值]
% y=zeros(8760,1);%预测后的标准风速[时序标号 对应时序的值]
%
%
%
% y(1)=e(1);%第一项等于第一个噪声
% for i=2:optn %2-n项通项公式一致
%     y(i)=-sys.A(2:i)*y(i-1:-1:1)+sys.C(1:i)*e(i:-1:1);
% end
%
% for i=optn+1:8760% n+1至8760项通项公式一致
%     y(i)=-sys.A(2:optn+1)*y(i-1:-1:i-optn)+sys.C(1:optn)*e(i:-1:i-optn+1);
% end
% % 得到预测值y后反变换为实际风速
% Simspeed=y*sigma+mu;
%
% % 计算arma预测风速的概率分布
% [count,x]=hist(Simspeed,ceil(max(Simspeed)/0.5));%x是区间中心数,组距0.5
% prob3=count/8760/0.5;%计算arma预测数据概率密度 ，频数除以数据种数，除以组距
%
% % 时序作图比较
% figure(4)
% plot(1:400,speed(1:400),'r-.',1:400,Simspeed(1:400),'b-');
% legend('实际风速','ARMA模拟风速');
% % 概率分布作图比较
% figure(5)
% bar(x,prob3,1)
% title('ARMA预测概率分布');
%
% save('result_arma.mat')
%% 2.ARMA模型预测风速
clc
clear
tic
load data
y=speed(1:300);
Data=y;              %共300个数据
SourceData=Data(1:250,1); %前250个训练集
step=50;                  %后50个测试
TempData=SourceData;
TempData=detrend(TempData);%去趋势线
TrendData=SourceData-TempData;%趋势函数
%--------差分，平稳化时间序列---------
H=adftest(TempData);
difftime=0;
SaveDiffData=[];
while ~H
SaveDiffData=[SaveDiffData,TempData(1,1)];
TempData=diff(TempData);%差分，平稳化时间序列
difftime=difftime+1;%差分次数
H=adftest(TempData);%adf检验，判断时间序列是否平稳化
end
%---------模型定阶或识别--------------
u = iddata(TempData);
test = [];
for p=1:5                       %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
for q=1:5                    %移动平均对应ACF
m = armax(u,[p q]);
AIC = aic(m);              %armax(p,q),计算AIC
test = [test;p q AIC];
end
end
for k=1:size(test,1)
if test(k,3) == min(test(:,3)) %选择AIC值最小的模型
p_test = test(k,1);
q_test = test(k,2);
break;
end
end
%------1阶预测-----------------
TempData=[TempData;zeros(step,1)];
n=iddata(TempData); %m = armax(u(1:ls),[p_test q_test]);        %armax(p,q),[p_test q_test]对应AIC值最小，自动回归滑动平均模型
m = armax(u,[p_test q_test]);% -------------------------------------------
P1=predict(m,n,1);
PreR=P1.OutputData;
PreR=PreR';Noise.std=sqrt(m.NoiseVariance);
e=normrnd(0,Noise.std,1,300);
for i=251:300PreR(i)=-m.A(2:p_test+1)*PreR(i-1:-1:i-p_test)'+m.C(1:q_test+1)*e(i:-1:i-q_test)';
end
% -------------------------------------------
%----------还原差分-----------------
if size(SaveDiffData,2)~=0
for index=size(SaveDiffData,2):-1:1
PreR=cumsum([SaveDiffData(index),PreR]);
end
end %-------------------预测趋势并返回结果----------------
mp1=polyfit([1:size(TrendData',2)],TrendData',1);
xt=[];
for j=1:step
xt=[xt,size(TrendData',2)+j];
end
TrendResult=polyval(mp1,xt);
PreData=TrendResult+PreR(size(SourceData',2)+1:size(PreR,2));
tempx=[TrendData',TrendResult]+PreR;    % tempx为预测结果
plot(tempx,'r-.');
hold on
plot(Data,'b');
legend('ARMA拟合时序曲线','实际时序风速');
save('resultarma.mat');
toc
%% 2.计及风速和元件故障的风电场出力
% 2.1得到N台机组M状态的风电场出力模型
% 风速单位m/s ,切出功率单位MW
%
clc
clear
load result_weibull
Generator.Wind=struct('vin',3,'vout',25,'vr',15,'Pr',1.5,'FOR',0.028,'lamda',5,'mu',175.2);
N=10;
M=6;
% 风电功率转换函数
[~,k1,k2]=Power([],Generator);StateWeibull=zeros(M,2);
for i=1:MStateWeibull(i,1)=(i-1)/(M-1)*Generator.Wind.Pr;if StateWeibull(i,1)==0
%         StateWeibull(i,2)=1-exp(-(Vin/c)^k)+exp(-(Vout/c)^k))+StateWeibull(i,2)=wblcdf((((2*i-1)/(2*M-2))*Generator.Wind.Pr-k2)/k1,c,k)+exp(-(Generator.Wind.vout/c)^k) ; endif StateWeibull(i,1)>0&&StateWeibull(i,1)<Generator.Wind.PrStateWeibull(i,2)=wblcdf((((2*i-1)/(2*M-2))*Generator.Wind.Pr-k2)/k1,c,k)-wblcdf((((2*i-3)/(2*M-2))*Generator.Wind.Pr-k2)/k1,c,k);endif StateWeibull(i,1)==Generator.Wind.PrStateWeibull(i,2)=wblcdf((Generator.Wind.Pr-k2)/k1,c,k)-wblcdf((((2*i-3)/(2*M-2))*Generator.Wind.Pr-k2)/k1,c,k);end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 门云阁.MATLAB物理计算与可视化[M].清华大学出版社，2013.

【物理应用】基于matlab非序贯蒙特卡洛法评估风电系统【含matlab源码 766期】相关推荐

【Matlab生物电信号】生物电信号仿真【含GUI源码 684期】
一.代码运行视频(哔哩哔哩) [Matlab生物电信号]生物电信号仿真[含GUI源码 684期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1]董兵,超于毅,李 ...
【Matlab语音分析】语音信号分析【含GUI源码 1718期】
一.代码运行视频(哔哩哔哩) [Matlab语音分析]语音信号分析[含GUI源码 1718期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1]韩纪庆,张磊,郑铁 ...
【Matlab身份证识别】身份证号码识别【含GUI源码 014期】
一.代码运行视频(哔哩哔哩) [Matlab身份证识别]身份证号码识别[含GUI源码 014期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1] 蔡利梅.MAT ...
【Matlab水果识别】自助水果超市【含GUI源码 594期】
一.代码运行视频(哔哩哔哩) [Matlab水果识别]自助水果超市[含GUI源码 594期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1]倪云峰,叶健,樊娇娇 ...
【Matlab验证码识别】遗传算法和最大熵优化+大津法（OTSU）+自定义阈值数字验证码识别【含GUI源码 1694期】
一.代码运行视频(哔哩哔哩) [Matlab验证码识别]遗传算法和最大熵优化+大津法(OTSU)+自定义阈值数字验证码识别[含GUI源码 1694期] 二.matlab版本及参考文献 1 matlab ...
【Matlab人脸识别】BP神经网络人脸识别（含识别率）【含GUI源码 891期】
一.代码运行视频(哔哩哔哩) [Matlab人脸识别]BP神经网络人脸识别(含识别率)[含GUI源码 891期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1] ...
【Matlab人脸识别】形态学教室人数统计（带面板）【含GUI源码 1703期】
一.代码运行视频(哔哩哔哩) [Matlab人脸识别]形态学教室人数统计(带面板)[含GUI源码 1703期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1]孟 ...
【Matlab人脸识别】人脸实时检测与跟踪【含GUI源码 673期】
一.代码运行视频(哔哩哔哩) [Matlab人脸识别]人脸实时检测与跟踪[含GUI源码 673期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1]孟逸凡,柳益君 ...
【Matlab图像融合】小波变换遥感图像融合【含GUI源码 744期】
一.代码运行视频(哔哩哔哩) [Matlab图像融合]小波变换遥感图像融合[含GUI源码 744期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1] 包子阳,余 ...
【Matlab语音加密】语音信号加密解密（带面板）【含GUI源码 181期】
一.代码运行视频(哔哩哔哩) [Matlab语音加密]语音信号加密解密(带面板)[含GUI源码 181期] 二.matlab版本及参考文献 1 matlab版本 2014a 2 参考文献 [1]韩纪庆 ...