CN102588211A - 一种风力发电机组全工况模型预测控制方法及*** - Google Patents

Abstract

一种风力发电机组全工况模型预测控制方法及***，***包括MPC控制器、反馈信息测量器、风轮、传动链、塔筒、发电单元、变桨驱动器、变流器，风轮、传动链、塔筒和发电单元的状态变量通过反馈信息测量器检测得到，检测结果传递给MPC控制器，通过MPC控制器计算出叶片桨距角和发电机转矩的目标，通过变桨驱动器和变流器对叶片桨距角和风力发电机转矩进行调整。本方法采用变桨控制预测模型和转矩控制预测模型计算控制增量，考虑了传动链扭转位移、传动链扭转速度、叶片平面外一阶挥舞位移、叶片平面外一阶挥舞速度、塔筒前后一阶摆动位移、塔筒前后一阶摆动速度、机组的机械载荷等状态变量，两个预测模型可在不同的工况下进行自动切换，实现全工况运行。

Description

一种风力发电机组全工况模型预测控制方法及***

技术领域

本发明涉及风力发电机组控制技术领域，具体涉及一种风力发电机组全工况模型预测控制方法及***。

背景技术

模型预测控制(MPC)又称基于模型的预测控制(MBPC)，最初由Richalet和Culter等人提出，从1968年模型预测控制的四个基本原则提出发展至今，已经在很多工业控制领域有所应用。模型预测控制包含四项原则，预测模型、参考响应轨迹、滚动优化、反馈校正。

模型预测控制的第一个原则是预测模型。预测控制是一种基于模型的算法，该模型称为预测模型。这是模型预测控制的基础，因为模型的好坏直接影响控制的质量，另外它也是模型预测控制不同于普通PID控制的根本所在。预测控制根据***的未来输入及其历史信息预测未来输出值，预测模型强调其预测功能而不论其形式根据预测模型展示的未来行为，依据某种参考轨迹及优化标准，结合反馈信息，及时调整，以求出未来输入序列作为有自适应机制的预测控制，预测模型不是一成不变的，它直接或间接地被改变预测控制的发展。

模型预测控制的第二个原则是参考响应轨迹。由于控制量的幅值受实际条件的限制，只能在一定的范围内变化，而且控制量的变化速度一般也在一定的范围之内，另外被控对象往往受到不可测干扰的影响，这些因素使得被控变量很难在一拍之内跟踪设定值，一般需要经历一段时间输出才能跟踪设定值，因此根据设定值建立参考响应轨迹，使被控变量沿着参考响应轨迹逐步趋向设定值，这样可防止直接跟踪可能带来的大幅超调而导致的长时间振荡。它是设定值经过柔化的产物。参考响应轨迹概念是Propoi于1963年提出的。参考响应轨迹在模型预测控制中的应用是控制领域的一大进步，也是它优于预估控制的一个重要方面，它使得输出预测的时间长度可以超过时滞长度。

模型预测控制的第三个原则是滚动优化。模型预测控制的优化是一种反馈控制律，其控制量的大小是根据预测输出偏差来调整的。预测输出偏差指输出预测值偏离参考输出值的大小，参考输出值根据参考响应轨迹确定。预测控制的优化不同于传统意义下的优化，它是一种变时段的滚动优化。控制量的预测同输出预测一样不只是下一个控制时刻的预测，在每一采样时刻，优化性能指标涉及从该时刻起到未来的某时段，这不同于预估控制只预测下一个控制时刻的控制量的大小。在下一采样时刻，优化时段同时向后推移，优化是反复在线进行的。模型预测控制在每一个控制时刻只把未来第一个控制时刻的预测控制值应用于控制，其它时刻的预测控制值用于目标函数的性能指标计算。滚动优化是预测控制区别于传统最优控制的根本点。传统的最优控制是用一个性能函数来判断全局最优化，预测控制的优化不是一次离线完成，而是反复进行的，因此适用于动态特性变化和存在不确定因素的复杂工业过程。

模型预测控制的第四个原则是反馈校正。预测控制是一种闭环控制方法，通过滚动优化，确定一系列未来输入后，预测控制仅把当前时刻的控制输入作用于对象，而不是逐一实施。各个输入序列到下一采样时刻，检测实际输出，并根据这一实时信息直接或间接地对原预测模型进行修正(比如，优化时加入预测值与实际值的误差分量或根据实时输出信息在线辨识修改模型)，然后进行新的优化。输入控制序列不像传统的最优控制一次算出，而是反复在线求出提高了对模型失配及环境干扰的鲁棒性。

预测控制***的另一个优越之处就在于通过对具有鲜明物理意义的参数调节(如控制时域、优化时域、参考轨迹时间常数、性能指标中的加权系数)，极易调节成好的鲁棒性这给操作实现带来极大的方便，因此，模型预测控制也是极方便的控制方法学。

模型预测控制已经被广泛地应用于催化裂化、原油蒸馏、延迟焦化、聚乙烯反应器等大型工业过程控制中，但在风力发电机组的控制中很少有人使用。风力发电机组在每个运行阶段的控制量都为了是***趋向于其控制目标，可以理解为快速的过程控制，因此，模型预测控制方法完全可以应用于风力发电机组。

预测控制的主要特征是：以预测模型为基础，采用二次在线滚动优化性能指标和反馈校正等控制策略来消除被控对象建模误差，克服结构、参数与环境等不确定性因素的影响，有效地弥补了现代控制理论对复杂受控对象本身无法避免的不足。针对风电的被控对象建模困难，结构与参数具有不确定，风资源大扰动等特点，更具备适用性。

目前风力发电机组的控制器存在的一个重要问题就是控制作用有一定量的延迟。现有控制器的控制作用大多是通过对前一时刻测量值的评估而做出的，当前的控制量直接基于上一时刻的测量值，因此，如果当前时刻外作用发生波动，控制***在下一时刻才能响应，从而使得控制作用有一定量的滞后。

中国专利公开号CN101560950A，公开日是2009年10月21日，名称为“基于预测控制风力涡轮机的方法和装置”中公开了一种通过预测叶片周围风场来进行风轮桨距控制的方法，该专利以控制叶片的载荷为目标，通过预测时间序列得到风场预测模型，进而利用空气动力学模型计算得到空气产生的载荷值，从而给出用于叶片调整的控制值，将风力涡轮机上的气流影响和相应的叶片调整之间的延迟时间降低。

该专利给出的方法能够弥补以往控制方法有控制延迟的缺陷，相比以往的控制方法，有很大的进步。但该控制方法本质上不是模型预测控制。首先，该控制的预测是对风速的预测，而不是对***未来输出状态的预测；其次，其反馈校正是对风俗预测模型的校正，而不是对控制量误差的校正；然后，该控制方法没有使用滚动优化和参考轨迹。

