CN116796387B - 基于全信息自适应合作博弈理论的风力机翼型优化方法 - Google Patents

Abstract

本发明公开了基于全信息自适应合作博弈理论的风力机翼型优化方法，涉及风力机叶片翼型优化设计技术领域，该方法包括：步骤1：基于parsec参数法在原始翼型上进行拟合，得到初始设计变量值，修改翼型的几何外形；步骤2：计算翼型的极惯性矩；步骤3：计算翼型年发电量；步骤4：建立翼型约束条件：对翼型厚度t进行约束，并保证优化后的翼型年发电量比原始翼型发电量高；基于初始设计变量值计算设计变量的上下限；步骤5：建立翼型的目标函数：基于步骤2，设置优化目标函数；步骤6：基于全信息自适应合作博弈方法，弈对设计变量和目标函数进行划分、收益函数计算、非劣解的筛选，得到pareto最优解集。本发明优化后的翼型具有良好的发电性能和颤振性能。

Description

基于全信息自适应合作博弈理论的风力机翼型优化方法

技术领域

本发明涉及风力机叶片翼型优化设计技术领域，特别是涉及一种基于全信息自适应合作博弈理论的风力机翼型优化方法。

背景技术

为了提高发电功率，降低发电成本，现代风力发电机叶片的尺寸显著增大。伴随风力机叶片长度的增加，叶片的刚度降低，柔性特征突出。在空气动力、惯性力和弹性力的相互作用下，风力机叶片模态之间会发生振动耦合，从而引发颤振现象。颤振会造成叶片疲劳，产生裂纹，甚至断裂，因此，为了使发电机正常运行，在设计风力机叶片时需要重视其颤振性能的设计。

捕获风能的关键因素与风力机叶片的气动外形有关，叶片截面年发电量与叶片的发电性能有关，通过对风机叶片截面翼型的优化设计，可以提高叶片截面翼型的颤振临界速度和发电性能，降低发电成本和颤振现象的发生。风能是新型可再生能源，应用广泛，伴随现代风力发电技术在不断升级，优化风力机叶片发电性能和颤振性能成为设计叶片的重点。

现有的多目标翼型优化算法，可分为三类：数学求解法；通过分组、分层和分类等方法将多目标转化为单目标的方法；群智能优化方法。其中，加权方法主要是构造加权函数，其中各目标的地位通过权重来体现，这类方法对问题的pareto前沿太敏感，存在很大局限性，难以合理对目标函数之间的影响和冲突进行智能调节。通过分组、分层和分类的方法存在一些缺点，它主要是通过决策者来决定次序和权重，存在一些主观因素，得到的只是可行域内的一个解集，并非最优解集。群智能算法也是应用比较广泛的算法，常见的有遗传算法、粒子群算法、模拟退火算法等，它在很大程度上可以改善多目标问题的求解结果，且具有良好的鲁棒性，但是具有容易在计算过程中陷入局部最优和收敛速度过慢等问题。

发明内容

本发明的目的是提供一种基于全信息自适应合作博弈理论的风力机翼型优化方法，能够实现对翼型上下表面的拟合，将复杂问题简单化，使得优化后的翼型具有更好的发电性能和抗颤振性能。

为实现上述目的，本发明提供了如下方案：

一种基于全信息自适应合作博弈理论的风力机翼型优化方法，包括：

步骤1：基于parsec参数法在原始翼型上进行拟合，得到初始设计变量值x₀，修改翼型的几何外形；

步骤2：计算极惯性矩；

在极惯性矩计算过程中，首先建立极惯性矩与几何外形的关系，然后计算翼型的面积和质心坐标；平移坐标系，使坐标原点和质心重合，根据平行轴原理，求出翼型相对于质心轴系的惯性矩，得到翼型的极惯性矩；

步骤3：计算翼型年发电量；

步骤4：建立翼型约束条件：对翼型厚度t进行约束，并保证优化后的翼型年发电量比原始翼型发电量高，基于初始设计变量值计算设计变量的上下限，即：

厚度约束条件 t₁≤t≤t₂

翼型年发电量约束条件 AEP_opt>AEP_ori

式中，t₁和t₂为翼型厚度的下限和上限，AEP_opt为优化后的截面年发电量，AEP_ori为原截面年发电量；

步骤5：建立翼型的优化目标函数：基于步骤2和步骤3，设置优化目标函数为：

步骤6：基于全信息自适应合作博弈方法进行设计变量和目标函数的优化；即运用全信息自适应合作博弈对设计变量和目标函数进行划分、收益函数计算、非劣解的筛选，最后得到一个pareto最优解集。

进一步的，所述步骤1基于parsec参数法在原始翼型上进行拟合，得到初始设计变量值x₀，修改翼型的几何外形，具体包括：

S101、采取parsec参数法进行原始翼型上下表面拟合，得到翼型上、下表面拟合的形状函数：

