[发明专利]一种风力发电机组全工况模型预测控制方法及系统

权利要求书

1.一种风力发电机组全工况模型预测控制系统，包括风轮、传动链、塔筒、发电单元、变桨驱动器、变流器，其特征在于：该控制系统还包括反馈信息测量器和MPC控制器；

发电单元的输出端接至电网，变桨驱动器输入端与MPC控制器输出端口相连，变流器输入端连接到MPC的输出端口，变流器输出端连接到发电单元；

所述反馈信息测量器包括安装于风力发电机主轴的转速传感器、安装于滑环尾部的风轮绝对位置编码器、安装于风力发电机尾部的发电机编码器、安装于机舱底盘的振动分析仪、安装于发电机定子侧的电压测量传感器、安装于发电机定子侧的电流测量传感器、装设与机舱外部支架上的风速传感器和安装于变距减速机旁的叶片桨距角编码器，各反馈信息测量器的输出端与MPC控制器的输入端口相连；

所述MPC控制器选用PLC作为硬件系统，PLC包括机舱子站和塔底子站，机舱子站和塔底子站之间通过FAST总线通信。

2.采用如权利要求1所述的风力发电机组全工况模型预测控制系统进行风力发电机组全工况模型预测控制的方法，其特征在于：具体按如下步骤进行：

步骤1：系统的离线准备；

MPC控制器的离线准备工作，包括：理论计算：状态矩阵A，控制输入增益矩阵B，扰动输入增益矩阵Γ；选择控制目标，确定测量状态输出矩阵C_y，为被控状态输出矩阵C_z；根据机型确定被控输入和被控输出的约束条件、被控输出设定值z_r随风速变化曲线；选择合适的采样周期T、预测时域N、控制时域N_c、参考轨迹时间常数T_g、校正向量h、输出误差权重Q、控制量变化幅度权重R；将以上数据存入存储模块，以便实时调用；

步骤2：反馈信息测量器实时检测风力发电机组的状态变量，通过风力发电机组全工况预测模型计算控制增量；

风力发电机组的状态变量包括：风轮转速、传动链扭转位移、传动链扭转速度、叶片平面外一阶挥舞位移、叶片平面外一阶挥舞速度、塔筒前后一阶摆动位移和塔筒前后一阶摆动速度；

计算控制增量采用两个预测模型：变桨控制预测模型和转矩控制预测模型，在不同工况下进行自动切换，在额定风速和额定风速以下时采用转矩控制预测模型，控制量为转矩，状态变量包括：风轮转速、传动链扭转位移、传动链扭转速度、叶片平面外一阶挥舞位移、叶片平面外一阶挥舞速度、塔筒前后一阶摆动位移和塔筒前后一阶摆动速度，并将风速增量作为扰动输入变量；在额定风速以上时采用变桨控制预测模型，控制量为桨距角，状态变量包括：风轮转速、传动链扭转位移和传动链扭转速度，并将风速增量作为扰动输入变量；

所述变桨控制预测模型，在离散化之前的状态空间模型如下：

x·1x·2x·3x·4x·5x·6x·7=Irot0003M400M46010000000103M40Irot0M46Irot0001000M04000M000M060000010M640003M600M66-1]]>

·γ-Kd-Cd0-3C400-C460010000γIrot-(KdIgen+KdIrot)-(CdIgen+CdIrot)0-3C40Irot0-C46Irot0000100-C0400-K00-C00000000001-C6400-K60-C60-K60-C66x1x2x3x4x5x6x7]]>

+Irot0003M400M46010000000103M40Irot0M46Irot0001000M04000M000M060000010M640003M600M66-13ζ103ζ1Irot0ζb0ζtδβ]]>

+Irot0003M400M46010000000103M40Irot0M46Irot0001000M04000M000M060000010M640003M600M66-13α103α1Irot0αb0αtδν---(1)]]>

y=C·x1x2x3x4x5x6x7]]>

所述转矩控制预测模型，离散化之前的状态空间模型如下：