因此该方法本身存在不足，如没有考虑动态特性变化和不确定因素的影响，同时可能带来大幅超调，导致长时间振荡。

风力发电机组的控制变量主要是桨距角和转矩，在不同的运行阶段，其控制目标不同，采用的控制变量也不同。一台变速型风力发电机组典型的运行阶段和其运行轨迹如图1所示。

阶段1是并网转速之前，控制***主要是监控风速状况是否满足启动条件，满足条件则执行启动过程，打开桨距角增加双馈风电机组转速到并网转速切入电网后转换到阶段12。

阶段12是阶段1和阶段2的过渡阶段，按一个设定的斜率增加转矩给定到阶段2的最优给定曲线上即可。

阶段2是低风速最优功率跟踪阶段，在阶段2控制发电机转矩从风中捕获最大的能量。

阶段23是阶段2和阶段3的过渡阶段，按一个设定的斜率增加转矩给定到阶段3。

阶段3是变桨控制限制功率阶段，控制变桨保证机组输出功率稳定在额定值。这个阶段双馈风电机组承受载荷比较大，对塔筒的晃动，叶片的挥舞以及主传动链的扭振都必须加以控制。

由上分析可以看出，风电机组控制器的控制目标是：1)风电机组在整个运行范围内稳定可靠地在上述5个阶段运行；2)优化机组的运行性能，提高机组的发电效率与质量，减小机组的机械载荷。

《华北电力大学学报》2008年第35期中“基于混合逻辑动态模型的变速风力机全工况运行的控制研究”中介绍了风力机全工况运行下的混杂***模型，该文结合风力机在不同工况点的非线性模型，引入混合逻辑动态模型的建模方法，建立风力机全工况运行下的混杂***模型，并在此模型基础上，引入广义模型预测控制进行最优控制，使在风速变化的情况下，风力机可以良好的跟踪设定的功率曲线。

其不足之处是其控制目标只考虑了使风电机组在整个运行范围内稳定可靠地运行，没有考虑风电机组的运行性能优化，如提高机组的发电质量，减小机组的机械载荷等。并且，其控制模型中没有考虑风速扰动的影响，无法保证控制策略在风速扰动时实现机组的稳定运行和快速跟随。

现有技术中，预测风速却没有考虑动态特性变化和不确定因素的影响；没有参考轨迹可能带来大幅超调，导致长时间振荡；未考虑风电机组的运行性能优化，如提高机组的发电质量，减小机组的机械载荷等；没有考虑风速扰动的影响，无法保证控制策略在风速扰动时实现机组的稳定运行和快速跟随。

发明内容

针对现有技术存在的不足，本发明提供一种风力发电机组全工况模型预测控制方法及***。

一种风力发电机组全工况模型预测控制***，包括MPC控制器、反馈信息测量器、风轮、传动链、塔筒、发电单元、变桨驱动器、变流器；

风轮、传动链、塔筒和发电单元为被控对象，他们的状态变量通过反馈信息测量器检测得到；检测结果传递给MPC控制器，通过MPC控制器计算出叶片桨距角和发电机转矩的调整值；变桨驱动器和变流器为执行装置，通过它们对风轮和发电单元进行调整。

发电单元的输出端接至电网，发电单元风轮配备的变桨驱动器，调整变桨角度，实现风能最大利用率；变桨驱动器输入端与MPC控制器输出端口相连，变流器输入端连接到MPC控制器的输出端口，变流器输出端连接到发电单元。

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

所述MPC控制器，MPC控制器选用PLC，PLC包括机舱子站和塔底子站，机舱子站和塔底子站之间通过FAST总线通信。内部嵌入状态估计模块、预测模块、在线校正模块、参考轨迹模块和滚动优化模块，预测模块中，定义两个预测模型：转矩控制预测模型和变桨控制预测模型。两个预测模型可在不同的工况下进行自动切换，在额定风速和额定风速以下时投入控制量为转矩的预测模型，额定风速以上时投入控制量为桨距角的预测模型；参考轨迹模块根据***的工况选择对应的参考轨迹；所述预测模型为计入风速扰动变量的离散状态空间模型；

本***的工作过程如下：风轮、传动链、塔筒和发电单元的状态变量通过监测装置检测得到，检测结果传递给MPC控制器，通过MPC控制器计算出叶片桨距角和发电机转矩的调整值，通过变桨驱动器和变流器对叶片和风力发电机进行调整。

采用上述控制***进行风力发电机组全工况模型预测控制的方法，具体按如下步骤进行：

步骤1：***的离线准备；

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

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

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

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

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

$\left[\begin{array}{c}{x^{\·}}_1\\ {x^{\·}}_2\\ {x^{\·}}_3\\ {x^{\·}}_4\\ {x^{\·}}_5\\ {x^{\·}}_6\\ {x^{\·}}_7\end{array}\right]={\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}$

$\·\left[\begin{array}{ccccccc}\mathrm{\γ}& {-K}_d& {-C}_d& 0& {-3C}_{40}& 0& -C_{46}\\ 0& 0& 1& 0& 0& 0& 0\\ \frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -(\frac{K_d}{I_{\mathrm{gen}}}+\frac{K_d}{I_{\mathrm{rot}}})& -(\frac{C_d}{I_{\mathrm{gen}}}+\frac{C_d}{I_{\mathrm{rot}}})& 0& -\frac{{3C}_{40}}{I_{\mathrm{rot}}}& 0& -\frac{C_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 0& 1& 0& 0\\ {-C}_{04}& 0& 0& {-K}_{00}& {-C}_{00}& 0& 0\\ 0& 0& 0& 0& 0& 0& 1\\ {-C}_{64}& 0& 0& {-K}_{60}& {-C}_{60}& {-K}_{60}& {-C}_{66}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\end{array}\right]$

$+{\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}\left[\begin{array}{c}{3\mathrm{\ζ}}_1\\ 0\\ \frac{{3\mathrm{\ζ}}_1}{I_{\mathrm{rot}}}\\ 0\\ {\mathrm{\ζ}}_b\\ 0\\ {\mathrm{\ζ}}_t\end{array}\right]\mathrm{\δ\β}$

$+{\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}\left[\begin{array}{c}{3\mathrm{\α}}_1\\ 0\\ \frac{{3\mathrm{\α}}_1}{I_{\mathrm{rot}}}\\ 0\\ {\mathrm{\α}}_b\\ 0\\ {\mathrm{\α}}_t\end{array}\right]\mathrm{\δ\ν}---\left(1\right)$

$y=C\·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\end{array}\right]$

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