上表面：

下表面：

上表面对应系数的取值由下面方程组所确定：

同样的，下表面对应系数的取值由下面方程组所确定：

式中，x_te为后缘位置；x_up为翼型上表面顶部位置；x_lo为翼型下表面底部位置；r_le,up为翼型上表面前缘半径；r_le,lo为翼型下表面前缘半径；z_te为翼型后缘处纵坐标；Δz_te为翼型后缘厚度；z_up为翼型上表面上峰峰值；z_lo为翼型下表面下峰峰值；z_xx,up为翼型上表面曲率；z_xx,lo为翼型下表面曲率；α_te为翼型后缘方向，β_te为翼型后缘楔角；这13个参数对应了优化问题的设计变量，由于x_te恒等于1，因此剩余12个参数组成的设计变量如下所示：

x＝(r_le,up,r_le,lo,x_up,z_up,z_xx,up,x_lo,z_lo,z_xx,lo,Δz_te,z_te,α_te,β_te)。

步骤102、使用最小二乘拟合方法得到原始翼型表面的12个初始设计变量值。

进一步的，所述步骤S102使用最小二乘法拟合得到原始翼型表面的12个设计变量值，得到修改后的翼型的几何外形，具体包括：

原始翼型表面的12个初始设计变量值的表达式为：

其中，z_u(i)代表翼型上翼面的第i个点的纵坐标，z_ou(i)代表原始翼型上翼面的第i个点的纵坐标，z_l(i)代表翼型下翼面的第i个点的纵坐标，z_ol(i)代表原始翼型下翼面的第i个点的纵坐标。

进一步的，所述步骤2中计算翼型的面积和质心坐标点，具体包括：

求翼型的面积：

其中，S为翼型面积；c为翼型弦长，取值为标准值1；x代表横坐标，y代表纵坐标；代表翼型上翼面曲线函数，为翼型上翼面曲线函数。

翼型的质心坐标点：

其中，ρ(x,y)代表翼型不同坐标点处的密度，假设翼型材质是均匀，即ρ(x,y)当做常数，可进一步化简，即：

进一步的，所述步骤2中得出翼型的极惯性矩，具体包括：

若将翼型前缘位置平移至坐标点处，在翼型坐标平面内任意一点的轴惯性矩表示为：

若将平面坐标系的原点平移到质心处，两坐标轴的方向保持不变，则对应的轴惯性矩为：

求得相对质心的惯性矩为：

I＝I_x1+I_y1。

进一步的，所述步骤3计算翼型截面年发电量，具体包括：

S301：设置初始轴向诱导因子a＝0，切向诱导因子b＝0。

S302：计算入流角φ和局部攻角α,计算公式如下：

α＝φ-θ

其中，λ_r为叶片展向r处的局部叶尖速比，B代表风力机的叶片数量，hub_r为风力机轮毂的半径，blade_r为叶片总长度，λ为叶片的叶尖速比，θ为叶片截面处翼型相对于叶尖翼型的扭角；

S303：用Xfoil软件计算该翼型在攻角α下的升力系数cl和阻力系数cd。

S304：计算出法向力系数C_n和切向力系数C_t，计算公式如下：

C_n＝cl cosφ+cd sinφ

C_t＝cl sinφ-cd cosφ

S305：计算此时的轴向诱导因子a，切向诱导因子b，计算公式如下：

其中，a_c＝0.2，σ′表示风轮实度，计算公式如下：

F表示损失系数，计算公式如下：

F＝F_tipF_hub

Prandtl叶尖修正因子

Prandtl叶根修正因子

S306：如果轴向诱导因子a和切向诱导因子b变化小于阈值δ＝1×10^-6，那么结束迭代，得到最终的轴向诱导因子a和切向诱导因子b，否则继续步骤302。

S307：计算出截面功率P，最后结合风速的威尔分布公式得出年发电量AEP：

进一步的，所述步骤6基于全信息自适应合作博弈方法进行设计变量和目标函数优化，具体包括：

S601：设置权重组数，假设共有Q组权重系数，q为当前的组数，权重系数的分配如下：

式中，w_ij代表权重系数；

S602：设置迭代次数，假设最大次数为M，k为当前迭代次数，随机生成初始策略S⁽⁰⁾；

S603：根据自适应策略集，对初始策略S⁽⁰⁾的设计变量划分给策略子集具体划分方法如下：

(1)对m个目标函数优化设计，得到m组最优解：

对应的策略集表示为：

(2)对任意设计变量x_j，对目标f_i的影响度为θ_ji表示为：

其中，Δx_j为步长；

影响度归一化处理：

基于“空间距离”和“力矩”大小对策略空间进行划分，根据设计变量到目标函数之间的“距离空间”衡量远近，在策略子集划分中还要依据力矩的大小来划分，计算公式如下：

力矩的计算公式：