x·1x·2x·3γIrot-KdIror-CdIrot001γIrot-(KdIgen+KdIrot)-(CdIgen+CdIrot)x1x2x3+001IgenδTgen+αIrot0αIrotδν---(2)]]>

y=C·x1x2x3]]>

变桨控制预测模型和转矩控制预测模型的统一状态空间标准形式为：

x·=Ax+BΔu+ΓΔυ---(3)]]>

y＝Cx

公式(1)～(3)的变量说明如下：

为状态向量，

其中q₁为叶片1平面外一阶模态位移，q₂为叶片2平面外一阶模态位移，q₃为叶片3平面外一阶模态位移，q₄为风轮方位角，q₅为发电机方位角，q₆为塔筒前后一阶模态位移，为q的微分；

为控制增量输入向量，Δυ∈R^O为扰动增量输入向量，y ∈R^P为测量输出向量；

A∈R^N×N为状态矩阵，B∈R^N×N为控制输入增益矩阵，Γ∈R^N×O为扰动输入增益矩阵，C∈R^P×N为状态输出矩阵，为x的微分；

I_rot为风轮转动惯量，I_gen为发电机转动惯量，T_gen为发电机电磁转矩；K_d是传动链扭振的扭转刚度，C_d是传动链扭振阻尼；M_ij为第i个状态变量和第j个状态变量之间的耦合质量系数，C_ij为第i个状态变量和第j个状态变量之间的耦合阻尼系数，K_ij为第i个状态变量和第j个状态变量之间的耦合刚度系数，其中，0＜＜i＜＜j＜＜6；

α＝(δT_aero)/δυ，γ＝(δT_aero)/(δω_r)，ζ＝(δT_aero)/δβ，其中T_aero为风轮的启动转矩，δυ，δω_r，δβ分别是风速、风轮转速和桨距角相对平衡点的增量值；

ζ_t＝(δThrust_aero)/δβ，α_t＝(δThrust_aero)/δυ，Thrust_aero为塔筒的气动推力；

ζi=δTaeroδβi,]]>αi=δTaeroδυi,]]>ζbi=δFaeroiδβi,]]>αbi=δFaeroiδυi,]]>其中δF_aeroi为叶片i平面外一阶挥舞推理增量，β_i，υ_i分别为叶片i的桨距角和叶片i处的风速，i＝1，2，3.ζ_b，α_b分别为ζ_bi和α_bi的平均值；α，γ，ζ，α_t，ζ_t，α_i，ζ_i，α_bi ，ζ_bi均为过程变量；

对公式(3)的状态空间标准形式离散化，得到的线性离散状态空间模型如下：

x(k+1)＝Ax(k)+BΔu(k)+ΓΔυ(k)                   (4)

y(k)＝C_yx(k)                                      (5)

z(k)＝C_zx(k)                                      (6)

其中，x为状态向量，Δu控制输入向量，Δυ为扰动输入向量，y为测量输出向量，z为被控输出变量，C_y∈R^P×N为测量状态输出矩阵，C_z∈R^P×N为被控状态输出矩阵；

步骤3：经风力发电机组全工况预测模型计算出的控制增量，作为下一时刻的桨距角增量或转矩增量，通过变桨驱动器和变流器对叶片和和发电机进行调整控制。

3.根据权利要求2所述的风力发电机组全工况模型预测控制的方法，其特征在于：步骤2所述的通过风力发电机组全工况预测模型计算控制增量，具体按如下步骤进行：

步骤2.1：初始化；

初始化是通过反馈信息测量器检测预测模型的状态变量x(k)，并把它设定为下一步的模型预测初值，从步骤2.2起即转入实时控制；

步骤2.2：采用卡尔曼滤波算法对不可测状态变量叶片一阶挥舞位移进行状态估计，将估计后的状态变量输出；

在每一个采样时刻，都需要预测未来一段时间的被控输出变量值；

设当前时刻为k，信号Δu(k+j|k)表示输入Δu的未来在(k+j)时刻的值，信号x(k+j|k)、z(k+j|k)的含义同Δu(k+j|k)一样；