$\left[\begin{array}{c}{x^{\·}}_1\\ {x^{\·}}_2\\ {x^{\·}}_3\end{array}\right]\left[\begin{array}{ccc}\frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -\frac{K_d}{I_{\mathrm{ror}}}& -\frac{C_d}{I_{\mathrm{rot}}}\\ 0& 0& 1\\ \frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -(\frac{K_d}{I_{\mathrm{gen}}}+\frac{K_d}{I_{\mathrm{rot}}})& -(\frac{C_d}{I_{\mathrm{gen}}}+\frac{C_d}{I_{\mathrm{rot}}})\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ \frac{1}{I_{\mathrm{gen}}}\end{array}\right]{\mathrm{\δT}}_{\mathrm{gen}}+\left[\begin{array}{c}\frac{\mathrm{\α}}{I_{\mathrm{rot}}}\\ 0\\ \frac{\mathrm{\α}}{I_{\mathrm{rot}}}\end{array}\right]\mathrm{\δ\ν}---\left(2\right)$

$y=C\·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$

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

$\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{B\Δu}+\mathrm{\Γ\Δ\υ}---\left(3\right)$

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×M为控制输入增益矩阵，Γ∈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为塔筒的气动推力。

${\mathrm{\ζ}}_i=\frac{{\mathrm{\δT}}_{\mathrm{aero}}}{{\mathrm{\δ\β}}_i},$ ${\mathrm{\α}}_i=\frac{{\mathrm{\δT}}_{\mathrm{aero}}}{{\mathrm{\δ\υ}}_i},$ ${\mathrm{\ζ}}_{\mathrm{bi}}=\frac{{\mathrm{\δF}}_{\mathrm{aeroi}}}{{\mathrm{\δ\β}}_i},$ ${\mathrm{\α}}_{\mathrm{bi}}=\frac{{\mathrm{\δF}}_{\mathrm{aeroi}}}{{\mathrm{\δ\υ}}_i},$ 其中δF_aeroi为叶片i平面外一阶挥舞推理增量，β_υυ_i分别为叶片i的桨距角和叶片i处的风速，i＝1，2，3.ζ_b，α_b分别为ζ_bi和α_bi的平均值。α，γ，ζ，α_t，ζ_c，α_t，ζ_t，α_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为被控状态输出矩阵。

通过风力发电机组全工况预测模型计算控制增量，具体按如下步骤进行：

步骤2.1：初始化；

模型预测的在线计算由初始化与实时控制组成。

初始化是通过反馈信息测量器检测预测模型的状态变量x(k)，并把它设定为下一步的模型预测初值，从步骤2.2起即转入实时控制，在每一个采样时刻的在线计算流程如图3所示。

步骤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)=A^{j}x\left(k\right)+\begin{array}{}\end{array}\left[\begin{array}{cccc}{A}^{j-1}& A^{j-2}& ...& I\end{array}\right]B\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+j-1|k)\end{array}\right]$

$+\left[\begin{array}{cccccc}{A}^{j-1}& A^{j-2}& ...& I& ...& 0\end{array}\right]\mathrm{\Γ}\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+j-1|k)\end{array}\right]---\left(7\right)$

当N_c≤j≤N时有

$x(k+j|k)=A^{j}x\left(x\right)+\left[\begin{array}{ccccc}{A}^{j-1}& ...& A^{j-N_c}& ...& I\end{array}\right]B\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c+1|k)\end{array}\right]$

$+\left[\begin{array}{ccccc}{A}^{j-1}& ...& I& ...& 0\end{array}\right]\mathrm{\Γ}\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N-1|k)\end{array}\right]---\left(8\right)$

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

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

则式(9)可写为：

$\left[\begin{array}{c}x(k+1|k)\\ \·\\ \·\\ \·\\ x(k+N_c|k)\\ x(k+N_c+1|k)\\ \·\\ \·\\ \·\\ x(k+N|k)\end{array}\right]=\mathrm{px}\left(k\right)+G_z\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right]+G_{\mathrm{\υ}}\left[\begin{array}{c}\mathrm{\Δ\υ}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δ\υ}(k+N-1|k)\end{array}\right]---\left(10\right)$

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

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

令 $Z\left(k\right)=\left[\begin{array}{c}z(k+1|k)\\ \·\\ \·\\ \·\\ z(k+N|k)\end{array}\right],$ $\mathrm{\ΔU}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right],$ $\mathrm{\Δ\υ}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ\υ}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δ\υ}(k+N-1|k)\end{array}\right]$

则公式(10)可以写作：

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

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

计算当前的估计状态：

$\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+M\left(k\right)[y\left(k\right)-C_z\hat{x}\left(k\right|k-1)]---\left(13\right)$

其中，

表示在时刻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)

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

$\hat{x}(k+1|k)=A\hat{x}\left(k\right|k)+\mathrm{Bu}\left(k\right)---\left(15\right)$

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)加权的方式修正对未来输出的预测：

$\tilde{Z}(k+1)=\tilde{Z}\left(k\right)+\mathrm{he}(k+1)---\left(19\right)$

其中 $\tilde{Z}(k+1)=\left[\begin{array}{c}\tilde{z}(k+1|k+1)\\ \·\\ \·\\ \·\\ \tilde{z}(k+N|k+1)\end{array}\right]$ 为校正后的输出预测向量，有权系数组成的N维向量 $h=\left[\begin{array}{c}{h}_1\\ \·\\ \·\\ \·\\ h_N\end{array}\right]$ 为校正向量。

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

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

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

其中，

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

步骤2.8：引入目标函数Jk惩罚校正后的被控输出

和参考轨迹z_r(k+j|k)之间的偏差，进行最优化求解，从而得到即时控制增量。

目标函数定义如下：

$\min J_k[x\left(k\right),u]=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|z_r(k+j|k)-\tilde{z}(k+j|k)\left|\right|}_{Q\left(j\right)}^2+\underset{j=0}{\overset{N_c-1}{\mathrm{\Σ}}}{\left|\right|\mathrm{\Δu}(k+j|k)\left|\right|}_{R\left(j\right)}^2---\left(20\right)$

约束条件为：

$\begin{array}{cc}{z}_{\min }\≤\tilde{z}(k+j|k)\≤z_{\max }& j=0,...,N\end{array}$

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

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

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

$\min J_k={\left|\right|Z_r\left(k\right)-\tilde{Z}\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_R^2---\left(22\right)$

其中

$Z_r\left(k\right)=\left[\begin{array}{c}{z}_r(k+1|k)\\ \·\\ \·\\ \·\\ z_r(k+N|k)\end{array}\right],$ $\tilde{Z}\left(k\right)=\left[\begin{array}{c}\tilde{z}(k+1|k)\\ \·\\ \·\\ \·\\ \tilde{z}(k+N|k)\end{array}\right],$ $\mathrm{\ΔU}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right]$

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