式中，d(i,j)表示设计变量到目标函数的空间距离；Mo(j)表示设计变量对于目标的力矩；

λ为力矩的阈值，定义如下:

(3)将每个目标函数即博弈方的所有设计变量的距离d(j,i)从最小到最大排序；具有相同距离的设计变量按影响度D(j,i)排序；影响度越大，排名越靠前；排名靠前的设计变量被划分给该博弈方，直到设计变量的力矩Mo(j)的累积值超过阈值λ；初步划分后，如果将相同的设计变量分配给多个博弈方，则该设计变量将以更大的影响度D(j,i)重新分配给博弈方；

S604：基于k，判断是否是第一次迭代，如果是第一次迭代，随机初始化为各博弈方策略子集；如果不是第一次迭代，进入步骤605其中，表示第i个博弈方的第N个策略子集；

S605：基于通过模拟二进制交叉和多项式变异得到N个新策略

S606：计算2N个策略对应的最终成本函数并将前N个较小的最终成本函数对应的策略赋值给成本函数的计算方法如下：

式中，表示第i个博弈方在采用策略S_i时各个博弈方的相对成本，即目标函数f_i规格化的无量纲值；和分别表示第i个博弈方的初始策略和对应的补集；w_ij代表权系数，权系数取值大小代表各博弈方之间合作的程度，w_ij取值越大，则说明合作度越高；“*”表示计算两个矩阵的哈达玛积，其元素定义为两个矩阵对应元素的乘积；

S607：得到每个博弈方策略子集，即

S608：将得到的各博弈方策略方子集放在集合B中，即B＝B∪{S⁽¹⁾,S⁽²⁾,...,S^(N)}；更新迭代次数；

S609：判断是否对所有博弈方都进行了博弈，如果没有，则进入S604，否则进入S610；

S610：对迭代当前次数下的所有博弈方的策略子集进行整理，即更新下一次迭代次数，进入S3，如果是最后一次迭代，则进入S611；

S611：更新下一组权重系数的取值，执行S2；如果是最后一组权重系数，则进入步骤612；

S612：把集合B中的非劣解放入集合A中，即A中为全信息自适应合作博弈的最终解集。

进一步的，所述步骤603，根据自适应策略集，对初始策略S⁽⁰⁾的设计变量划分给策略子集中，具体划分方法如下：

(1)对m个目标函数优化设计，得到m组最优解：

对应的策略集表示为：

(2)对任意设计变量x_j，对目标f_i的影响度为θ_ji表示为：

其中，Δx_j为步长；

影响度归一化处理：

基于“空间距离”和“力矩”大小对策略空间进行划分，根据设计变量到目标函数之间的“距离空间”衡量远近，在策略子集划分中依据力矩的大小来划分，计算公式如下：

力矩的计算公式：

式中，d(j，i)表示设计变量到目标函数的空间距离；Mo(j)表示设计变量对于目标的力矩；

λ为力矩的阈值，定义如下:

(3)将每个目标函数即博弈方的所有设计变量的距离d(j,i)从最小到最大排序；具有相同距离的设计变量按影响度D(j,i)排序；距离越大，排名越靠前；排名靠前的设计变量被划分给该博弈方，直到设计变量的力矩Mo(j)的累积值超过阈值λ；初步划分后，如果将相同的设计变量分配给多个博弈方，则该设计变量以更大的影响度D(j,i)重新分配给博弈方。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明提供的基于全信息自适应合作博弈理论的风力机翼型优化方法，前期使用parsec参数法在原始翼型上进行拟合，修改翼型的几何外形，运用全信息自适应合作博弈方法对12个设计变量和2个目标函数进行策略子集和策略空间的划分，通过运算，可以得到逼近pareto边界的优化解集，提高其抗颤振稳定性能并增强其发电性能，使得叶片截面发电量更高。

1.本申请在采取parsec参数法进行翼型上下表面的拟合时，使用12个表面特征函数，不需要基线形状，可以生成较宽的翼型范围，对翼型进行厚度约束时，可以使用简单的定界约束或线性约束来近似表示，更容易预测单个参数对翼型气动性能的影响。

2.本申请采用全信息自适应合作博弈方法，对于目标函数和设计变量之间的对应关系，博弈方法可以通过对博弈策略空间的划分来揭示，而且可以根据对手的策略调整自身策略，将高维、复杂的问题简单化，为多目标优化问题提供一个新的思路。

3.本申请建立了翼型极惯性矩与翼型几何外形的关系，通过全信息自适应合作博弈方法不断优化，能够找到pareto最优解集。

4.相较于传统的优化方法，比如加权法，遗传算法、粒子群优化算法、多目标进化算法等，这类方法是最常见的也是使用比较多的算法，但是其缺点也很明显，为了获得一组pareto最优解，往往计算周期长，甚至难以求解。本申请采用的全信息自适应合作博弈方法，可以清晰的看出目标函数和设计变量之间的关系，将目标函数转化为博弈方，设计变量当做博弈方的策略子集，通过一定的方法把策略子集划分到每一博弈方所对应的策略空间里，还可以通过对方策略来调整自身的策略，这降低了问题求解的难度，缩减了求解时间。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例的原始翼型图；

