风速Weibull分布和光伏Beta分布的参数拟合方法

在风光场景生成、随机优化调度等研究中，常常假设风速服从Weibull分布，太阳辐照度服从Beta分布。那我们如何得到两个分布的参数呢？文本首先介绍了风速Weibull分布和辐照度Beta分布的基本概率模型及其性性质，之后以MATLAB代码为例阐述了如何根据历史观测数据对两种分布的参数进行估计。

Weibull分布

风机出力的不确定性主要来源于风速固有的间歇性、随机性和波动性，学术界及工业均普遍认为风速的概率分布服从经典的两参数“Weibull分布”。风速的Weibull分布的概率密度函数(Probability Density Function, PDF)为：
f(v;c,k)=(kc)(vc)k−1exp⁡[−(vc)k]f(v ; c, k)= \left(\frac{k}{c}\right)\left(\frac{v}{c}\right)^{k-1} \exp \left[-\left(\frac{v}{c}\right)^{k}\right] f(v;c,k)=(ck)(cv)k−1exp[−(cv)k]其中，vvv 为实际风速，ccc 是尺度参数（scale parameter），kkk是形状参数（shape parameter）。kkk 是最重要的参数，决定分布密度曲线的基本形状，ccc 起放大或缩小曲线的作用，但不影响分布的形状。

Weibull分布的累积分布函数(Cumulative Distribution Function, CDF)为
F(v;c,k)=1−exp⁡[−(vc)k]F(v; c, k)=1-\exp \left[-\left(\frac{v}{c}\right)^{k}\right] F(v;c,k)=1−exp[−(cv)k]

Beta 分布的均值与方差：

Weibull分布的均值为：

μ=cΓ(1+1k)\mu=c \,\Gamma\left(1+\frac{1}{k}\right) μ=cΓ(1+k1)其中：Γ\GammaΓ 是gamma函数。

Weibull分布的方差为：

σ2=c2[Γ(1+2k)−Γ(1+1k)2]\sigma^2=c^{2}\left[\Gamma\left(1+\frac{2}{k}\right)-\Gamma\left(1+\frac{1}{k}\right)^{2}\right] σ2=c2[Γ(1+k2)−Γ(1+k1)2]

参数估计

两种方法

利用MATLAB自带的 fitdist 函数直接对历史观测数据进行拟合

% 分别对NT个时刻的风速数据进行Weibull分布拟合
% 首先要准备一个风速数据windspeed，维度为(N, NT)
% 其中，N为样本个数，NTime为时刻数如24for i = 1 : NTimewif = fitdist(windspeed(:, i), 'weibull');k_wif(i) = wif.B;c_wif(i) = wif.A;
end

利用随机变量的期望和方差计算参数：包括SDM和MOM两种

% 计算风速均值
wind_speed_avg = mean(wind_speed);
% 计算风速标准差
wind_speed_std = std(wind_speed);% 标准差法参数计算 standard deviation method
k_wstdm = (wind_speed_std/wind_speed_avg)^-1.086;
c_wstdm = (wind_speed_avg)/(gamma(1+(1/k_wstdm)));% 矩量法参数计算 Method of Moments Method Factor Calculation
k_mom = (0.9874/(wind_speed_std/wind_speed_avg))^1.0983;
c_mom = (wind_speed_avg)/(gamma(1+(1/k_mom)));

Beta分布

光伏输出功率主要受太阳辐照度的影响，这与太阳的位置、光伏发电设备的地理位置和天气条件有关。研究表明，一定时间段内的太阳辐照度近似为Beta分布。Beta分布，也称B分布，是一组定义在 (0,1)(0,1)(0,1) 区间的连续型概率密度分布，有两个参数：α,β\alpha, \betaα,β，称为形状参数。Beta 分布的PDF来描述太阳辐照度的概率性质：
f(s;α,β)=Γ(α+β)Γ(α)Γ(β)s(α−1)(1−s)(β−1),0≤s≤1,α≥0,β≥0\begin{aligned} f(s;\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)}s^{(\alpha-1)}(1-s)^{(\beta-1)}, \quad 0 \leq s \leq 1, \alpha \geq 0, \beta \geq 0 \end{aligned} f(s;α,β)=Γ(α)Γ(β)Γ(α+β)s(α−1)(1−s)(β−1),0≤s≤1,α≥0,β≥0
其中 sss 是太阳辐照度。

Beta 分布的均值与方差：

Beta 分布的期望

μ=αα+β\begin{aligned} \mu=\frac{\alpha}{\alpha+\beta} \end{aligned} μ=α+βα

Beta 分布的方差

σ2=αβ(α+β+1)(α+β)2\begin{aligned} \sigma^2 = \frac{\alpha \beta}{(\alpha+\beta+1)(\alpha+\beta)^{2}} \end{aligned} σ2=(α+β+1)(α+β)2αβ

因此，结合上述两个公式可以根据随机变量的历史观测数据的期望和方差求出Beta分布的两个参数：
α=(1−μσ2−1μ)μ2=μ(μ(1−μ)σ2−1)β=α(1μ−1)=(1−μ)(μ(1−μ)σ2−1)\alpha=\left(\frac{1-\mu}{\sigma^2}-\frac{1}{\mu}\right)\mu^2=\mu(\frac{\mu(1-\mu)}{\sigma^2}-1)\\ \beta =\alpha\left(\frac{1}{\mu}-1\right)=(1-\mu)(\frac{\mu(1-\mu)}{\sigma^2}-1) α=(σ21−μ−μ1)μ2=μ(σ2μ(1−μ)−1)β=α(μ1−1)=(1−μ)(σ2μ(1−μ)−1)

参数估计

两种方法

利用MATLAB自带的 fitdist 函数直接对历史观测数据进行拟合

% NTime为时刻数，如24
% NS 为观测数据维度
% solardat为历史观测数据，维度为(NS, NTime)for i = 1 : NTimeparas = fitdist(solardata(:, i), 'beta');a(i) = paras.a;b(i) = paras.b;
end

利用历史观测数据的期望和方差求出Beta的两个参数

% solardata为历史观测数据
% NTime为时刻数，如24
% 计算均值
si_avg = mean(solardata);
% 计算标准差
si_std = std(solardata)^2;a = zeros(1, NTime);
b = zeros(1, NTime);
for i = 1 : NTimea(i) = si_avg(i) * ( si_avg(i) * (1 - si_avg(i))/si_std(i) - 1 );b(i)  = (1 - si_avg(i)) * ( si_avg(i) * (1 - si_avg(i))/si_std(i) - 1);
end