定义：

$E\left(k\right)=Z_r\left(k\right)-\tilde{Z}\left(k\right)=Z_r\left(k\right)-C_z\mathrm{Px}\left(k\right)-C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)---\left(23\right)$

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

$\min J_k={\left|\right|C_z{G}_z\mathrm{\ΔU}\left(k\right)-E\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_R^2$

$=[{\mathrm{\ΔU}}^T\left(k\right){\left(C_z{G}_z\right)}^T-E{\left(k\right)}^T]Q[{\left(C_z{G}_z\right)}^T\mathrm{\ΔU}\left(k\right)-E\left(k\right)]+{\mathrm{\ΔU}}^T\left(k\right)\mathrm{R\ΔU}\left(k\right)$

$={\mathrm{\ΔU}}^T\left(k\right)[{\left(C_z{G}_z\right)}^{t}Q\left(C_z{G}_z\right)+R]\mathrm{\ΔU}\left(k\right)-{2E\left(k\right)}^T{\mathrm{QC}}_z{G}_z\mathrm{\ΔU}\left(k\right)+{E\left(k\right)}^T\mathrm{QE}\left(k\right)---\left(24\right)$

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

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

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

$\left[\begin{array}{c}I\\ -I\\ C_z{G}_z\\ -C_z{G}_z\end{array}\right]\mathrm{\ΔU}\left(k\right)\≤\left[\begin{array}{c}{b}_1\\ -b_2\\ z_1-C_z\mathrm{Px}\left(k\right)-C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)\\ -z_1+C_z\mathrm{Px}\left(k\right)+C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)\end{array}\right]---\left(26\right)$

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

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

${\min }_{\mathrm{\ΔU}\left(k\right)}\frac{1}{2}{\mathrm{\ΔU}}^T\left(k\right)\mathrm{H\ΔU}\left(k\right)+\mathrm{f\ΔU}\left(k\right)---\left(27\right)$

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

${\mathrm{\Δu}}_{\mathrm{opt}}\left(k\right)=\left[\begin{array}{cccc}{I}_{n_u}& 0& ...& 0\end{array}\right]{\mathrm{\ΔU}}_{\mathrm{opt}}\left(k\right)---\left(28\right)$

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

有益效果：

1.本发明控制***采用的模型预测控制方法，由于具有滚动优化模块和参考轨迹模块，***本身考虑了动态特性变化和不确定因素的影响；同时避免了***大幅超调，导致长时间振荡的现象；

2.本发明的预测模型包含传动链扭转位移、传动链扭转速度、叶片平面外一阶挥舞位移、叶片平面外一阶挥舞速度、塔筒前后一阶摆动位移、塔筒前后一阶摆动速度，减小机组的机械载荷等状态变量，同时在优化模型中增加了各状态变量的约束条件，可以有效控制机组的机械载荷；

3.本发明的预测模型包含风速扰动变量，将每一采样时刻的风速扰动都计入了控制策略中，可以保证在风速扰动时实现机组的稳定运行和快速跟随；

4.本发明方法有两个预测模型，变桨控制预测模型和转矩控制预测模型，两个预测模型可在不同的工况下进行自动切换，实现了全工况运行。

附图说明

图1是一台变速型风力发电机组典型的运行阶段和其运行轨迹图；

图2是本发明实施例风力发电机组全工况模型预测控制原理图；

图3是本发明实施例模型MPC控制器原理框图；

图4为本发明实施例在线计算流程图；

图5为本发明实施例模型预测控制器PLC配置方案构成图；

图6为本发明实施例比较阶跃风输入下***的响应曲线图；

图7为本发明实施例比较阶跃风输入下叶片挥舞位移曲线图；

图8为本发明实施例比较阶跃风输入下塔筒顶部前后摆动位移曲线图；

图9为本发明实施例试验湍流风输入曲线图；

图10为本发明实施例比较湍流风输入***响应曲线图；

图11为本发明实施例比较湍流风输入下叶片挥舞位移曲线图；

图12为本发明实施例比较湍流风输入塔筒顶部前后摆动位移曲线图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

本实施例采用以3MW双馈风力发电机组的全工况模型预测控制***，包括MPC控制器1、反馈信息测量器6、风轮3、传动链4、塔筒5、发电单元8、变桨驱动器2、变流器7、控制原理如图2所示，发电单元的输出端接至电网，发电单元风轮配备的变桨驱动器，调整叶片桨距角，实现风能最大利用率；变桨驱动器输入端与MPC控制器输出端口相连，变流器输入端连接到MPC控制器的输出端口，变流器输出端连接到发电单元。

所述MPC控制器，内部嵌入状态估计模块、预测模块、在线校正模块、参考轨迹模块和滚动优化模块，原理如图3所示。预测模块中，定义两个预测模型：转矩控制预测模型和变桨控制预测模型。两个预测模型可在不同的工况下进行自动切换，在额定风速和额定风速以下时投入控制量为转矩的预测模型，额定风速以上时投入控制量为桨距角的预测模型；参考轨迹模块根据***的工况选择对应的参考轨迹；所述预测模型为计入风速扰动变量的离散状态空间模型。

MPC控制器选用PLC，PLC包括机舱子站和塔底子站，分别置于机舱控制柜和塔底控制柜中，机舱子站和塔底子站之间通过FAST总线通信，如图5所示。

机舱子站配置如下：FS211/N为FAST总线子模块；DIO280为数字输入/输出模块；DI232为数字输入模块；PTAI216为温度记录模块；AIO288为模拟输入输出模块；CNT204为计数模块；RS204为串口通信模块；CM202为CAN总线模块；EM203为以太网主站模块；LM201为占位模块；WT310为显示屏。

图中塔底子站配置如下：SWI205为工业以太网交换机模块；NT250为电源模块；MPC240为处理器；PCC201为存储卡模块；FM211为FAST总线主模块；DIO248为数字输入/输出模块；PTAI216为温度记录模块；LM201为占位模块。

AIO288模块、RS204模块是MPC控制器与反馈信息测量器6的通信接口。反馈信息测量器6包括：安装于风力发电机主轴的转速传感器(普通接近开关)、安装于滑环尾部的绝对位置编码器(型号为HENGSTLER-AC58)、安装于风力发电机尾部的发电机编码器(型号为LeineLinde-RHI503)、安装于机舱底盘的振动分析仪(型号为PCH1024)、安装于发电机定子侧的电压测量传感器(型号为SIEMENS-3RT)、安装于发电机定子侧的电流测量传感器(型号为SIEMENS-ESM500-9920)、装设与机舱外部支架上的风速传感器(型号为FT702)和安装于变距减速机旁的叶片位置编码器(型号为Baumer-SSI)。