图2为本发明实施例的叶片截面诱导速度的原理示意图；

图3为本发明实施例的叶片截面上的局部载荷的原理示意图；

图4为本发明实施例的基于自适应全信息合作博弈算法优化设计的pareto最优解集(未加约束)；

图5a-5c为本发明实施例中优化后新翼型与原始翼型的比较图；

图6为本发明实施例中符合约束要求的pareto图；

图7为本发明所述方法的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种基于全信息自适应合作博弈理论的风力机翼型优化方法，能够实现对翼型上下表面的拟合，将复杂问题简单化，使得优化后的翼型具有更好的发电性能和抗颤振性能。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

由经典抗颤振理论的研究，得出翼型的极惯性矩对其抗颤振性能具有突出的影响，其值可通过高斯积分计算得出。叶片的发电性能与叶片截面年发电量有关，其值可通过叶素动量理论(BEM)求得。将叶片截面最大年发电量和极惯性矩作为优化目标，运用全信息自适应博弈算法，优化得出的翼型比原有翼型具有更好的发电性能和抗颤振性能。

实施例1

如图1所示，本实施例中，选择5MV风力机中的DU99-W-350为原始翼型。

本发明实施例提供的基于全信息自适应合作博弈理论的风力机翼型优化方法，包括：

步骤1：基于parsec参数法在原始翼型上进行拟合，得到初始设计变量值x₀，修改翼型的几何外形；

翼型几何外形直接影响极惯性矩和叶片截面年发电量，通过对原始拟合进行参数化拟合可以修改翼型的形状。这样做的优点是：更新后的翼型可以继承原有翼型的大部分优点，也可以节省时间。

步骤2：计算极惯性矩；

在对极惯性矩计算过程中，首先建立极惯性矩与几何外形的关系，然后计算出翼型的面积和质心坐标；平移坐标系，使坐标原点和质心重合，根据平行轴原理，可求出翼型相对于质心轴系的惯性矩，最终得到翼型的极惯性矩。

步骤3：计算翼型年发电量。

在计算翼型年发电量，需要结合叶素动量理论公式，入流角φ时，使用经典叶素动量方程在计算的过程中会发现在大部分的情况下并不收敛，故而引用松弛因子策略，可以使算例均能收敛，更新的入流角公式为φ_i+1＝(1-w)φ_i+wf(φ_i)。

步骤4：建立翼型约束条件：对翼型厚度t进行约束，并保证优化后的翼型年发电量比原始翼型发电量高，基于步骤1的初始设计变量值计算设计变量的上下限，即：

厚度约束条件t₁≤t≤t₂

翼型年发电量约束条件 AEP_opt>AEP_ori

设计变量约束 0.85x₀≤x≤1.15x₀

式中，t₁和t₂为翼型厚度的下限和上限，AEP_opt为优化后的截面年发电量，AEP_ori为原截面年发电量；

步骤5：建立翼型的目标函数：基于步骤2和步骤3得到目标函数：

f₁(x)＝1/I

f₂(x)＝10⁶/AEP

其中，I为翼型的极惯性矩，AEP为年发电量。

步骤6：基于全信息自适应合作博弈方法，对设计变量和目标函数进行划分、收益函数计算、非劣解的筛选，最后得到一个pareto最优解集，如图4所示。

本实施例中，所述步骤1，基于parsec参数法在原始翼型上进行拟合，得到设计变量，修改翼型的几何外形，具体包括：

S101：翼型上下表面拟合的形状函数：

上表面：

下表面：

上表面对应系数的取值由下面方程组所确定：

同样的，下表面对应系数的取值由下面方程组所确定：

式中，x_te为后缘位置；x_up为翼型上表面顶部位置；x_lo为翼型下表面底部位置；r_le,up为翼型上表面前缘半径；r_le,lo为翼型下表面前缘半径；z_te为翼型后缘处纵坐标；Δz_te为翼型后缘厚度；z_up为翼型上表面上峰峰值；z_lo为翼型下表面下峰峰值；z_xx,up为翼型上表面曲率；z_xx,lo为翼型下表面曲率；α_te为翼型后缘方向，β_te为翼型后缘楔角；这13个参数对应了优化问题的设计变量，由于x_te恒等于1，因此剩余12个参数组成的设计变量如下所示：

x＝(r_le,up,r_le,lo,x_up,z_up,z_xx,up,x_lo,z_lo,z_xx,lo,Δz_te,z_te,α_te,β_te)。

S102、使用最小二乘拟合方法得到原始翼型表面的12个初始设计变量值。

进一步的，所述S102具体包括：