在每一时刻k，要确定从该时起的起N_c个控制增量，使被控对象在其作用下未来的N个时刻的输出预测值尽可能接近期望，N为预测时域，N_c为控制时域，总是假设N_c＜N，即当j＞N_c时Δu(k+j|k)＝0，即u(k+j|k)＝u(k+N_c-1|k)，j＞N_c，也就是说，当j＞N_c时，一个零阶保持器作用在输入上；

利用公式(4)进行迭代计算，

当j≤N_c时有

x(k+j|k)=Ajx(k)+Aj-1Aj-2...IBΔu(k|k)···Δu(k+j-1|k)]]>

+Aj-1Aj-2...I...0ΓΔu(k|k)···Δu(k+j-1|k)---(7)]]>

当N_c≤j≤N时有

x(k+j|k)=Ajx(x)+Aj-1...Aj-Nc...IBΔu(k|k)···Δu(k+Nc+1|k)]]>

+Aj-1...I...0ΓΔu(k|k)···Δu(k+N-1|k)---(8)]]>

在每一时刻k的预测时域N内总的预测方程为：

定义P为状态变量的系数矩阵，G_z为控制输入增量的系数矩阵，G_υ为控制输入增量的系数矩阵，

则式(9)可写为：

x(k+1|k)···x(k+Nc|k)x(k+Nc+1|k)···x(k+N|k)=px(k)+GzΔu(k|k)···Δu(k+Nc-1|k)+GυΔυ(k|k)···Δυ(k+N-1|k)---(10)]]>

控制输出z的预测方程(6)可以写作：

z(k+j|k)＝C_zx(k+j|k)j＝1，...，N                   (11)

令Z(k)=z(k+1|k)···z(k+N|k),]]>ΔU(k)=Δu(k|k)···Δu(k+Nc-1|k),]]>Δυ(k)=Δυ(k|k)···Δυ(k+N-1|k)]]>则公式(10)可以写作：

Z(k)＝C_zPx(k)+C_zG_zΔU(k)+C_zG_uΔυ(k)                (12)

预测模型中的不可直接测量的状态变量：叶片平面外一阶挥舞位移，需要进行状态估计，状态估计采用卡尔曼滤波器算法，具体如下：

计算当前的估计状态：

x^(k|k)=x^(k|k-1)+M(k)[y(k)-Czx^(k|k-1)]---(13)]]>

其中，表示在时刻k的状态变量的估计，表示在时刻k-1对下一时刻状态变量的估计值，M(k)是反馈信息测量器增益，

M(k)＝P(k|k-1)C_z^T[C_zP(k|k-1)C_z^T+R]^-1                (14)

其中P(k|k-1)是误差方差阵：P(k|k)＝[I-M(k)C_z]P(k|k-1)

计算下一步更新状态和误差方差阵：

x^(k+1|k)=Ax^(k|k)+Bu(k)---(15)]]>

P(k+1|k)＝AP(k|k)A^T+BQB^T              (16)

其中R，Q是高斯白噪声的协方差阵；

使用代替x(k|k)计算预测输出；

步骤2.3：根据风速增量和估计后的状态变量，判断系统所处的运行工况，选择匹配的预测模型，即式(1)或式(2)；

若风速在额定风速和额定风速以下时投入转矩控制预测模型，控制量为转矩，状态变量包括：风轮转速、传动链扭转位移、传动链扭转速度、叶片平面外一阶挥舞位移、叶片平面外一阶挥舞速度、塔筒前后一阶摆动位移和塔筒前后一阶摆动速度，并将风速增量作为扰动输入变量；若风速在额定风速以上时投入变桨控制预测模型，控制量为桨距角，状态变量包括：风轮转速、传动链扭转位移和传动链扭转速度，并将风速增量作为扰动输入变量；

步骤2.4：利用检测到的风力发电机组的状态变量和所选择的预测模型参数，根据式(12)计算控制时域内的控制输出变量；

步骤2.5：将预测值与上一时刻的实际输出进行比较，得出输出误差，如式(18)；