CM202模块是MPC控制器与变桨控制器2、变流器7的通信接口，它们之间通过CAN总线通信。

FS211模块与FM211模块是机舱柜和塔底柜通信的接口模块，通过FAST总线通信。

SWI205模块是MPC控制器与塔底编程/维护设备、风场SCADA、Mobile SMS的通信接口。

PCC201为存储卡模块，存储的数据有状态矩阵A、控制输入增益矩阵B、扰动输入增益矩阵Γ、测量状态输出矩阵C_y、被控状态输出矩阵C_z、被控输入和被控输出的约束条件、被控输出设定值z_r随风速变化曲线、采样周期T、预测时域N、控制时域N_c、参考轨迹时间常数T_g、校正向量h、输出误差权重Q、控制量变化幅度权重R、历史风速、过程变量，以便处理器模块MPC240随时调用。

***的工作过程如下：风轮、传动链、塔筒和发电单元的状态变量通过监测装置检测得到，检测结果传递给MPC控制器，通过MPC控制器计算出叶片桨距角和发电机转矩的调整值，通过变桨驱动器和变流器对叶片和风力发电机进行调整。

采用上述控制***进行风力发电机组全工况模型预测控制的方法，具体按如下步骤进行：

步骤1：***的离线准备；

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

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

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

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

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

$\left[\begin{array}{c}{x^{\·}}_1\\ {x^{\·}}_2\\ {x^{\·}}_3\\ {x^{\·}}_4\\ {x^{\·}}_5\\ {x^{\·}}_6\\ {x^{\·}}_7\end{array}\right]={\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}$

$\·\left[\begin{array}{ccccccc}\mathrm{\γ}& {-K}_d& {-C}_d& 0& {-3C}_{40}& 0& -C_{46}\\ 0& 0& 1& 0& 0& 0& 0\\ \frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -(\frac{K_d}{I_{\mathrm{gen}}}+\frac{K_d}{I_{\mathrm{rot}}})& -(\frac{C_d}{I_{\mathrm{gen}}}+\frac{C_d}{I_{\mathrm{rot}}})& 0& -\frac{{3C}_{40}}{I_{\mathrm{rot}}}& 0& -\frac{C_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 0& 1& 0& 0\\ {-C}_{04}& 0& 0& {-K}_{00}& {-C}_{00}& 0& 0\\ 0& 0& 0& 0& 0& 0& 1\\ {-C}_{64}& 0& 0& {-K}_{60}& {-C}_{60}& {-K}_{60}& {-C}_{66}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\end{array}\right]$

$+{\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}\left[\begin{array}{c}{3\mathrm{\ζ}}_1\\ 0\\ \frac{{3\mathrm{\ζ}}_1}{I_{\mathrm{rot}}}\\ 0\\ {\mathrm{\ζ}}_b\\ 0\\ {\mathrm{\ζ}}_t\end{array}\right]\mathrm{\δ\β}$

$+{\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}\left[\begin{array}{c}{3\mathrm{\α}}_1\\ 0\\ \frac{{3\mathrm{\α}}_1}{I_{\mathrm{rot}}}\\ 0\\ {\mathrm{\α}}_b\\ 0\\ {\mathrm{\α}}_t\end{array}\right]\mathrm{\δ\ν}---\left(1\right)$

$y=C\·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\end{array}\right]$

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

$\left[\begin{array}{c}{x^{\·}}_1\\ {x^{\·}}_2\\ {x^{\·}}_3\end{array}\right]\left[\begin{array}{ccc}\frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -\frac{K_d}{I_{\mathrm{ror}}}& -\frac{C_d}{I_{\mathrm{rot}}}\\ 0& 0& 1\\ \frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -(\frac{K_d}{I_{\mathrm{gen}}}+\frac{K_d}{I_{\mathrm{rot}}})& -(\frac{C_d}{I_{\mathrm{gen}}}+\frac{C_d}{I_{\mathrm{rot}}})\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ \frac{1}{I_{\mathrm{gen}}}\end{array}\right]{\mathrm{\δT}}_{\mathrm{gen}}+\left[\begin{array}{c}\frac{\mathrm{\α}}{I_{\mathrm{rot}}}\\ 0\\ \frac{\mathrm{\α}}{I_{\mathrm{rot}}}\end{array}\right]\mathrm{\δ\ν}---\left(2\right)$

$y=C\·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$

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

$\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{B\Δu}+\mathrm{\Γ\Δ\υ}---\left(3\right)$

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×M为控制输入增益矩阵，Γ∈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为塔筒的气动推力。 ${\mathrm{\ζ}}_i=\frac{{\mathrm{\δT}}_{\mathrm{aero}}}{{\mathrm{\δ\β}}_i},$ ${\mathrm{\α}}_i=\frac{{\mathrm{\δT}}_{\mathrm{aero}}}{{\mathrm{\δ\υ}}_i},$ ${\mathrm{\ζ}}_{\mathrm{bi}}=\frac{{\mathrm{\δF}}_{\mathrm{aeroi}}}{{\mathrm{\δ\β}}_i},$ ${\mathrm{\α}}_{\mathrm{bi}}=\frac{{\mathrm{\δF}}_{\mathrm{aeroi}}}{{\mathrm{\δ\υ}}_i},$ 其中δF_aeroi为叶片i平面外一阶挥舞推理增量，β_i，υ_u分别为叶片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_x∈R^P×N为被控状态输出矩阵。

通过风力发电机组全工况预测模型计算控制增量，具体按如下步骤进行：

步骤2.1：初始化；

模型预测的在线计算由初始化与实时控制组成，流程如图4所示。

初始化是通过反馈信息测量器检测预测模型的状态变量x(k)，并把它设定为下一步的模型预测初值，从步骤2.2起即转入实时控制，在每一个采样时刻的在线计算流程如图3所示。

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

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

假设当前时刻为k，信号Δu(k+j|k)表示输入Δu的未来在(k+j)时刻的值，信号x(k+j|k)、z(k+j|k)的含义同Δ(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)=A^{j}x\left(k\right)+\begin{array}{}\end{array}\left[\begin{array}{cccc}{A}^{j-1}& A^{j-2}& ...& I\end{array}\right]B\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+j-1|k)\end{array}\right]$

$+\left[\begin{array}{cccccc}{A}^{j-1}& A^{j-2}& ...& I& ...& 0\end{array}\right]\mathrm{\Γ}\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+j-1|k)\end{array}\right]---\left(7\right)$

当N_c≤j≤N时有