将拟合问题转化为误差的平方和最小值问题，求解该问题即可得到原始翼型表面的12个初始设计变量值，表达式如下所示：

计算得到初始设计变量：

x₀＝(0.0965242795453330,0.0990891355391492,0.307013671746253,

0.168646859306725,-1.46228122574921,0.293885008920728,

-0.181492235141888,2.18992410687365,0.00558424803211468,

0.000208439579293796,14.2288001044516,4.24220613257765)

本实施例中，所述步骤2具体包括：

S201、计算出翼型的面积和质心坐标点；

由于极惯性矩I取决于翼型几何外形，而I的求解不能被直接体现出，需要建立几何外形与极惯性矩之间的关系。首先需要计算出翼型的面积和质心坐标点，进而得到相对质心的惯性矩。

S202、得出翼型极惯性矩。

进一步的实施例中，所述S201具体包括：

求翼型的面积：

其中，S为翼型的面积；c为翼型弦长，取值为标准值1；x代表横坐标，y代表纵坐标；代表翼型上翼面曲线函数，为翼型下翼面曲线函数。

翼型质心坐标：

ρ(x,y)代表翼型不同坐标点处的密度，假设翼型材质是均匀，即ρ(x,y)当做常数，可进一步化简，即：

进一步的实施例中，所述S201具体包括：

若将翼型前缘位置平移至坐标点处，在翼型坐标平面内任意一点的轴惯性矩表示为：

若将平面坐标系的原点平移到质心处，两坐标轴的方向保持不变，则对应的轴惯性矩为：

根据上两个公式计算可求得相对质心的惯性矩为：

I＝I_x1+I_y1。

通过对风力机叶片二元翼段进线性行参数敏感性分析，发现临界颤振速度与叶片的扭转刚度、重心位置、弹性轴位置及当地空气密度密切相关。研究结果表明提高叶片扭转刚度、将叶片重心位置前移、弹性轴位置向前缘移动或者将风机置于低密度区域均可以提高抗颤振性能，但是后三者很难实现。所以要提高翼型的抗颤振性能需要提高叶片翼型的扭转刚度。翼型的扭转刚度与三个因素有关：

为扭转刚度，G为剪切模量，与材料的选择有关，选取剪切模量越高的材料其扭转刚度越好，因此抗颤振性能越高；R为翼型所在截面距离叶根的半径，由公式可知，半径和扭转刚度成反比，即R越大，扭转刚度越小，即抗颤振性能减小；I为翼型极惯性矩，与翼型外形密切相关，翼型抗颤振性能随I增大而增加。

本实施例中，所述步骤3计算翼型截面年发电量，具体包括：

S301：设置初始轴向诱导因子a＝0，切向诱导因子b＝0。

S302：计算入流角φ和局部攻角α，计算公式如下：

α＝φ-θ

其中，λ_r为叶片展向r处的局部叶尖速比，B代表风力机的叶片数量，hub_r为风力机轮毂的半径，blade_r为叶片总长度，λ为叶片的叶尖速比，θ为叶片截面处翼型相对于叶尖翼型的扭角；

S303：用Xfoil软件计算该翼型在攻角α下的升力系数cl和阻力系数cd。

S304：计算出法向力系数C_n和切向力系数C_t，计算公式如下：

C_n＝cl cosφ+cd sinφ

C_t＝cl sinφ-cd cosφ

S305：计算此时的轴向诱导因子a，切向诱导因子b，计算公式如下：

其中，a_c＝0.2，σ′表示风轮实度，计算公式如下：

F表示损失系数，计算公式如下：

F＝F_tipF_hub

Prandtl叶尖修正因子

Prandtl叶根修正因子

S306：如果轴向诱导因子a和切向诱导因子b变化小于阈值δ＝1×10^-6，那么结束迭代，得到最终的轴向诱导因子a和切向诱导因子b，否则继续步骤302。

S307：计算出截面功率P，最后结合风速的威尔分布公式得出年发电量AEP：

用叶素动量理论进行分析时，其一运行叶片截面诱导速度和叶素上的局部载荷分析图分别如图2和图3所示。

本实施例中，步骤4中，翼型沿着风力机叶片展向分布时，其在不同叶片处存在着厚度要求，所以要对翼型厚度t进行约束，同时要保证优化后的翼型年发电量比原始翼型发电量高；基于步骤1的初始设计变量值计算设计变量的上下限，以DU99-W350翼型为例：

厚度约束条件0.345≤t≤0.355

翼型年发电量约束条件AEP_opt>AEP_ori

设计变量约束0.85x₀≤x≤1.15x₀

在追求极惯性矩最优的同时，会影响翼型的发电性能。翼型的发电量作为衡量翼型发电性能的一个重要指标，将优化后翼型的发电量作为另一个约束条件，要求优化后的翼型发电量不小于原始翼型发电量。

本实施例中，如图7所示，所述步骤6基于全信息自适应合作博弈方法进行设计变量和目标函数优化，具体包括：

