本发明属于能源技术领域，具体涉及一种应用于风机的非线性经济模型预测控制方法。
背景技术：
现如今风力发电技术是依靠风机来实现的。风机是一个将风的动能转化为不同形式的电能的能源转换装置。目前，有四种类型的风机：定速定桨距、定速变桨距、变速定桨距、变速变桨距。定桨距是指风机的桨叶与轮毂是刚性连接，风机的桨距角保持不变。变桨距则是柔性连接，叶片的桨距角允许改变。定桨距风机与变桨距风机在叶片设计上是不一样的。定桨距风机的桨叶设计本身有失速特性，当风速超过额定风速时通过桨叶的被动失速调节，使风机转速不随风速增大而一直增大，从而限定了功率。变桨距风机当风速超过额定风速时根据风速的变化而调节桨叶的迎风角度，在一定范围内使风机的转速不随风速的大小而急剧变化。定桨距风机依赖于叶片独特的翼型结构。叶片本身结构较复杂导致成型工艺难度也较大。随着功率增大，叶片加长，所承受的气动推力大，使得叶片的刚度减弱，失速动态特性不易控制，所以很少应用在兆瓦级以上的大型风力发电机组的功率控制上。定速变桨距风机风速范围较小，所以对于兆瓦级以上的大型风力发电机的控制一般选用变速变桨距类型的风机。
变速变桨距风力发电机的控制目标一般为确保输出功率最大化的同时疲劳损失最小。风机是一个多变量、多约束、高度耦合且非仿射型的非线性系统，因此对于风机的控制存在一定挑战。目前已有很多线性或非线性的控制方法应用到风机控制中。如对于有时延的风机模型，采用图形方法设计了pi桨距角控制器；或者根据增益调度l1，在不同的操作点线性化风机模型，采用非凸优化方法设计一个输出反馈控制器。此外还有反步法控制、滑膜控制、h∞控制、鲁棒控制、自适应控制以及预测控制。其中预测控制因其能简单处理多变量多约束的系统，十分适合于兆瓦级以上的风力发电机。在考虑风机实际的操作约束和稳定性时，大多数风机的预测控制设计方案是双层控制结构，上层通过优化经济性能指标获得稳态设定值,传递到下层实现设定值跟踪。这种传统的跟踪式预测控制的目标函数可能无法直接体现经济性，忽略了动态跟踪过程中的经济性能，控制效果欠佳。
技术实现要素：
本发明的目的是提供一种应用于风机的非线性经济模型预测控制方法，解决了现有技术中存在的跟踪式预测控制效果欠佳、较难实现经济性的问题。
本发明所采用的技术方案是，应用于风机的非线性经济模型预测控制方法，具体按照以下步骤实施：
步骤1、基于风机原理进行机理建模；
步骤2、选取发电机转子角速度、风轮转子角速度、轴扭转角、塔架偏移量、塔架偏移速度作为风机模型的状态量，发电机转矩和桨距角作为风机模型的控制量，发电机转子角速度和发电机功率作为风机模型的输出量，以状态空间方程的形式来表示风机非线性模型；
步骤3、利用四阶龙格库塔方法对上述连续风机非线性模型离散化得到离散的风机非线性预测模型；
步骤4、控制发电机转矩和桨距角实现风力发电系统的控制任务要求：低风速时以最大风能捕获模式运行，高风速时发电功率稳定在额定功率不变，选取代表性的经济指标组成经济目标函数，进而体现风力发电系统的经济性能；
步骤5、设定状态量、控制量和输出量的约束条件；
步骤6、结合步骤3离散风机非线性预测模型，步骤4目标函数以及步骤5的各物理量约束条件，利用优化函数进行求解，得到最优控制量，至此完成风机的非线性经济模型预测控制方法的设计。
本发明的特点还在于，
步骤1具体为：
步骤1.1、自然风作用到叶片上，分别产生了空气动力学转矩tr和推力ft，其中tr作用于风轮转子上推动风轮旋转，ft作用于塔架上；tr、ft计算公式如下：
其中，ρ是空气密度，r是叶片半径，π是圆周率，vr是作用在叶片上的有效风速，cp和ct分别是风能利用系数和推力系数，均与叶尖速比λ和桨距角β有关，其中叶尖速比λ的定义为：
步骤1.2、将塔架抽象化为具有弹簧和阻尼器的刚体，建模时简化为一个二阶微分方程：
步骤1.3、用vt表示
步骤1.4、设风机的低速轴、高速轴和齿轮箱之间柔性连接，低速轴、高速轴、齿轮箱组成风机传动系统，将风机传动系统建模成弹簧阻尼质量系统，对柔性传动轴，扭转角定义为风轮转子角位置与发电机转子角位置偏差的绝对值，因此其导数为
步骤1.5、风轮转子由空气动力学转矩tr驱动并由低速轴制动，因此：
发电机转子由高速轴驱动并由发电机转矩tg制动，因此：
以上公式中，jr是风轮转子和低速轴的整体转动惯量，wg是发电机的角速度，jg是高速轴、齿轮箱和发电机的整体转动惯量，tg是发电机转矩，ng是齿轮比，kθ是风机传动系统整体刚度，bθ是风机传动系统整体阻尼系数，风机从自然风中捕获风能，最终由发电机产生的电功率为pe＝ηgwgtg，其中ηg表示发电机效率。
步骤1.1中采用拟合得到风能利用系数，其表达式为：
步骤2具体为：
从所述步骤1中选取发电机转子角速度、风轮转子角速度、轴扭转角、塔架偏移量和塔架偏移速度作为风机模型的状态量x，发电机转矩和桨距角作为风机模型的控制量u，发电机角速度和发电机功率作为风机模型的输出量y，将上述信息写成以下状态空间方程的形式：
y＝gc(x，u，w)，
其中，
fc(x,u,w)为反映状态量数学关系的函数，与状态量、控制量和风速有关；gc(x,u,w)为反映输出量数学关系的函数，与状态量、控制量和风速有关。
步骤3具体为：
步骤3.1、风机的状态表达式为：
步骤3.2、设一个采样周期t内风机的控制量u和自然风速w保持不变，则
步骤3.3、离散后的风机非线性模型表示为x(k+1)＝f(x(k)，u(k),w(k))，风机在k时刻的状态量x(k)是由上一时刻的状态量、控制量以及自然风速决定，k时刻的输出量y(k)也可以表示为y(k)＝g(x(k),u(k),w(k))，由当前时刻的状态量、控制量和自然风速决定，至此，得到离散的风机非线性预测模型。
步骤4中经济指标具体如下：
第一个经济指标le1为功率最大化le1(x,u)＝-pe；
第二个经济指标le2为轴扭转角上的疲劳损失和塔架上的疲劳损失最小化
第三个经济指标le3减小桨叶上的疲劳损害le3(x,u)＝δβ2+α2β2；
第四个经济指标为
对上述四个经济指标线性加权，得到要优化的经济目标函数le(x,u)为：
le(x,u)＝λ1le1+λ2le2+λ3le3+λ4le4；
以上公式中α1表示塔架偏移速度vt相对于轴扭转角θ的重视程度，α2表示桨距角β相对于桨距角变化率δβ的重视程度，同理，λ1,λ2,λ3,λ4表示四个经济指标le1,le2,le3,le4相对于其他指标的重要程度，经济目标中的各个经济指标均是指标值越小表示经济性越好。
步骤5中状态量、控制量和输出量的约束条件具体为：
发电机角速度、转子角速度和轴扭转角满足：
wgmin≤wg≤wgmax,wrmin≤wr≤wrmax，θmin≤θ≤θmax，
发电机转矩和转矩变化率满足：
0≤tg≤tgmax,-δtgmax≤δtg≤δtgmax，
桨距角和桨距角变化率满足：
βmin≤β≤βmax，-δβmax≤δβ≤δβmax，
发电机功率满足：
0≤pe≤perated；
以上公式中凡是带下标min表示此物理量的最小值，带下标max表示此物理量的最大值，带下标rated表示此物理量的额定值，以上最大最小值和额定值是风机的固有属性，与风机容量有关，一旦风机容量确定，具体数值也就确定。
步骤6具体如下：
步骤6.1、经济模型预测控制问题的标准形式具体为：
其中，x(k+i|k)表示在k时刻状态量的第k+i步预测值，u(k+i-1|k)表示在k时刻控制量的k+i步预测值，y(k+i-1|k)表示在k时刻输出量的k+i-1步预测值，w(k+i-1)表示在k+i-1时刻的实际风速值，δu表示控制量的变化率；
步骤6.2、采用内点法对以上控制问题进行优化，得到最优发电机转矩和桨距角实现风机的控制。
本发明的有益效果是，一种应用于风机的非线性经济模型预测控制方法，首先基于风机原理进行机理建模；选取发电机转子角速度、风轮转子角速度、轴扭转角、塔架偏移量、塔架偏移速度作为风机模型的状态量，发电机转矩和桨距角作为风机模型的控制量，发电机转子角速度和发电机功率作为风机模型的输出量，以状态空间方程的形式来表示风机非线性模型；然后离散化风机非线性预测模型；选取代表性的经济指标组成经济目标函数，进而体现风力发电系统的经济性能；设定状态量、控制量和输出量的约束条件；最后利用优化函数进行求解，得到最优控制量。
具体实施方式
下面结合具体实施方式对本发明进行详细说明。
本发明应用于风机的非线性经济模型预测控制方法，具体按照以下步骤实施：
步骤1、基于风机原理进行机理建模；
步骤2、选取发电机转子角速度、风轮转子角速度、轴扭转角、塔架偏移量、塔架偏移速度作为风机模型的状态量，发电机转矩和桨距角作为风机模型的控制量，发电机转子角速度和发电机功率作为风机模型的输出量，以状态空间方程的形式来表示风机非线性模型；
步骤3、利用四阶龙格库塔方法对上述连续风机非线性模型离散化得到离散的风机非线性预测模型；
步骤4、控制发电机转矩和桨距角实现风力发电系统的控制任务要求：低风速时以最大风能捕获模式运行，高风速时发电功率稳定在额定功率不变，选取代表性的经济指标组成经济目标函数，进而体现风力发电系统的经济性能；
步骤5、设定状态量、控制量和输出量的约束条件；
步骤6、结合步骤3离散风机非线性预测模型，步骤4目标函数以及步骤5的各物理量约束条件，利用优化函数进行求解，得到最优控制量，至此完成风机的非线性经济模型预测控制方法的设计。
其中，步骤1具体为：
步骤1.1、自然风作用到叶片上，分别产生了空气动力学转矩tr和推力ft，其中tr作用于风轮转子上推动风轮旋转，ft作用于塔架上；tr、ft计算公式如下：
其中，ρ是空气密度，r是叶片半径，π是圆周率，vr是作用在叶片上的有效风速，cp和ct分别是风能利用系数和推力系数，均与叶尖速比λ和桨距角β有关，其中叶尖速比λ的定义为：
步骤1.2、将塔架抽象化为具有弹簧和阻尼器的刚体，建模时简化为一个二阶微分方程：
步骤1.3、用vt表示
步骤1.4、设风机的低速轴、高速轴和齿轮箱之间柔性连接，低速轴、高速轴、齿轮箱组成风机传动系统，将风机传动系统建模成弹簧阻尼质量系统，对柔性传动轴，扭转角定义为风轮转子角位置与发电机转子角位置偏差的绝对值，因此其导数为
步骤1.5、风轮转子由空气动力学转矩tr驱动并由低速轴制动，因此：
发电机转子由高速轴驱动并由发电机转矩tg制动，因此：
以上公式中，jr是风轮转子和低速轴的整体转动惯量，wg是发电机的角速度，jg是高速轴、齿轮箱和发电机的整体转动惯量，tg是发电机转矩，ng是齿轮比，kθ是风机传动系统整体刚度，bθ是风机传动系统整体阻尼系数，风机从自然风中捕获风能，最终由发电机产生的电功率为pe＝ηgwgtg，其中ηg表示发电机效率。
步骤1.1中采用拟合得到风能利用系数，其表达式为：
步骤2具体为：
从所述步骤1中选取发电机转子角速度、风轮转子角速度、轴扭转角、塔架偏移量和塔架偏移速度作为风机模型的状态量x，发电机转矩和桨距角作为风机模型的控制量u，发电机角速度和发电机功率作为风机模型的输出量y，将上述信息写成以下状态空间方程的形式：
y＝gc(x，u，w)，
其中，
x＝[wgwrθytvt]t,u＝[tgβ]t,y＝[wgpe]t
fc(x,u，w)为反映状态量数学关系的函数，与状态量、控制量和风速有关；gc(x，u,w)为反映输出量数学关系的函数，与状态量、控制量和风速有关。
步骤3具体为：
步骤3.1、风机的状态表达式为：
步骤3.2、设一个采样周期t内风机的控制量u和自然风速w保持不变，则
步骤3.3、离散后的风机非线性模型表示为x(k+1)＝f(x(k)，u(k)，w(k))，风机在k时刻的状态量x(k)是由上一时刻的状态量、控制量以及自然风速决定，k时刻的输出量y(k)也可以表示为y(k)＝g(x(k)，u(k),w(k))，由当前时刻的状态量、控制量和自然风速决定，至此，得到离散的风机非线性预测模型。
步骤4中经济指标具体如下：
第一个经济指标le1为功率最大化le1(x,u)＝-pe；
第二个经济指标le2为轴扭转角上的疲劳损失和塔架上的疲劳损失最小化
第三个经济指标le3减小桨叶上的疲劳损害le3(x,u)＝δβ2+α2β2；
第四个经济指标为
对上述四个经济指标线性加权，得到要优化的经济目标函数le(x,u)为：
le(x,u)＝λ1le1+λ2le2+λ3le3+λ4le4；
以上公式中α1表示塔架偏移速度vt相对于轴扭转角θ的重视程度，α2表示桨距角β相对于桨距角变化率δβ的重视程度，同理，λ1,λ2,λ3,λ4表示四个经济指标le1,le2,le3,le4相对于其他指标的重要程度，经济目标中的各个经济指标均是指标值越小表示经济性越好。
步骤5中状态量、控制量和输出量的约束条件具体为：
发电机角速度、转子角速度和轴扭转角满足：
wgmin≤wg≤wgmax,wrmin≤wr≤wrmax,θmin≤θ≤θmax，
发电机转矩和转矩变化率满足：
0≤tg≤tgmax,-δtgmax≤δtg≤δtgmax，
桨距角和桨距角变化率满足：
βmin≤β≤βmax,-δβmax≤δβ≤δβmax，
发电机功率满足：
0≤pe≤perated；
以上公式中凡是带下标min表示此物理量的最小值，带下标max表示此物理量的最大值，带下标rated表示此物理量的额定值，以上最大最小值和额定值是风机的固有属性，与风机容量有关，一旦风机容量确定，具体数值也就确定。
步骤6具体如下：
步骤6.1、经济模型预测控制问题的标准形式具体为：
其中，x(k+i|k)表示在k时刻状态量的第k+i步预测值，u(k+i-1|k)表示在k时刻控制量的k+i步预测值，y(k+i-1|k)表示在k时刻输出量的k+i-1步预测值，w(k+i-1)表示在k+i-1时刻的实际风速值，δu表示控制量的变化率；
步骤6.2、采用内点法对以上控制问题进行优化，得到最优发电机转矩和桨距角实现风机的控制。
本发明应用于风机的非线性经济模型预测控制方法，采用美国国家可再生能源实验室(nrel)5mw风力发电机模型，经典的风机模型有3个自由度，附加考虑塔架的影响，使风机模型增加了2个自由度。假设自然风速已知，自然风对叶片产生压力，继而产生空气动力学转矩和推力分别作用到风机转子和塔架上，齿轮箱连接风机转子侧和发电机侧，在发电机侧产生转矩带动发电机工作，进而产生电能，通过适当调整接入电网。将风机模型离散化，基于风机实际运行时存在操作约束，同时为了保证经济性，选择合适的目标函数，据此形成了关于风机的非线性经济模型预测控制问题，应用内点法进行问题求解。
由于推力会使塔架有所偏移，造成作用到叶片上的有效风速发生变化，因此上述的风机模型附加了与塔架相关的状态量，整理得到一个包含5个状态量的非线性风机模型，这5个状态量分别为转子角速度、发电机角速度和因两者惯性不同导致的转子角位置与发电机角位置之间的绝对差值，塔架偏移量、塔架偏移速度。利用四阶龙格库塔方法对上述模型离散化得到离散的非线性预测模型。
上述的目标函数的选取是将功率最大化、轴扭转角上的疲劳损失和塔架上的疲劳损失、桨叶上的负载以及平滑控制直接作为优化指标，对这些指标线性加权得到优化目标，所得的目标函数是一个完全体现经济性的目标。
上述的约束选取很全面。除了塔架偏移量和塔架偏移速度外，对于风机模型中的状态量：转子角速度、发电机角速度、轴扭转角，控制量：发电机转矩、桨距角，输出量：发电机功率，都可通过查找nrel5mw风机的定义文档获得其上下限，将该上下限作为控制问题的硬约束。
本发明应用于风机的非线性经济模型预测控制方法，仿真验证如下：
一、阶跃变化的低风速：
采样时间为0.5s，预测时域为50，风速是阶跃变化的，在0-50s时为7.5m/s，在50-100s时为5m/s,在100-150s时为9m/s，目标函数中各个指标前的权重λ1＝-20,λ2＝λ3＝100,λ4＝0.1,α1＝α2＝1。当风速为7.5m/s时，基于在低风速时风机遵循最大风能捕获的原则，可计算得到此时风机模型的状态量，将其作为初值。在matlab中运行改进的非线性经济模型预测控制方法的程序，与传统的预测控制方法比较控制性能更好。
二、特定范围下的风速：
工程定义下的额定风速为11.4m/s，低于额定风速和高于额定风速，风机的工作情况完全不相同。为了验证低于额定风速和高于额定风速风机的工作情况是否符合实际情况，选择了在额定风速附近变化的风速。采样时间、预测时域、模型初值和目标函数中的权重同一，风速以每50s增加0.1m/s的增幅从10.5m/s变化到12m/s，在matlab中运行改进的非线性经济模型预测控制方法的程序，发现工程定义的额定风速符合了经济性。
三、随机变化的风速：
风能软件fast能产生任意变化的风速且fast内包含实际的风机模型，fast子程序可联结到simulink中的s-function模块，因此，fast可与matlabsimulink连接，将nrel5mw风机模型的模型参数写入fast的基本输入文件，在本发明中对风机的控制只涉及到了转矩/功率控制和桨距角控制，且转矩和桨距角控制都来源于外界simulink，因此设置pcmode、vscontrl分别为2和3，将改进的非线性经济模型预测控制方法写入simulink的matlabfunction模块，分别采用两类风速：一是均值为7m/s，从5m/s变化到10m/s的低风速，二是均值为20m/s，从16m/s变化到24m/s的高风速，两类风速都变化600s。运行simulink，验证了改进的非线性经济模型预测控制方法的有效性和实用性。
本发明的方法考虑实时过程经济优化，控制效果更好，具体表现为超调量更小，响应时间更短，变化曲线更平滑。而且，本发明的方法验证了工程上选择的风机额定风速考虑了经济性。