$x(k+j|k)=A^{j}x\left(x\right)+\left[\begin{array}{ccccc}{A}^{j-1}& ...& A^{j-N_c}& ...& I\end{array}\right]B\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c+1|k)\end{array}\right]$

$+\left[\begin{array}{ccccc}{A}^{j-1}& ...& I& ...& 0\end{array}\right]\mathrm{\Γ}\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N-1|k)\end{array}\right]---\left(8\right)$

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

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

则式(9)可写为：

$\left[\begin{array}{c}x(k+1|k)\\ \·\\ \·\\ \·\\ x(k+N_c|k)\\ x(k+N_c+1|k)\\ \·\\ \·\\ \·\\ x(k+N|k)\end{array}\right]=\mathrm{px}\left(k\right)+G_z\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right]+G_{\mathrm{\υ}}\left[\begin{array}{c}\mathrm{\Δ\υ}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δ\υ}(k+N-1|k)\end{array}\right]---\left(10\right)$

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

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

令 $Z\left(k\right)=\left[\begin{array}{c}z(k+1|k)\\ \·\\ \·\\ \·\\ z(k+N|k)\end{array}\right],$ $\mathrm{\ΔU}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right],$ $\mathrm{\Δ\υ}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ\υ}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δ\υ}(k+N-1|k)\end{array}\right]$

则公式(10)可以写作：

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

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

计算当前的估计状态：

$\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+M\left(k\right)[y\left(k\right)-C_z\hat{x}\left(k\right|k-1)]---\left(13\right)$

其中，

表示在时刻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)

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

$\hat{x}(k+1|k)=A\hat{x}\left(k\right|k)+\mathrm{Bu}\left(k\right)---\left(15\right)$

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)加权的方式修正对未来输出的预测：

$\tilde{Z}(k+1)=\tilde{Z}\left(k\right)+\mathrm{he}(k+1)---\left(19\right)$

其中 $\tilde{Z}(k+1)=\left[\begin{array}{c}\tilde{z}(k+1|k+1)\\ \·\\ \·\\ \·\\ \tilde{z}(k+N|k+1)\end{array}\right]$ 为校正后的输出预测向量，有权系数组成的N维向量 $h=\left[\begin{array}{c}{h}_1\\ \·\\ \·\\ \·\\ h_N\end{array}\right]$ 为校正向量。

步骤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)之间的偏差，进行最优化求解，从而得到即时控制增量。

目标函数定义如下：

$\min J_k[x\left(k\right),u]=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|z_r(k+j|k)-\tilde{z}(k+j|k)\left|\right|}_{Q\left(j\right)}^2+\underset{j=0}{\overset{N_c-1}{\mathrm{\Σ}}}{\left|\right|\mathrm{\Δu}(k+j|k)\left|\right|}_{R\left(j\right)}^2---\left(20\right)$

约束条件为：

$\begin{array}{cc}{z}_{\min }\≤\tilde{z}(k+j|k)\≤z_{\max }& j=0,...,N\end{array}$

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

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

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

$\min J_k={\left|\right|Z_r\left(k\right)-\tilde{Z}\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_R^2---\left(22\right)$

其中

$Z_r\left(k\right)=\left[\begin{array}{c}{z}_r(k+1|k)\\ \·\\ \·\\ \·\\ z_r(k+N|k)\end{array}\right],$ $\tilde{Z}\left(k\right)=\left[\begin{array}{c}\tilde{z}(k+1|k)\\ \·\\ \·\\ \·\\ \tilde{z}(k+N|k)\end{array}\right],$ $\mathrm{\ΔU}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right]$

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

定义：

$E\left(k\right)=Z_r\left(k\right)-\tilde{Z}\left(k\right)=Z_r\left(k\right)-C_z\mathrm{Px}\left(k\right)-C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)---\left(23\right)$

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

$\min J_k={\left|\right|C_z{G}_z\mathrm{\ΔU}\left(k\right)-E\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_R^2$

$=[{\mathrm{\ΔU}}^T\left(k\right){\left(C_z{G}_z\right)}^T-E{\left(k\right)}^T]Q[{\left(C_z{G}_z\right)}^T\mathrm{\ΔU}\left(k\right)-E\left(k\right)]+{\mathrm{\ΔU}}^T\left(k\right)\mathrm{R\ΔU}\left(k\right)$

$={\mathrm{\ΔU}}^T\left(k\right)[{\left(C_z{G}_z\right)}^{t}Q\left(C_z{G}_z\right)+R]\mathrm{\ΔU}\left(k\right)-{2E\left(k\right)}^T{\mathrm{QC}}_z{G}_z\mathrm{\ΔU}\left(k\right)+{E\left(k\right)}^T\mathrm{QE}\left(k\right)---\left(24\right)$

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

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

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

$\left[\begin{array}{c}I\\ -I\\ C_z{G}_z\\ -C_z{G}_z\end{array}\right]\mathrm{\ΔU}\left(k\right)\≤\left[\begin{array}{c}{b}_1\\ -b_2\\ z_1-C_z\mathrm{Px}\left(k\right)-C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)\\ -z_1+C_z\mathrm{Px}\left(k\right)+C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)\end{array}\right]---\left(26\right)$

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

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

${\min }_{\mathrm{\ΔU}\left(k\right)}\frac{1}{2}{\mathrm{\ΔU}}^T\left(k\right)\mathrm{H\ΔU}\left(k\right)+\mathrm{f\ΔU}\left(k\right)---\left(27\right)$

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

${\mathrm{\Δu}}_{\mathrm{opt}}\left(k\right)=\left[\begin{array}{cccc}{I}_{n_u}& 0& ...& 0\end{array}\right]{\mathrm{\ΔU}}_{\mathrm{opt}}\left(k\right)---\left(28\right)$

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

本***变桨控制预测模型和转矩控制预测模型切换风速为14m/s，即风速在14m/s以下时使用转矩控制预测模型，在14m/s以上时使用变桨控制预测模型。

其变桨控制线性离散预测模型是将变桨控制预测模型在阶段3的稳态工作点按控制周期0.02s离散化后获得，其阶段3的稳态工作点条件为：

风速υ＝18(m/s)；风轮转速ω_r＝14.25(rpm)；叶片桨距角β＝16.66(deg)。

其变桨控制状态线性离散预测模型参数如表1所示。

表13 MW风力发电机组变桨控制7状态线性离散预测模型参数

其转矩控制线性离散预测模型是将转矩控制预测模型在阶段2的稳态工作点按控制周期0.02s离散化后获得，其阶段2的稳态工作点条件为：

风速υ＝13(m/s)；风轮转速ω_r＝14.25(rpm)；叶片桨距角β＝6.941(deg)。

其转矩控制状态线性离散预测模型参数如表2所示。