S601：设置权重组数，假设共有Q组权重系数，q为当前的组数，权重系数的分配如下：

式中，w_ij代表权重系数。

S602：设置迭代次数，假设最大次数为M，k为当前迭代次数，随机生成初始策略S⁽⁰⁾。

S603：根据自适应策略集，对初始策略S⁽⁰⁾的设计变量划分给策略子集(表示第m个博弈方)，具体划分方法如下：

(1)对m个目标函数优化设计，得到m组最优解：

对应的策略集表示为：

(2)对任意设计变量x_j，对目标f_i的影响度为θ_ji表示为：

其中Δx_j为步长；

影响度归一化处理：

基于“空间距离”和“力矩”大小对策略空间进行划分，根据设计变量到目标函数之间的“距离空间”衡量远近，在策略子集划分中还要依据力矩的大小来划分，计算公式如下：

力矩的计算公式：

式中，d(i,j)表示设计变量到目标函数的空间距离；Mo(j)表示设计变量对于目标的力矩。

λ为力矩的阈值，定义如下:

(3)每个目标函数(博弈方)的所有设计变量的距离d(j,i)从最小到最大排序。具有相同距离的设计变量按影响度D(j,i)排序。影响度越大，排名越靠前。排名靠前的设计变量被划分给该博弈方，直到设计变量的力矩Mo(j)的累积值超过阈值λ。初步划分后，如果将相同的设计变量分配给多个博弈方，则该设计变量将以更大的影响度D(j,i)重新分配给博弈方。

步骤604：判断是否是第一次迭代，如果是第一次迭代，随机初始化为各博弈方策略子集；如果不是第一次迭代，进入步骤605，其中，表示第i个博弈方的第N个策略子集。

S605：基于通过模拟二进制交叉和多项式变异得到N个新策略

S606：计算2N个策略对应的最终成本函数并将前N个较小的最终成本函数对应的策略赋值给成本函数的计算方法如下：

式中，表示第i个博弈方在采用策略S_i时各个博弈方的相对成本，即目标函数f_i规格化的无量纲值。和分别表示第i个博弈方的初始策略和对应的补集。w_ij代表权系数，权系数取值大小代表各博弈方之间合作的程度，w_ij取值越大，则说明合作度越高。“*”表示计算两个矩阵的哈达玛积，其元素定义为两个矩阵对应元素的乘积。

S607：得到每个博弈方策略子集，即

S608：将得到的各博弈方策略方子集放在集合B中，即B＝B∪{S⁽¹⁾,S⁽²⁾,...,S^(N)}；更新迭代次数。

S609：判断是否对所有博弈方都进行了博弈，如果没有，则进入S604，否则进入S610。

S610：对迭代当前次数下的所有博弈方的策略子集进行整理，即更新下一次迭代次数，进入S603，如果是最后一次迭代，则进入S611。

S611：更新下一组权重系数的取值，执行S602；如果是最后一组权重系数，则进入S612。

S612：把集合B中的非劣解放入集合A中，即A中为全信息自适应合作博弈的最终解集。

如图5a，5b，5c所示，在pareto最优解集中选择三组翼型与原始翼型进行比较的结果。

实施例2

本实施例中，以DU99-W350翼型为原始翼型，该翼型最大厚度为0.35，叶片长度为61.5m，轮毂半径为1.5m，选取叶片展向处为24.05m的位置，扭角为9.011°，弦长为4.249，通过步骤1～步骤3，对翼型进行拟合和极惯性矩计算得到初始翼型的极惯性矩和年发电量分别为90.6939和20.1346，计算得出，通过步骤六中对全信息自适应合作博弈方法的解释，划分两个博弈方分别为f₁与f₂，本次博弈最大迭代次数20代，策略数量为100，策略保留数量为50，每一次迭代均在10组不同的权系数下进行，每进行一次迭代时，都会根据上次的收益去重新划分策略集。

第一次策略划分为：

S₂＝{x₂,x₅,x₈,x₁₀,x₁₁,x₁₂}；

经过20次迭代后，对计算出来的结果进行厚度和年发电量进行约束，得出符合条件的点，如图6所示。

综上，本发明提供的基于全信息自适应合作博弈理论的风力机翼型优化方法，前期使用parsec参数法在原始翼型上进行拟合，修改翼型的几何外形，运用全信息自适应合作博弈理论知识对12个设计变量和2个目标函数进行策略子集和策略空间的划分，通过运算，可以得到逼近pareto边界的优化解集，提高其抗颤振稳定性能并增强其发电性能(叶片截面发电量更高)。由经典抗颤振理论的研究，得出翼型的极惯性矩对其抗颤振性能具有突出的影响，其值可通过高斯积分计算得出。叶片的发电性能与叶片截面年发电量有关，其值可通过叶素动量理论(BEM)求得。将叶片截面最大年发电量和极惯性矩作为优化目标，运用全信息自适应博弈算法，优化得到的翼型比原有翼型具有更好的发电性能和抗颤振性能。

