DFIG控制3： 风机模型和最基础的MPPT

DFIG3控制3： 风机模型和最基础的MPPT
基于教程的第3部分。（第二段对比了仿真和根据模型计算的稳态值，以后再学）。
DFIM Tutorial 3 – Wind Turbine Model based on Doubly Fed Induction Generator in MATLAB-Simulink

这个教程提供的一套基本的风机模型参数，还是挺有帮助的。

风机的一些公式

记录一下风机的几个基本公式。
风机的转矩$T_t$和风速$V_w$的关系：$T_t=\frac{1}{2}\rho \pi R^3{V}_w^2{C}_t$
风机的风能利用系数$C_p$和转矩系数$C_t$的关系：$C_p(\lambda )=\lambda C_t(\lambda )$
风机的风能利用系数$C_p(\lambda ,\beta )$常用表达式：
$$\begin{array}{rl}{C}_p(\lambda ,\beta )& =c_1(\frac{c_2}{{\lambda }_i}-c_3\beta -c_4{\beta }^5-c_6)e^{c_7/{\lambda }_i}\\ \frac{1}{{\lambda }_i}& =\frac{1}{\lambda +c_8\beta }-\frac{c_9}{{\beta }^3+1}\end{array}$$
其中，
叶尖速比，$\lambda =\frac{R{\Omega }_t}{V_w}$
桨距角 $\beta$（在这个仿真中始终为0）
C1-C9是一些系数。
基于上面的公式，可以画出$C_t$和$\lambda$曲线如下图（左），可见有一个最佳的$C_t$值（前提是风速在一定范围内），也就是说，当风速变化时，可以通过调节风机转速（也就是叶尖速比$\lambda =\frac{R{\Omega }_t}{V_w}$）来达到最优$C_t$值，实现最大功率跟踪MPPT（条件是风速在一定范围内）。

右图表示风机输出的功率和风速的关系：

1. 风速过小时，风机不转，无功率输出
2. 风速适中时，通过MPPT达到当前风速的最大功率，随着风速的增大，输出功率上升
3. 风速高于设定值，通过控制桨距角 $\beta$等，输出功率不再上升

MPPT

在最大功率点，$C_t=C_{topt}$，$C_p=C_{pmax}$，对应于最佳的叶尖速比 $\lambda ={\lambda }_{opt}$
因此有：
$$\begin{array}{rl}{T}_{topt}& =\frac{1}{2}\rho \pi R^3{V}_w^2{C}_{topt}\\ & =\frac{1}{2}\rho \pi R^3\frac{R^2{\Omega }_t^2}{{\lambda }_{opt}^2}\frac{C_{pmax}}{{\lambda }_{opt}}\\ & =\frac{1}{2}\rho \pi R^5\frac{C_{pmax}}{{\lambda }_{opt}^3}{\Omega }_t^2\\ & =k_{opt}{\Omega }_t^2\end{array}$$
而这个系数可以直接计算得到（其中${\lambda }_{opt}$和$C_{pmax}$通过实际测量获取）：$k_{opt}=\frac{1}{2}\rho \pi R^5\frac{C_{pmax}}{{\lambda }_{opt}^3}$
间接转速控制，就是根据风机转速，直接计算最优的转矩，然后换算后得到DFIG的q轴电流参考值，然后去控制DFIG的电流： ${\Omega }_t\to T_{topt}\to i_q$
如下图：

变速箱

风机（低速，高转矩）<—变速箱1:N—>DFIG（高速，低转矩）
这个仿真中，N=100。
DFIG（下标m）和风机（下标t）的转速和转矩关系：
$$\begin{array}{rl}{\omega }_m& =N{\Omega }_t\\ T_m& =\frac{1}{N}{T}_t\end{array}$$
这和后面的一些计算有关。

仿真模型

顶层

1. 转速由MPPT控制，所以转速PI控制不再生效
2. 新增风机和变速器模型，输入风速和DFIG转速omega_m，输出DFIG的转矩Tm。（上一个仿真中，转矩Tm是常数）

风机模型

1. 计算$\lambda$
2. 提供$C_t-\lambda$关系，使用查找表和插值法，输入$\lambda$，输出$C_t$
3. 根据上文的公式，计算风机的转矩和DFIG的转矩。
   

控制器

1. damping，Dt_m，影响很小，仿真中没添加。
2. 其实就是上文的$k_{opt}$相关公式，注意变速箱对转速和转矩的影响。
   

控制系统中，除了不使用速度PI控制，其他不变。

仿真结果

和视频一样，选两个风速做仿真。

根据功率的公式，$P=T_e{\omega }_m$

1. 风速7m/s，转速约120rad/s，转矩约4250N・m，功率约0.5MW
2. 风速8.5m/s，转速约145rad/s，转矩约6250N・m，功率约0.91MW