表23 MW风力发电机组转矩控制3状态线性离散预测模型参数

试验阶跃风速输入下风电机组响应，风速输入为17m/s在250s处跳跃为18m/s。

控制目标(被控输出设定值z_r)为：稳定风轮转速为14.25rpm、增加传动链模态阻尼、增加叶片挥舞模态阻尼和塔筒前后模态阻尼，传动链扭振、叶片一阶挥舞和塔筒前后一阶摆动偏移为零。

模型的输入输出约束为：

-9(deg/s)≤Δβ≤9(deg/m)

8.9(rpm)≤ω_r≤16(rpm)

采样周期T为0.02秒，参考轨迹时间常数T_g为0.1，预测时域N、控制时域N_c选择为10和1，校正向量 $h=\left[\begin{array}{c}1\\ 0.5\\ 0.5\\ 0.5\\ 0\\ 0\\ 0\end{array}\right],$ 输出误差权重 $Q=\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 0& 0.5& 0& 0& 0& 0& 0\\ 0& 0& 0.5& 0& 0& 0& 0\\ 0& 0& 0& 0.5& 0& 0& 0\\ 0& 0& 0& 0& 0.1& 0& 0\\ 0& 0& 0& 0& 0& 0.1& 0\\ 0& 0& 0& 0& 0& 0& 0.1\end{array}\right],$ 控制量变化幅度权重R为300。

在MATLAB中进行仿真，变增益PI控制和模型预测控制采用相同的输入，图6为阶跃风输入下风轮转速响应，图7为叶片平面外一阶挥舞位移，图8为塔筒顶部前后一阶摆动位移，可以看出模型预测控制在稳定传动链扭振、减小叶片挥舞载荷和减小塔筒前后摆动载荷方面效果明显。

对湍流风输入测试，结果如下：图9是测试湍流风，图10是湍流风输入下风轮转速响应，图11是湍流风输入下的叶片一阶挥舞位移，图12是湍流风输入下塔筒顶部前后一阶摆动位移，可以看出湍流风速输入下模型预测控制在载荷控制方面效果也是很明显的。

Claims (3)

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：反馈信息测量器实时检测风力发电机组的状态变量，通过风力发电机组全工况预测模型计算控制增量；

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

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

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

$\left[\begin{array}{c}{x^{\·}}_1\\ {x^{\·}}_2\\ {x^{\·}}_3\\ {x^{\·}}_4\\ {x^{\·}}_5\\ {x^{\·}}_6\\ {x^{\·}}_7\end{array}\right]={\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}$

$\·\left[\begin{array}{ccccccc}\mathrm{\γ}& {-K}_d& {-C}_d& 0& {-3C}_{40}& 0& -C_{46}\\ 0& 0& 1& 0& 0& 0& 0\\ \frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -(\frac{K_d}{I_{\mathrm{gen}}}+\frac{K_d}{I_{\mathrm{rot}}})& -(\frac{C_d}{I_{\mathrm{gen}}}+\frac{C_d}{I_{\mathrm{rot}}})& 0& -\frac{{3C}_{40}}{I_{\mathrm{rot}}}& 0& -\frac{C_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 0& 1& 0& 0\\ {-C}_{04}& 0& 0& {-K}_{00}& {-C}_{00}& 0& 0\\ 0& 0& 0& 0& 0& 0& 1\\ {-C}_{64}& 0& 0& {-K}_{60}& {-C}_{60}& {-K}_{60}& {-C}_{66}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\end{array}\right]$

$+{\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}\left[\begin{array}{c}{3\mathrm{\ζ}}_1\\ 0\\ \frac{{3\mathrm{\ζ}}_1}{I_{\mathrm{rot}}}\\ 0\\ {\mathrm{\ζ}}_b\\ 0\\ {\mathrm{\ζ}}_t\end{array}\right]\mathrm{\δ\β}$

$+{\left[\begin{array}{ccccccc}{I}_{\mathrm{rot}}& 0& 0& 0& {3M}_{40}& 0& M_{46}\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& \frac{{3M}_{40}}{I_{\mathrm{rot}}}& 0& \frac{M_{46}}{I_{\mathrm{rot}}}\\ 0& 0& 0& 1& 0& 0& 0\\ M_{04}& 0& 0& 0& M_{00}& 0& M_{06}\\ 0& 0& 0& 0& 0& 1& 0\\ M_{64}& 0& 0& 0& {3M}_{60}& 0& M_{66}\end{array}\right]}^{-1}\left[\begin{array}{c}{3\mathrm{\α}}_1\\ 0\\ \frac{{3\mathrm{\α}}_1}{I_{\mathrm{rot}}}\\ 0\\ {\mathrm{\α}}_b\\ 0\\ {\mathrm{\α}}_t\end{array}\right]\mathrm{\δ\ν}---\left(1\right)$

$y=C\·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\\ x_7\end{array}\right]$

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

$\left[\begin{array}{c}{x^{\·}}_1\\ {x^{\·}}_2\\ {x^{\·}}_3\end{array}\right]\left[\begin{array}{ccc}\frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -\frac{K_d}{I_{\mathrm{ror}}}& -\frac{C_d}{I_{\mathrm{rot}}}\\ 0& 0& 1\\ \frac{\mathrm{\γ}}{I_{\mathrm{rot}}}& -(\frac{K_d}{I_{\mathrm{gen}}}+\frac{K_d}{I_{\mathrm{rot}}})& -(\frac{C_d}{I_{\mathrm{gen}}}+\frac{C_d}{I_{\mathrm{rot}}})\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ \frac{1}{I_{\mathrm{gen}}}\end{array}\right]{\mathrm{\δT}}_{\mathrm{gen}}+\left[\begin{array}{c}\frac{\mathrm{\α}}{I_{\mathrm{rot}}}\\ 0\\ \frac{\mathrm{\α}}{I_{\mathrm{rot}}}\end{array}\right]\mathrm{\δ\ν}---\left(2\right)$

$y=C\·\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$

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

$\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{B\Δu}+\mathrm{\Γ\Δ\υ}---\left(3\right)$

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为塔筒的气动推力；

${\mathrm{\ζ}}_i=\frac{{\mathrm{\δT}}_{\mathrm{aero}}}{{\mathrm{\δ\β}}_i},$ ${\mathrm{\α}}_i=\frac{{\mathrm{\δT}}_{\mathrm{aero}}}{{\mathrm{\δ\υ}}_i},$ ${\mathrm{\ζ}}_{\mathrm{bi}}=\frac{{\mathrm{\δF}}_{\mathrm{aeroi}}}{{\mathrm{\δ\β}}_i},$ ${\mathrm{\α}}_{\mathrm{bi}}=\frac{{\mathrm{\δF}}_{\mathrm{aeroi}}}{{\mathrm{\δ\υ}}_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)=A^{j}x\left(k\right)+\begin{array}{}\end{array}\left[\begin{array}{cccc}{A}^{j-1}& A^{j-2}& ...& I\end{array}\right]B\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+j-1|k)\end{array}\right]$