在本实施例中的其余技术特征，本领域技术人员均可以根据实际情况进行灵活选用以满足不同的具体实际需求。然而，对于本领域普通技术人员显而易见的是：不必采用这些特定细节来实行本发明。在其他实例中，为了避免混淆本发明，未具体描述公知的组成，结构或部件，均在本发明的权利要求书请求保护的技术方案限定技术保护范围之内。

本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。在以上描述中，为了提供对本发明的透彻理解，阐述了大量特定细节。然而，对于本领域普通技术人员显而易见的是：不必采用这些特定细节来实行本发明。在其他实例中，为了避免混淆本发明，未具体描述公知的技术，例如具体的施工细节，作业条件和其他的技术条件等。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，包括：

步骤1：基于parsec参数法在原始翼型上进行拟合，得到初始设计变量值x₀，修改翼型的几何外形；

步骤2：计算翼型的极惯性矩；

在极惯性矩计算过程中，首先建立极惯性矩与几何外形的关系，然后计算翼型的面积和质心坐标；平移坐标系，使坐标原点和质心重合，根据平行轴原理，求出翼型相对于质心轴系的惯性矩，得到翼型的极惯性矩；

步骤3：计算翼型年发电量；

步骤4：建立翼型约束条件：对翼型厚度t进行约束，并保证优化后的翼型年发电量比原始翼型发电量高，基于初始设计变量值x₀计算设计变量的上下限；即：

厚度约束条件t₁≤t≤t₂

翼型年发电量约束条件AEP_opt>AEP_ori

式中，t₁为翼型厚度的下限，t₂为翼型厚度的上限，AEP_opt为优化后的截面年发电量，AEP_ori为原截面年发电量；

步骤5：建立翼型的目标函数：基于步骤2和步骤3，设置目标函数为：

f₁(x)＝1/I

f₂(x)＝10⁶/AEP

其中，I为翼型的极惯性矩，AEP为年发电量；

步骤6：基于全信息自适应合作博弈方法，对设计变量和目标函数进行划分、收益函数计算、非劣解的筛选，最后得到一个pareto最优解集。

2.根据权利要求1所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述步骤1，基于parsec参数法在原始翼型上进行拟合，得到初始设计变量值x₀，修改翼型的几何外形，具体包括以下步骤：

S101、采取parsec参数法进行原始翼型上下表面拟合，得到翼型上、下表面拟合的形状函数：

上表面：

下表面：

上表面对应的系数的取值由下面方程组所确定：

下表面对应系数的取值由下面方程组所确定：

式中，x_te为后缘位置，x_te恒等于1；x_up为翼型上表面位置，x_lo为翼型下表面底部位置，r_le,up为翼型上表面前缘半径；r_le,lo为翼型下表面前缘半径，z_te为翼型后缘处纵坐标，Δz_te为翼型后缘厚度，z_up为翼型上表面上峰峰值，z_lo为翼型下表面下峰峰值，z_xx,up为翼型上表面曲率，z_xx,lo为翼型下表面曲率，α_te为翼型后缘方向，β_te为翼型后缘楔角；

得到设计变量：

x＝(r_le,up,r_le,lo,x_up,z_up,z_xx,up,x_lo,z_lo,z_xx,lo,Δz_te,z_te,α_te,β_te)；

S102、使用最小二乘法拟合得到原始翼型表面的12个设计变量值，得到修改后的翼型的几何外形。

3.根据权利要求2所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述步骤S102使用最小二乘法拟合得到原始翼型表面的12个设计变量值，得到修改后的翼型的几何外形，具体包括：

原始翼型表面的12个初始设计变量值的表达式为：

其中，z_u(i)代表翼型上翼面的第i个点的纵坐标，z_ou(i)代表原始翼型上翼面的第i个点的纵坐标，z_l(i)代表翼型下翼面的第i个点的纵坐标，z_ol(i)代表原始翼型下翼面的第i个点的纵坐标。

4.根据权利要求1所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述步骤2中计算翼型的面积和质心坐标点，具体包括：

求翼型的面积：

其中，S为翼型面积；c为翼型弦长，取值为标准值1；x代表横坐标，y代表纵坐标；代表翼型上翼面曲线函数，为翼型下翼面曲线函数；

翼型的质心坐标点：

其中，ρ(x,y)代表翼型不同坐标点处的密度，假设翼型材质是均匀，即ρ(x,y)当做常数，可进一步化简，即：

5.根据权利要求1所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述步骤2中得出翼型的极惯性矩，具体包括：

若将翼型前缘位置平移至坐标点处，在翼型坐标平面内任意一点的轴惯性矩表示为：

若将平面坐标系的原点平移到质心处，两坐标轴的方向保持不变，则对应的轴惯性矩为：