在下一采样时刻首先检测对象的实际输出y(k+1)，并把它与式(11)算出的模型预测输出z(k+1)相比较，构成输出误差

e(k+1)＝y(k+1)-z(k+1|k)              (18)

步骤2.6：将下一时刻的预测模型输出与本时刻的实际输出进行比较，得出输出误差，输出误差与本时刻的预测模型输出值叠加，作为在线校正后的预测值，如式(19)；

采用对e(k+1)加权的方式修正对未来输出的预测：

Z~(k+1)=Z~(k)+he(k+1)---(19)]]>

其中Z~(k+1)=z~(k+1|k+1)···z~(k+N|k+1)]]>为校正后的输出预测向量，有权系数组成的N维向量h=h1···hN]]>为校正向量；

步骤2.7：根据在线校正过的预测值，确定参考轨迹；

采用的参考轨迹是一阶指数变化形式

z_r(k+j)＝a^jz(k)+(1-a^j)z_r           (17)

其中，T_g为参考轨迹的时间常数，下标r表示参考值，z(k)为当前时刻的控制输出反馈值，z_r为设定值，根据预先设定的风速和控制目标的关系曲线得到；

步骤2.8：引入目标函数Jk惩罚校正后的被控输出和参考轨迹z_r(k+j|k)之间的偏差，进行最优化求解，从而得到即时控制增量；

目标函数定义如下：

minJk[x(k),u]=Σj=1N||zr(k+j|k)-z~(k+j|k)||Q(j)2+Σj=0Nc-1||Δu(k+j|k)||R(j)2---(20)]]>

约束条件为：

zmin≤z~(k+j|k)≤zmaxj=0,...,N]]>

Δu_min≤Δu(k+j|k)≤Δu_max        j＝0，…，N_c-1        (21)

其中Q(j)，R(j)是权系数，分别表示对跟踪误差和控制量变化的抑制；

将目标函数(20)写成向量形式：

minJk=||Zr(k)-Z~(k)||Q2+||ΔU(k)||R2---(22)]]>

其中

Zr(k)=zr(k+1|k)···zr(k+N|k),]]>Z~(k)=z~(k+1|k)···z~(k+N|k),]]>ΔU(k)=Δu(k|k)···Δu(k+Nc-1|k)]]>

由权系数构成的对角阵Q和R分别称为误差权矩阵和控制权矩阵；

定义：

E(k)=Zr(k)-Z~(k)=Zr(k)-CzPx(k)-CzGυΔυ(k)---(23)]]>

这个向量可以认为是轨迹误差，也就是目标轨迹和系统自由响应之间的偏差；这样，将目标函数改写为：

minJk=||CzGzΔU(k)-E(k)||Q2+||ΔU(k)||R2]]>

=[ΔUT(k)(CzGz)T-E(k)T]Q[(CzGz)TΔU(k)-E(k)]+ΔUT(k)RΔU(k)]]>

=ΔUT(k)[(CzGz)tQ(CzGz)+R]ΔU(k)-2E(k)TQCzGzΔU(k)+E(k)TQE(k)---(24)]]>

最终可以写成如下的标准形式：

其中H＝2((C_zG_z)^TQ(C_zG_z)+R)，f＝-2E(k)^TQC_zG_z，

约束条件(18)可以写成一个矩阵的形式

I-ICzGz-CzGzΔU(k)≤b1-b2z1-CzPx(k)-CzGυΔυ(k)-z1+CzPx(k)+CzGυΔυ(k)---(26)]]>

其中b₁，b₂的维数为N_c×n_u，分别包括了N_c个Δu_max，Δu_min；同样，z₁，z₂的维数为N，包含了N个z_max，z_min；

这样，问题转化为求解如下的有约束优化问题

minΔU(k)12ΔUT(k)HΔU(k)+fΔU(k)---(27)]]>

并且满足公式(23)，滚动优化要求在每个控制周期仅仅使用优化解得第一列即控制时域中的第一步，这样即时控制增量可写做：

Δuopt(k)=Inu0...0ΔUopt(k)---(28).]]>