$+\left[\begin{array}{cccccc}{A}^{j-1}& A^{j-2}& ...& I& ...& 0\end{array}\right]\mathrm{\Γ}\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+j-1|k)\end{array}\right]---\left(7\right)$

当N_c≤j≤N时有

$x(k+j|k)=A^{j}x\left(x\right)+\left[\begin{array}{ccccc}{A}^{j-1}& ...& A^{j-N_c}& ...& I\end{array}\right]B\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c+1|k)\end{array}\right]$

$+\left[\begin{array}{ccccc}{A}^{j-1}& ...& I& ...& 0\end{array}\right]\mathrm{\Γ}\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N-1|k)\end{array}\right]---\left(8\right)$

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

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

则式(9)可写为：

$\left[\begin{array}{c}x(k+1|k)\\ \·\\ \·\\ \·\\ x(k+N_c|k)\\ x(k+N_c+1|k)\\ \·\\ \·\\ \·\\ x(k+N|k)\end{array}\right]=\mathrm{px}\left(k\right)+G_z\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right]+G_{\mathrm{\υ}}\left[\begin{array}{c}\mathrm{\Δ\υ}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δ\υ}(k+N-1|k)\end{array}\right]---\left(10\right)$

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

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

令 $Z\left(k\right)=\left[\begin{array}{c}z(k+1|k)\\ \·\\ \·\\ \·\\ z(k+N|k)\end{array}\right],$ $\mathrm{\ΔU}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right],$ $\mathrm{\Δ\υ}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δ\υ}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δ\υ}(k+N-1|k)\end{array}\right]$ 则公式(10)可以写作：

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

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

计算当前的估计状态：

$\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+M\left(k\right)[y\left(k\right)-C_z\hat{x}\left(k\right|k-1)]---\left(13\right)$

其中，

表示在时刻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)

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

$\hat{x}(k+1|k)=A\hat{x}\left(k\right|k)+\mathrm{Bu}\left(k\right)---\left(15\right)$

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)加权的方式修正对未来输出的预测：

$\tilde{Z}(k+1)=\tilde{Z}\left(k\right)+\mathrm{he}(k+1)---\left(19\right)$

其中 $\tilde{Z}(k+1)=\left[\begin{array}{c}\tilde{z}(k+1|k+1)\\ \·\\ \·\\ \·\\ \tilde{z}(k+N|k+1)\end{array}\right]$ 为校正后的输出预测向量，有权系数组成的N维向量 $h=\left[\begin{array}{c}{h}_1\\ \·\\ \·\\ \·\\ h_N\end{array}\right]$ 为校正向量；

步骤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)之间的偏差，进行最优化求解，从而得到即时控制增量；

目标函数定义如下：

$\min J_k[x\left(k\right),u]=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left|\right|z_r(k+j|k)-\tilde{z}(k+j|k)\left|\right|}_{Q\left(j\right)}^2+\underset{j=0}{\overset{N_c-1}{\mathrm{\Σ}}}{\left|\right|\mathrm{\Δu}(k+j|k)\left|\right|}_{R\left(j\right)}^2---\left(20\right)$

约束条件为：

$\begin{array}{cc}{z}_{\min }\≤\tilde{z}(k+j|k)\≤z_{\max }& j=0,...,N\end{array}$

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

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

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

$\min J_k={\left|\right|Z_r\left(k\right)-\tilde{Z}\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_R^2---\left(22\right)$

其中

$Z_r\left(k\right)=\left[\begin{array}{c}{z}_r(k+1|k)\\ \·\\ \·\\ \·\\ z_r(k+N|k)\end{array}\right],$ $\tilde{Z}\left(k\right)=\left[\begin{array}{c}\tilde{z}(k+1|k)\\ \·\\ \·\\ \·\\ \tilde{z}(k+N|k)\end{array}\right],$ $\mathrm{\ΔU}\left(k\right)=\left[\begin{array}{c}\mathrm{\Δu}\left(k\right|k)\\ \·\\ \·\\ \·\\ \mathrm{\Δu}(k+N_c-1|k)\end{array}\right]$

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

定义：

$E\left(k\right)=Z_r\left(k\right)-\tilde{Z}\left(k\right)=Z_r\left(k\right)-C_z\mathrm{Px}\left(k\right)-C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)---\left(23\right)$

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

$\min J_k={\left|\right|C_z{G}_z\mathrm{\ΔU}\left(k\right)-E\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_R^2$

$=[{\mathrm{\ΔU}}^T\left(k\right){\left(C_z{G}_z\right)}^T-E{\left(k\right)}^T]Q[{\left(C_z{G}_z\right)}^T\mathrm{\ΔU}\left(k\right)-E\left(k\right)]+{\mathrm{\ΔU}}^T\left(k\right)\mathrm{R\ΔU}\left(k\right)$

$={\mathrm{\ΔU}}^T\left(k\right)[{\left(C_z{G}_z\right)}^{t}Q\left(C_z{G}_z\right)+R]\mathrm{\ΔU}\left(k\right)-{2E\left(k\right)}^T{\mathrm{QC}}_z{G}_z\mathrm{\ΔU}\left(k\right)+{E\left(k\right)}^T\mathrm{QE}\left(k\right)---\left(24\right)$

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

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

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

$\left[\begin{array}{c}I\\ -I\\ C_z{G}_z\\ -C_z{G}_z\end{array}\right]\mathrm{\ΔU}\left(k\right)\≤\left[\begin{array}{c}{b}_1\\ -b_2\\ z_1-C_z\mathrm{Px}\left(k\right)-C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)\\ -z_1+C_z\mathrm{Px}\left(k\right)+C_z{G}_{\mathrm{\υ}}\mathrm{\Δ\υ}\left(k\right)\end{array}\right]---\left(26\right)$

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

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

${\min }_{\mathrm{\ΔU}\left(k\right)}\frac{1}{2}{\mathrm{\ΔU}}^T\left(k\right)\mathrm{H\ΔU}\left(k\right)+\mathrm{f\ΔU}\left(k\right)---\left(27\right)$

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

${\mathrm{\Δu}}_{\mathrm{opt}}\left(k\right)=\left[\begin{array}{cccc}{I}_{n_u}& 0& ...& 0\end{array}\right]{\mathrm{\ΔU}}_{\mathrm{opt}}\left(k\right)---\left(28\right).$

Images