求得相对质心的惯性矩为：

I＝I_x1+I_y1。

6.根据权利要求1所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述步骤3，计算翼型截面年发电量，具体包括：

S301：设置初始轴向诱导因子a＝0，切向诱导因子b＝0；

S302：计算入流角φ和局部攻角α，计算公式如下：

α＝φ-θ

其中，λ_r为叶片展向r处的局部叶尖速比，B代表风力机的叶片数量，hub_r为风力机轮毂的半径，blade_r为叶片总长度，λ为叶片的叶尖速比，θ为叶片截面处翼型相对于叶尖翼型的扭角；

S303：用Xfoil软件计算该翼型在攻角α下的升力系数cl和阻力系数cd；

S304：计算出法向力系数C_n和切向力系数C_t，计算公式如下：

C_n＝cl cosφ+cd sinφ

C_t＝cl sinφ-cd cosφ

S305：计算此时的轴向诱导因子a，切向诱导因子b，计算公式如下：

其中，a_c＝0.2，σ′表示风轮实度，计算公式如下：

F表示损失系数，计算公式如下：

F＝F_tipF_hub

Prandtl叶尖修正因子

Prandtl叶根修正因子

S306：如果轴向诱导因子a和切向诱导因子b变化小于阈值δ＝1×10^-6，那么结束迭代，得到最终的轴向诱导因子a和切向诱导因子b，否则继续步骤302；

S307：计算出截面功率P，最后结合风速的威尔分布公式得出年发电量AEP：

7.根据权利要求1所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述步骤6基于全信息自适应合作博弈方法进行设计变量和目标函数优化，具体包括：

S601：设置权重组数，假设共有Q组权重系数，q为当前的组数，权重系数的分配如下：

式中，w_ij代表权重系数；

S602：设置迭代次数，假设最大次数为M，k为当前迭代次数，随机生成初始策略S⁽⁰⁾；

S603：根据自适应策略集，对初始策略S⁽⁰⁾的设计变量划分给策略子集

S604：基于k，判断是否是第一次迭代，如果是第一次迭代，随机初始化为各博弈方策略子集；如果不是第一次迭代，进入步骤605其中，表示第i个博弈方的第N个策略子集；

S605：基于通过模拟二进制交叉和多项式变异得到N个新策略

S606：计算2N个策略对应的最终成本函数并将前N个较小的最终成本函数对应的策略赋值给成本函数的计算方法如下：

式中，表示第i个博弈方在采用策略S_i时各个博弈方的相对成本，即目标函数f_i规格化的无量纲值；和分别表示第ⁱ个博弈方的初始策略和对应的补集；w_ij代表权系数，权系数取值大小代表各博弈方之间合作的程度，w_ij取值越大，则说明合作度越高；“*”表示计算两个矩阵的哈达玛积，其元素定义为两个矩阵对应元素的乘积；

S607：得到每个博弈方策略子集，即

S608：将得到的各博弈方策略方子集放在集合B中，即B＝B∪{S⁽¹⁾,S⁽²⁾,...,S^(N)}；更新迭代次数；

S609：判断是否对所有博弈方都进行了博弈，如果没有，则进入S604，否则进入S610；

S610：对迭代当前次数下的所有博弈方的策略子集进行整理，即更新下一次迭代次数，进入S603，如果是最后一次迭代，则进入S11；

S611：更新下一组权重系数的取值，执行S602；如果是最后一组权重系数，则进入S612；

S612：把集合B中的非劣解放入集合A中，即A中为全信息自适应合作博弈的最终解集。

8.根据权利要求7所述的基于全信息自适应合作博弈理论的风力机翼型优化方法，其特征在于，所述S603，根据自适应策略集，对初始策略S⁽⁰⁾的设计变量划分给策略子集中，具体划分方法如下：

(1)对m个目标函数优化设计，得到m组最优解：

对应的策略集表示为：

(2)对任意设计变量x_j，对目标f_i的影响度为θ_ji表示为：

其中，Δx_j为步长；

影响度归一化处理：

基于“空间距离”和“力矩”大小对策略空间进行划分，根据设计变量到目标函数之间的“距离空间”衡量远近，在策略子集划分中依据力矩的大小来划分，计算公式如下：

力矩的计算公式：

其中，d(j，i)表示设计变量到目标函数的空间距离；Mo(j)表示设计变量对于目标的力矩；

λ为力矩的阈值，定义如下：

(3)将每个目标函数即博弈方的所有设计变量的距离d(j,i)从最小到最大排序；具有相同距离的设计变量按影响度D(j,i)排序；影响度越大，排名越靠前；排名靠前的设计变量被划分给该博弈方，直到设计变量的力矩Mo(j)的累积值超过阈值λ；初步划分后，如果将相同的设计变量分配给多个博弈方，则该设计变量以更大的影响度D(j,i)重新分配给博弈方。