与之前的功率-风速曲线对应，说明的确基本实现了MPPT：

初始化代码

基于上一个仿真，新增了风机部分，初始化代码如下：

clear
% DFIG parameters -> Rotor parameters referred to the stator side
f = 50; % Stator frequency (Hz)
Ps = 2e6; % Rated stator power (W)
n = 1500; % Rated rotational speed (rev/min)
Vs = 690; % Rated stator voltage (V)
Is = 1760; % Rated stator current (A)
Tem = 12732; % Rated torque (N.m)

p=2; % Pole pair

u = 1/3; % Stator/rotor turns ratio
Vr = 2070;    % Rated rotor voltage (non-reached) (V) 
smax = 1/3;   % Maximum slip
Vr_stator = (Vr*smax) *u;    % Rated rotor voltage referred to stator (V)
Rs = 2.6e-3;    % Stator resistance(ohm)                
Lsi = 0.087e-3; %  Leakage inductance (stator & rotor) (H)                 
Lm = 2.5e-3; %  Magnetizing inductance (H)
Rr = 2.9e-3; % Rotor resistance referred to stator (ohm)
Ls = Lm + Lsi; % Stator inductance (H)
Lr = Lm + Lsi; % Rotor inductance (H)
Vbus = Vr_stator*sqrt(2); % DC de bus voltge referred to stator (V)

sigma = 1-Lm^2/(Ls*Lr);

Fs = Vs*sqrt(2/3)/(2*pi*f);  % Stator Flux (aprox.) (Wb)
J = 127/2; % Inertia, originally 127, reduced by 2 to make the response faster
D = 1e-3; %damping 

fsw = 4e3; 
Ts =  1/fsw/50;

kT = -1.5*p*(Lm/Ls)*Fs; % kT, coef of output of the speed controller  

tau_i = (sigma*Lr)/Rr;
tau_n = 0.05; 
wni = 100*(1/tau_i);
wnn = 1/tau_n;

kp_id = (2*wni*sigma*Lr)-Rr;
kp_iq = kp_id;

ki_id = (wni^2)*Lr*sigma;
ki_iq = ki_id;

kp_n = (2*wnn*J)/p; %kp_n = (2*wnn*J)/p;
ki_n = (wnn^2)*J/p; %ki_n = (wnn^2)*J/p;


% Three blade wind turbine mode
N = 100; % Gearbox ratio
Radio= 42; % blade radius 
ro= 1.225; % Air density

% Cp and Ct curves
beta=0;  % Pitch angle
ind2=1;

for lambda=0.1:0.01:11.8
    lambdai(ind2)= (1./((1./(lambda-0.02.*beta) + (0.003./ (beta^3+1)))));
    Cp(ind2)=0.73*(151./lambdai(ind2) - 0.58*beta - 0.002.*beta^2.14-13.2).*(exp(-18.4./lambdai (ind2))); 
    Ct(ind2)=Cp(ind2)/lambda;
    ind2=ind2+1;
end

tab_lambda=[0.1:0.01:11.8];


% Kopt for MPPT
Cp_max = 0.44; 
lambda_opt = 7.2;
Kopt = ((0.5*ro*pi* (Radio^5) *Cp_max)/(lambda_opt^3));
% Power curve in fucntion of wind speed

P = 1.0e+06 *[0,0,0,0,0,0,0,0.0472,0.1097,0.1815,0.2568,0.3418, ...
            0.4437,0.5642, 0.7046, 0.8667,1.0518,1.2616, 1.4976, 1.7613,2.0534,...
            2.3513,2.4024,2.4024,2.4024, 2.4024,2.4024,2.4024];
V = [0.0000,0.5556,1.1111,1.6667,2.2222,2.7778,3.3333,3.8889, 4.4444,... 
    5.0000,5.5556,6.1111,6.6667,7.2222,7.7778,8.3333,8.8889,9.4444, ...
    10.0000,10.5556,11.1111, 11.6667,12.2222,12.7778,13.3333, 13.8889,...
    14.4444,15.0000];
figure
subplot(1,2,1)
plot(tab_lambda, Ct, 'linewidth',1.5)
xlabel('\lambda', 'fontsize',14)
ylabel('Ct', 'fontsize',14)
subplot(1,2,2)
plot (V, P, 'linewidth', 1.5)
grid
xlabel('Wind speed (m/s)', 'fontsize',14)
ylabel('Power (W)', 'fontsize',14)

一些问题

1. 对于风机，暂时只了解了这几个公式
2. 这个MPPT是直接用公式算了一下理论的最佳转矩。还是要再看下相关的公式和模型。和其他的MPPT策略。