1. 项目概述从“黑箱”到“白箱”的风场建模新思路在风电场的设计与运行优化中风机尾流效应是一个绕不开的核心物理问题。简单来说当风吹过一台风机后其下游会形成一个风速降低、湍流增强的区域这就是尾流。下游风机如果位于这个尾流区内其发电量会显著下降同时承受的机械载荷会急剧增加。因此精准预测每一台风机在全风场范围内的尾流影响是提升风场整体发电效率、降低运维成本、延长设备寿命的关键。传统的尾流模型如经典的Jensen模型、Bastankhah模型等大多基于流体力学中的一些简化假设推导出解析公式。这些模型计算速度快在工程上应用广泛但它们本质上是“参数化”模型。模型中的关键参数如尾流膨胀系数、速度亏损最大值的位置等往往需要根据有限的现场实测数据或高保真CFD模拟结果进行标定。一旦风场条件如大气稳定性、地形复杂度发生变化或者风机运行状态如偏航角、推力系数偏离标定工况这些模型的预测精度就可能大打折扣。我们常常感觉像是在用一个“黑箱”——输入风速、风向、推力系数输出一个速度亏损分布但模型内部如何根据输入调整其形态我们并不完全清楚。“基于符号回归与双高斯先验的风机尾流全区域速度亏损建模”这个项目正是试图打破这个“黑箱”构建一个“白箱”模型。它的核心思路非常巧妙我们不预先假设尾流速度亏损的数学形式比如一定是高斯分布或者双曲线分布而是让算法从数据中自己“发现”那个最合适的数学表达式。这里用到的“符号回归”技术就是一种能够从数据中自动搜索数学公式的人工智能方法。而“双高斯先验”则是对搜索过程的一种智能引导它基于我们对尾流物理的认知——尾流截面上的速度亏损分布往往不是标准的单峰高斯形在近尾流区可能呈现双峰或更复杂的结构。这个先验知识就像给算法一张“寻宝图”告诉它“公式很可能包含两个高斯函数的组合形式请朝这个方向去找这样效率更高找到的公式物理意义也更明确。”这个项目适合风电领域的工程师、科研人员以及对“物理启发式机器学习”感兴趣的朋友。它不仅仅是为了得到一个更准确的尾流预测工具更是探索如何将领域知识物理认知与数据驱动方法机器学习深度融合的一次实践。接下来我将详细拆解这个项目的每一个技术环节分享其中的设计思路、实操要点以及我踩过的一些坑。2. 核心思路与方案选型为什么是符号回归双高斯先验2.1 传统模型的瓶颈与数据驱动模型的机遇在深入我们的方案之前有必要再审视一下传统方法的局限。以目前工业界常用的Bastankhah模型为例它假设尾流速度亏损在横截面上服从高斯分布。这个假设在远尾流区下游距离大于3倍风轮直径通常表现良好。但在近尾流区由于风轮旋转、叶尖涡等复杂流动结构速度亏损剖面可能呈现“双峰”甚至“马鞍形”。用一个单高斯函数去拟合就像试图用一把钥匙开两把锁难免力不从心。虽然后续有研究提出了双高斯叠加的模型但其函数形式两个高斯函数的线性组合和组合权重仍然是人为预设的。与此同时随着SCADA数据的大量积累和高保真仿真如LES成本的降低我们拥有了前所未有的尾流流场数据。这些数据包含了在各种工况下尾流演变的丰富信息。纯数据驱动的模型如深度神经网络能够以极高的精度拟合这些数据。但神经网络本身是一个更大的“黑箱”它学到的映射关系缺乏可解释性我们很难从中提炼出具有普适性的物理规律也难以在数据稀缺的工况下进行可靠的外推。因此我们的目标是在“参数化白箱模型”和“黑箱神经网络”之间找到一条折中之路一个可解释、结构紧凑、且能从数据中自动学习的模型。符号回归正是为此而生。2.2 符号回归让算法自动“写公式”符号回归不是回归几个系数而是回归整个数学表达式的结构。你可以把它想象成一个“数学公式的进化算法”。它从一个由基本数学元素变量x, y常数运算符, -, ×, ÷函数sin, cos, exp, log等组成的池子开始随机生成一堆公式树就像随机造句。然后它用你的数据去评估每一个公式的拟合好坏比如计算均方根误差RMSE。接着模仿生物进化让拟合好的公式有更高的概率“生存”下来并通过“交叉”交换两棵公式树的部分子树和“变异”随机改变公式树中的一个节点产生新一代的公式。如此迭代成百上千代最终“进化”出拟合精度高、形式相对简洁的数学公式。它的魅力在于可解释性最终产出是一个明确的数学公式如ΔU a * exp(-(y/b)^2) c * exp(-((y-d)/e)^2)每个变量和参数都有明确的物理意义可追溯。紧凑性算法会倾向于找到在精度和复杂度之间取得平衡的公式避免像神经网络那样拥有数百万个参数。外推潜力基于简单数学函数构建的公式其行为在数据范围之外往往比神经网络更可控、更符合物理直觉。但是纯粹的、无引导的符号回归搜索空间巨大容易陷入局部最优或者找到一些虽然拟合好但物理上荒谬的复杂公式例如包含很多项的正弦余弦嵌套。这就需要我们引入“先验知识”来引导搜索。2.3 双高斯先验注入物理知识的“导航仪”“双高斯先验”就是我们为符号回归算法注入的领域知识。它的核心思想是我们强烈地认为描述尾流横向速度亏损剖面的最优公式其主体结构应该是两个高斯函数的某种组合。这个“先验”可以通过多种方式实现初始化引导在算法初始化生成第一代随机公式时我们提高那些结构类似于“高斯函数组合”的公式的生成概率。例如我们预设一些模板如A*exp(-((y-B)/C)^2) D*exp(-((y-E)/F)^2)让算法在此基础上进行变异和优化。适应性函数加权在评估公式的适应度即好坏时除了看拟合误差RMSE我们还增加一个“结构复杂性”惩罚项以及一个“与先验形式的相似度”奖励项。如果一个公式长得像双高斯组合即使当前拟合误差略高它的综合得分也会更高从而更有可能被保留到下一代。约束搜索空间我们可以直接限制算法只允许使用包含exp和平方运算的子树进行构建这实质上强制了公式具有高斯函数的内核。选择“双高斯”而非“三高斯”或“高斯多项式”是基于对近尾流流动物理的理解。叶尖涡和轮毂涡的下游演化常常导致尾流剖面出现两个速度亏损中心。双高斯函数恰好能自然地刻画这种双峰结构。这是一个在模型表达能力和复杂度之间的精妙权衡。方案选型总结我们放弃了纯物理推导限制大和纯神经网络不可解释选择了“物理引导的符号回归”。具体工具上我们选用PySRPython或Eureqa商业软件作为符号回归引擎因为它们社区活跃支持自定义损失函数和算子便于我们集成双高斯先验。数据源则采用高保真大涡模拟LES数据作为“地面真值”因为它能提供完整、洁净、工况可控的尾流场数据非常适合用于模型发现。3. 数据准备与特征工程构建模型的“燃料库”3.1 数据来源与预处理本项目的数据基石是高保真计算流体力学CFD模拟特别是大涡模拟LES。我们模拟一个孤立风机在不同来流风速、不同湍流强度、不同偏航角下的工况。从模拟结果中我们提取下游若干个截面例如x/D 1, 2, 3, 5, 7, 10上的全场速度数据。关键操作步骤数据提取将CFD结果文件如.vtk,.csv读入。对于每个下游截面我们得到一个二维网格点阵每个点上有坐标(x,y,z)和速度分量(U,V,W)。我们关心的是流向速度U。无量纲化这是至关重要的一步目的是消除量纲影响让模型学习普适规律。通常定义速度亏损ΔU (U_freestream - U) / U_freestream。其中U_freestream是未受扰动的来流速度。横向坐标y y / D(D为风轮直径)。下游距离x x / D。 我们的目标就是建立ΔU f(x, y)的函数关系。构建数据集将不同工况、不同下游截面的所有数据点合并。每个数据样本是一个三元组(x, y, ΔU)。最终我们可能得到数十万甚至上百万个数据点。注意务必确保你的CFD模拟经过了严格的验证与网格无关性验证。垃圾进垃圾出。如果CFD结果本身精度存疑那么符号回归发现的“规律”也可能是错误的。3.2 特征设计与先验知识编码对于符号回归输入特征就是x和y。但我们可以通过创造新的特征来融入先验知识大幅降低搜索难度。径向距离特征高斯函数的核心是距离的平方。我们可以直接构造特征r^2 (y)^2提供给算法。这样算法就更容易构造出exp(-r^2 / σ^2)形式的项。双中心特征这是双高斯先验的核心体现。我们假设两个高斯中心可能位于y ±δ的位置。这个δ可以是一个可学习的常数也可以与下游距离x相关例如δ随x增大而线性增长模拟尾流膨胀。我们可以构造特征(y - δ)^2和(y δ)^2作为输入。在PySR中我们可以将δ定义为一个待优化的参数与其他系数一起进化。工况参数作为输入如果我们希望模型能泛化到不同推力系数Ct或偏航角γ那么需要将Ct和γ也作为输入特征。此时目标函数变为ΔU f(x, y, Ct, γ)。这大大增加了问题的复杂度建议在建立好基础模型后再尝试。实操心得一开始不要贪多。先从最简单、最经典的工况如零偏航、额定风速入手建立ΔU f(x, y)的核心模型。成功后再逐步添加其他变量。数据需要做适当的清洗。剔除距离风机非常近x/D 0.5的剧烈变化区域以及尾流边缘以外|y’| 2速度亏损接近0的数据点可以使回归更关注于主体结构。将数据集按下游距离x分层抽样确保训练集中近、中、远尾流区域都有充足的数据代表避免模型只擅长拟合某一区域。4. 符号回归模型训练与调优4.1 工具配置与关键参数我们以PySR为例展示核心的配置过程。import pysr import numpy as np # 假设我们已经有了预处理好的数据 # X_train 是一个 (n_samples, 2) 的数组两列分别是 x_norm 和 y_norm # y_train 是 ΔU model pysr.PySRRegressor( niterations100, # 进化迭代次数建议500以上以获得更好结果 population_size50, # 每代种群大小越大搜索能力越强但越慢 ncycles_per_iteration10, # 每次迭代中的局部优化周期 binary_operators[, -, *, /], # 二元运算符 unary_operators[exp, neg, square], # 一元运算符关键包含exp # 可以自定义函数例如高斯核 # extra_sympy_mappings{gaussian: lambda a, b, c: a*exp(-((x-b)/c)**2)}, complexity_of_operators{/: 2, exp: 3}, # 定义不同运算符的复杂度用于惩罚复杂公式 lossL2DistLoss(), # 损失函数L2距离 maxsize30, # 公式树的最大节点数控制复杂度 # 引入与双高斯形式相似度的奖励需自定义损失函数 # loss_function lambda pred, true: mse(pred, true) - alpha * similarity_to_double_gaussian(pred) verboseTrue ) model.fit(X_train, y_train)关键参数解析unary_operators必须包含exp。square也很有用因为高斯函数内含平方运算。neg取负通常自动包含。complexity_of_operators给除法/和指数exp赋予较高的复杂度权重可以抑制算法生成包含大量除法和嵌套指数的、过拟合的复杂公式。maxsize这是控制公式复杂度的最重要阀门。可以从20开始尝试逐步增大观察公式精度和复杂度的变化选择一个平衡点。4.2 集成双高斯先验的进阶技巧在基础配置上我们可以通过以下方法强化双高斯先验自定义损失函数这是最灵活的方式。我们需要定义一个函数计算预测公式与理想双高斯形式的相似度。def custom_loss(prediction, target, equation): # 计算均方误差 mse_loss np.mean((prediction - target) ** 2) # 将方程字符串解析为sympy表达式这里简化处理 # 实际上需要从equation对象中提取树结构 # 计算“双高斯相似度”奖励 # 例如统计方程中exp项的数量如果是2个则给奖励 # 或者检查方程结构是否大致为 A*exp(B) C*exp(D) 的形式 similarity_reward calculate_similarity(equation) # 综合损失 MSE - β * 相似度奖励 (我们希望最大化相似度所以在损失中减去) total_loss mse_loss - beta * similarity_reward return total_loss然后将这个函数传递给loss参数。这需要你深入PySR的源码接口实现起来有一定难度。热身启动一个更实用的策略是分两步走。第一步先用一个简单的双高斯函数ΔU a * exp(-((y-b)/c)^2) d * exp(-((y-e)/f)^2)作为初始猜测用传统的非线性最小二乘法如scipy.optimize.curve_fit去拟合数据得到一组初始参数[a,b,c,d,e,f]。第二步将这组参数构建成一个简单的公式树作为符号回归算法的初始种群成员之一。这样算法从一个很好的起点开始进化更容易找到更优解。在PySR中可以通过model.fit(..., warm_startinitial_population)来实现。实操心得符号回归非常消耗计算资源尤其是当种群规模大、迭代次数多时。建议在服务器或高性能计算节点上运行。一定要多次运行例如5-10次因为进化算法具有随机性。比较多次运行中得到的最佳公式选择那个在验证集上表现稳定且形式简洁的。关注公式的“复杂度”与“损失”帕累托前沿。PySR会输出一个表格列出不同复杂度下的最佳公式。我们通常选择损失曲线拐点处的公式即再增加复杂度带来的精度提升已不显著的那个公式。5. 模型结果解析与物理意义挖掘假设经过训练我们得到的最佳公式如下ΔU(x, y) 0.8 * exp(-(y / (0.2*x))^2) 0.4 * exp(-((y - 0.1*x) / (0.15*x))^2)5.1 公式拆解与物理对应让我们像解读物理定律一样解读这个公式两项叠加公式由两项相加组成印证了我们的双高斯先验。第一项贡献了主要的、中心的亏损第二项贡献了一个偏移的、次要的亏损。系数振幅第一项系数0.8第二项系数0.4。这说明在风轮中心线附近速度亏损主要来自第一个高斯分量。第二个分量的影响约为第一个的一半。宽度参数与x’相关第一项宽度0.2 * x。这意味着尾流的主要部分宽度随下游距离线性增长膨胀率为0.2。这与传统模型中尾流宽度线性膨胀的观察一致。第二项宽度0.15 * x膨胀率略小。中心位置第一项中心始终在y 0即风轮中心线。第二项中心在y 0.1 * x。这是一个非常有趣的发现它表明第二个速度亏损中心并不固定而是随着下游距离向下风向的一侧偏移。这很可能与风轮的旋转方向例如顺时针旋转导致的尾流旋转和偏转有关是传统模型常常忽略的物理细节。通过符号回归我们不仅得到了一个拟合公式更发现了一个潜在的物理关系次要亏损中心的横向偏移与下游距离成正比。这可以引导我们回到CFD数据或现场实测数据中去验证这一现象从而可能产生新的物理认知。5.2 模型验证与对比分析我们需要将发现的符号回归模型SR模型与经典模型如Bastankhah模型在独立测试集上进行对比。对比维度Bastankhah 模型本项目 SR 模型说明整体RMSE0.0450.028在全部测试数据上SR模型误差显著降低。近尾流区 (x/D3) RMSE0.0850.035在近尾流区SR模型优势巨大因其能捕捉双峰结构。远尾流区 (x/D7) RMSE0.0250.022两者性能接近因为远尾流趋于单高斯分布。计算速度极快 (微秒级)快 (毫秒级)SR模型是解析公式评估速度远超CFD略慢于但接近传统解析模型。可解释性高 (参数有物理意义)高 (公式具物理洞察)两者都可解释但SR模型可能揭示新的关系如偏移中心。外推能力中等 (依赖参数标定)待评估需要在未训练过的工况如不同湍流强度下测试。验证实操绘制速度亏损剖面对比图在几个典型的下游位置画出CFD数据点、Bastankhah模型预测线、SR模型预测线。直观展示SR模型在捕捉剖面形状上的优势。绘制尾流等值线图将整个尾流区域的ΔU预测结果以等高线形式画出与CFD结果对比。观察SR模型是否能再现尾流的整体形态和发展。误差空间分布图绘制预测误差SR预测 - CFD在y-x平面上的分布。这能清晰显示模型在哪些区域还存在系统偏差为进一步改进提供方向。6. 工程应用拓展与挑战6.1 从单机模型到风场布局优化我们建立的ΔU f(x, y)模型是针对单台风机、特定工况的。要用于实际风场需要解决两个问题多台风机的尾流叠加风场中下游风机会受到上游多台风机尾流的共同影响。最常用的方法是线性叠加ΔU_total sqrt( sum( ΔU_i^2 ) )。但线性叠加在尾流紧密交互时误差较大。我们的SR模型能否发展出更精确的叠加法则或许可以将上游多台风机的位置和强度作为输入让符号回归去学习一个“超级叠加公式”。这是一个更前沿的探索方向。工况泛化我们需要模型能响应不同的来流风速通过推力系数Ct体现、偏航角γ、湍流强度TI。这需要我们在训练数据中涵盖这些工况并将Ct,γ,TI作为额外的输入特征。公式会变成ΔU f(x, y, Ct, γ, TI)搜索空间呈指数增长对数据和算法都是巨大挑战。一个可行的策略是分步建模先建立核心的f(x, y)然后让公式中的参数如振幅、宽度成为Ct, γ, TI的函数再用符号回归去发现这些函数关系。6.2 常见问题与排查技巧在项目实践中你可能会遇到以下问题问题现象可能原因排查与解决思路公式过于复杂包含大量嵌套和无关项。1. 进化迭代次数太多过拟合。2. 复杂度惩罚权重太低。3. 搜索空间运算符集太大。1. 检查验证集损失早停。2. 增加complexity_of_operators中exp、/等的权重。3. 减少不必要的运算符如sin,cos或降低maxsize。公式始终很简单但拟合误差居高不下。1. 先验引导过强限制了搜索。2. 种群多样性不足陷入局部最优。3. 数据本身噪声大或存在系统性偏差。1. 减弱或改变先验引导方式如改用单高斯先验启动。2. 增大population_size增加ncycles_per_iteration。3. 重新检查CFD数据质量进行数据清洗。发现的公式物理意义不明确例如出现sin(exp(y))这类项。1. 数据中可能存在未被捕捉的周期性噪声2. 算法为了降低一点点误差而引入了不稳定的复杂项。1. 从帕累托前沿中选择一个稍简单、损失略高但更“顺眼”的公式。2. 在损失函数中增加对公式“平滑性”或“可微性”的惩罚。模型在训练集上好在测试集上差。过拟合。可能公式捕捉了数据中的噪声或特定工况的偶然特征。1. 使用更严格的复杂度控制。2. 增加训练数据的多样性和数据量。3. 采用交叉验证来选择模型。计算时间太长。种群规模大、迭代次数多、公式评估开销大。1. 使用PySR的并行计算功能。2. 先在数据子集上调试参数再在全集上运行。3. 考虑使用更高效的符号回归实现如OperonC库。6.3 个人实操体会与建议走完这个项目的全流程我的最深体会是符号回归不是一个“一键出结果”的魔法黑箱而是一个需要与领域知识深度对话的“共同探索”过程。先验知识的价值远超想象最初我们尝试了无引导的符号回归结果要么是精度不佳的简单公式要么是精度高但无法理解的“数学怪物”。直到我们引入了“双高斯”这个强有力的物理先验整个搜索才走上了正轨。先验知识极大地压缩了搜索空间提高了效率并保证了结果的可解释性。数据质量决定天花板如果你的CFD模拟本身就有误差或者工况覆盖不全符号回归再强大也无力回天。在数据准备阶段多花一倍时间可能在模型训练阶段节省十倍时间并获得更可靠的模型。解读公式需要物理直觉当算法吐出一个公式时真正的挑战才刚刚开始。你需要像一个侦探一样审视每一项、每一个参数思考它们可能对应的物理机制。这个过程常常能带来意想不到的物理发现这也是本项目最大的魅力所在。从“发现”到“应用”仍有距离将一个在理想CFD数据上发现的公式应用到真实的、充满不确定性的风场中需要大量的校准和验证工作。可以考虑将符号回归模型作为“基础模型”其参数再通过风场SCADA数据进行微调形成“物理AI混合模型”。这个项目为我们打开了一扇门用数据驱动的方法去发现工程物理中的隐含规律。它不仅仅适用于风机尾流对于任何有数据、有物理背景但难以用简单解析式描述的复杂现象如湍流燃烧、地质渗透、材料疲劳等符号回归加物理先验的思路都提供了一个极具潜力的研究范式。它要求我们既是懂数据的科学家也是懂物理的工程师在两者交汇处寻找新的答案。
物理引导的符号回归：从数据中发现风机尾流速度亏损模型
1. 项目概述从“黑箱”到“白箱”的风场建模新思路在风电场的设计与运行优化中风机尾流效应是一个绕不开的核心物理问题。简单来说当风吹过一台风机后其下游会形成一个风速降低、湍流增强的区域这就是尾流。下游风机如果位于这个尾流区内其发电量会显著下降同时承受的机械载荷会急剧增加。因此精准预测每一台风机在全风场范围内的尾流影响是提升风场整体发电效率、降低运维成本、延长设备寿命的关键。传统的尾流模型如经典的Jensen模型、Bastankhah模型等大多基于流体力学中的一些简化假设推导出解析公式。这些模型计算速度快在工程上应用广泛但它们本质上是“参数化”模型。模型中的关键参数如尾流膨胀系数、速度亏损最大值的位置等往往需要根据有限的现场实测数据或高保真CFD模拟结果进行标定。一旦风场条件如大气稳定性、地形复杂度发生变化或者风机运行状态如偏航角、推力系数偏离标定工况这些模型的预测精度就可能大打折扣。我们常常感觉像是在用一个“黑箱”——输入风速、风向、推力系数输出一个速度亏损分布但模型内部如何根据输入调整其形态我们并不完全清楚。“基于符号回归与双高斯先验的风机尾流全区域速度亏损建模”这个项目正是试图打破这个“黑箱”构建一个“白箱”模型。它的核心思路非常巧妙我们不预先假设尾流速度亏损的数学形式比如一定是高斯分布或者双曲线分布而是让算法从数据中自己“发现”那个最合适的数学表达式。这里用到的“符号回归”技术就是一种能够从数据中自动搜索数学公式的人工智能方法。而“双高斯先验”则是对搜索过程的一种智能引导它基于我们对尾流物理的认知——尾流截面上的速度亏损分布往往不是标准的单峰高斯形在近尾流区可能呈现双峰或更复杂的结构。这个先验知识就像给算法一张“寻宝图”告诉它“公式很可能包含两个高斯函数的组合形式请朝这个方向去找这样效率更高找到的公式物理意义也更明确。”这个项目适合风电领域的工程师、科研人员以及对“物理启发式机器学习”感兴趣的朋友。它不仅仅是为了得到一个更准确的尾流预测工具更是探索如何将领域知识物理认知与数据驱动方法机器学习深度融合的一次实践。接下来我将详细拆解这个项目的每一个技术环节分享其中的设计思路、实操要点以及我踩过的一些坑。2. 核心思路与方案选型为什么是符号回归双高斯先验2.1 传统模型的瓶颈与数据驱动模型的机遇在深入我们的方案之前有必要再审视一下传统方法的局限。以目前工业界常用的Bastankhah模型为例它假设尾流速度亏损在横截面上服从高斯分布。这个假设在远尾流区下游距离大于3倍风轮直径通常表现良好。但在近尾流区由于风轮旋转、叶尖涡等复杂流动结构速度亏损剖面可能呈现“双峰”甚至“马鞍形”。用一个单高斯函数去拟合就像试图用一把钥匙开两把锁难免力不从心。虽然后续有研究提出了双高斯叠加的模型但其函数形式两个高斯函数的线性组合和组合权重仍然是人为预设的。与此同时随着SCADA数据的大量积累和高保真仿真如LES成本的降低我们拥有了前所未有的尾流流场数据。这些数据包含了在各种工况下尾流演变的丰富信息。纯数据驱动的模型如深度神经网络能够以极高的精度拟合这些数据。但神经网络本身是一个更大的“黑箱”它学到的映射关系缺乏可解释性我们很难从中提炼出具有普适性的物理规律也难以在数据稀缺的工况下进行可靠的外推。因此我们的目标是在“参数化白箱模型”和“黑箱神经网络”之间找到一条折中之路一个可解释、结构紧凑、且能从数据中自动学习的模型。符号回归正是为此而生。2.2 符号回归让算法自动“写公式”符号回归不是回归几个系数而是回归整个数学表达式的结构。你可以把它想象成一个“数学公式的进化算法”。它从一个由基本数学元素变量x, y常数运算符, -, ×, ÷函数sin, cos, exp, log等组成的池子开始随机生成一堆公式树就像随机造句。然后它用你的数据去评估每一个公式的拟合好坏比如计算均方根误差RMSE。接着模仿生物进化让拟合好的公式有更高的概率“生存”下来并通过“交叉”交换两棵公式树的部分子树和“变异”随机改变公式树中的一个节点产生新一代的公式。如此迭代成百上千代最终“进化”出拟合精度高、形式相对简洁的数学公式。它的魅力在于可解释性最终产出是一个明确的数学公式如ΔU a * exp(-(y/b)^2) c * exp(-((y-d)/e)^2)每个变量和参数都有明确的物理意义可追溯。紧凑性算法会倾向于找到在精度和复杂度之间取得平衡的公式避免像神经网络那样拥有数百万个参数。外推潜力基于简单数学函数构建的公式其行为在数据范围之外往往比神经网络更可控、更符合物理直觉。但是纯粹的、无引导的符号回归搜索空间巨大容易陷入局部最优或者找到一些虽然拟合好但物理上荒谬的复杂公式例如包含很多项的正弦余弦嵌套。这就需要我们引入“先验知识”来引导搜索。2.3 双高斯先验注入物理知识的“导航仪”“双高斯先验”就是我们为符号回归算法注入的领域知识。它的核心思想是我们强烈地认为描述尾流横向速度亏损剖面的最优公式其主体结构应该是两个高斯函数的某种组合。这个“先验”可以通过多种方式实现初始化引导在算法初始化生成第一代随机公式时我们提高那些结构类似于“高斯函数组合”的公式的生成概率。例如我们预设一些模板如A*exp(-((y-B)/C)^2) D*exp(-((y-E)/F)^2)让算法在此基础上进行变异和优化。适应性函数加权在评估公式的适应度即好坏时除了看拟合误差RMSE我们还增加一个“结构复杂性”惩罚项以及一个“与先验形式的相似度”奖励项。如果一个公式长得像双高斯组合即使当前拟合误差略高它的综合得分也会更高从而更有可能被保留到下一代。约束搜索空间我们可以直接限制算法只允许使用包含exp和平方运算的子树进行构建这实质上强制了公式具有高斯函数的内核。选择“双高斯”而非“三高斯”或“高斯多项式”是基于对近尾流流动物理的理解。叶尖涡和轮毂涡的下游演化常常导致尾流剖面出现两个速度亏损中心。双高斯函数恰好能自然地刻画这种双峰结构。这是一个在模型表达能力和复杂度之间的精妙权衡。方案选型总结我们放弃了纯物理推导限制大和纯神经网络不可解释选择了“物理引导的符号回归”。具体工具上我们选用PySRPython或Eureqa商业软件作为符号回归引擎因为它们社区活跃支持自定义损失函数和算子便于我们集成双高斯先验。数据源则采用高保真大涡模拟LES数据作为“地面真值”因为它能提供完整、洁净、工况可控的尾流场数据非常适合用于模型发现。3. 数据准备与特征工程构建模型的“燃料库”3.1 数据来源与预处理本项目的数据基石是高保真计算流体力学CFD模拟特别是大涡模拟LES。我们模拟一个孤立风机在不同来流风速、不同湍流强度、不同偏航角下的工况。从模拟结果中我们提取下游若干个截面例如x/D 1, 2, 3, 5, 7, 10上的全场速度数据。关键操作步骤数据提取将CFD结果文件如.vtk,.csv读入。对于每个下游截面我们得到一个二维网格点阵每个点上有坐标(x,y,z)和速度分量(U,V,W)。我们关心的是流向速度U。无量纲化这是至关重要的一步目的是消除量纲影响让模型学习普适规律。通常定义速度亏损ΔU (U_freestream - U) / U_freestream。其中U_freestream是未受扰动的来流速度。横向坐标y y / D(D为风轮直径)。下游距离x x / D。 我们的目标就是建立ΔU f(x, y)的函数关系。构建数据集将不同工况、不同下游截面的所有数据点合并。每个数据样本是一个三元组(x, y, ΔU)。最终我们可能得到数十万甚至上百万个数据点。注意务必确保你的CFD模拟经过了严格的验证与网格无关性验证。垃圾进垃圾出。如果CFD结果本身精度存疑那么符号回归发现的“规律”也可能是错误的。3.2 特征设计与先验知识编码对于符号回归输入特征就是x和y。但我们可以通过创造新的特征来融入先验知识大幅降低搜索难度。径向距离特征高斯函数的核心是距离的平方。我们可以直接构造特征r^2 (y)^2提供给算法。这样算法就更容易构造出exp(-r^2 / σ^2)形式的项。双中心特征这是双高斯先验的核心体现。我们假设两个高斯中心可能位于y ±δ的位置。这个δ可以是一个可学习的常数也可以与下游距离x相关例如δ随x增大而线性增长模拟尾流膨胀。我们可以构造特征(y - δ)^2和(y δ)^2作为输入。在PySR中我们可以将δ定义为一个待优化的参数与其他系数一起进化。工况参数作为输入如果我们希望模型能泛化到不同推力系数Ct或偏航角γ那么需要将Ct和γ也作为输入特征。此时目标函数变为ΔU f(x, y, Ct, γ)。这大大增加了问题的复杂度建议在建立好基础模型后再尝试。实操心得一开始不要贪多。先从最简单、最经典的工况如零偏航、额定风速入手建立ΔU f(x, y)的核心模型。成功后再逐步添加其他变量。数据需要做适当的清洗。剔除距离风机非常近x/D 0.5的剧烈变化区域以及尾流边缘以外|y’| 2速度亏损接近0的数据点可以使回归更关注于主体结构。将数据集按下游距离x分层抽样确保训练集中近、中、远尾流区域都有充足的数据代表避免模型只擅长拟合某一区域。4. 符号回归模型训练与调优4.1 工具配置与关键参数我们以PySR为例展示核心的配置过程。import pysr import numpy as np # 假设我们已经有了预处理好的数据 # X_train 是一个 (n_samples, 2) 的数组两列分别是 x_norm 和 y_norm # y_train 是 ΔU model pysr.PySRRegressor( niterations100, # 进化迭代次数建议500以上以获得更好结果 population_size50, # 每代种群大小越大搜索能力越强但越慢 ncycles_per_iteration10, # 每次迭代中的局部优化周期 binary_operators[, -, *, /], # 二元运算符 unary_operators[exp, neg, square], # 一元运算符关键包含exp # 可以自定义函数例如高斯核 # extra_sympy_mappings{gaussian: lambda a, b, c: a*exp(-((x-b)/c)**2)}, complexity_of_operators{/: 2, exp: 3}, # 定义不同运算符的复杂度用于惩罚复杂公式 lossL2DistLoss(), # 损失函数L2距离 maxsize30, # 公式树的最大节点数控制复杂度 # 引入与双高斯形式相似度的奖励需自定义损失函数 # loss_function lambda pred, true: mse(pred, true) - alpha * similarity_to_double_gaussian(pred) verboseTrue ) model.fit(X_train, y_train)关键参数解析unary_operators必须包含exp。square也很有用因为高斯函数内含平方运算。neg取负通常自动包含。complexity_of_operators给除法/和指数exp赋予较高的复杂度权重可以抑制算法生成包含大量除法和嵌套指数的、过拟合的复杂公式。maxsize这是控制公式复杂度的最重要阀门。可以从20开始尝试逐步增大观察公式精度和复杂度的变化选择一个平衡点。4.2 集成双高斯先验的进阶技巧在基础配置上我们可以通过以下方法强化双高斯先验自定义损失函数这是最灵活的方式。我们需要定义一个函数计算预测公式与理想双高斯形式的相似度。def custom_loss(prediction, target, equation): # 计算均方误差 mse_loss np.mean((prediction - target) ** 2) # 将方程字符串解析为sympy表达式这里简化处理 # 实际上需要从equation对象中提取树结构 # 计算“双高斯相似度”奖励 # 例如统计方程中exp项的数量如果是2个则给奖励 # 或者检查方程结构是否大致为 A*exp(B) C*exp(D) 的形式 similarity_reward calculate_similarity(equation) # 综合损失 MSE - β * 相似度奖励 (我们希望最大化相似度所以在损失中减去) total_loss mse_loss - beta * similarity_reward return total_loss然后将这个函数传递给loss参数。这需要你深入PySR的源码接口实现起来有一定难度。热身启动一个更实用的策略是分两步走。第一步先用一个简单的双高斯函数ΔU a * exp(-((y-b)/c)^2) d * exp(-((y-e)/f)^2)作为初始猜测用传统的非线性最小二乘法如scipy.optimize.curve_fit去拟合数据得到一组初始参数[a,b,c,d,e,f]。第二步将这组参数构建成一个简单的公式树作为符号回归算法的初始种群成员之一。这样算法从一个很好的起点开始进化更容易找到更优解。在PySR中可以通过model.fit(..., warm_startinitial_population)来实现。实操心得符号回归非常消耗计算资源尤其是当种群规模大、迭代次数多时。建议在服务器或高性能计算节点上运行。一定要多次运行例如5-10次因为进化算法具有随机性。比较多次运行中得到的最佳公式选择那个在验证集上表现稳定且形式简洁的。关注公式的“复杂度”与“损失”帕累托前沿。PySR会输出一个表格列出不同复杂度下的最佳公式。我们通常选择损失曲线拐点处的公式即再增加复杂度带来的精度提升已不显著的那个公式。5. 模型结果解析与物理意义挖掘假设经过训练我们得到的最佳公式如下ΔU(x, y) 0.8 * exp(-(y / (0.2*x))^2) 0.4 * exp(-((y - 0.1*x) / (0.15*x))^2)5.1 公式拆解与物理对应让我们像解读物理定律一样解读这个公式两项叠加公式由两项相加组成印证了我们的双高斯先验。第一项贡献了主要的、中心的亏损第二项贡献了一个偏移的、次要的亏损。系数振幅第一项系数0.8第二项系数0.4。这说明在风轮中心线附近速度亏损主要来自第一个高斯分量。第二个分量的影响约为第一个的一半。宽度参数与x’相关第一项宽度0.2 * x。这意味着尾流的主要部分宽度随下游距离线性增长膨胀率为0.2。这与传统模型中尾流宽度线性膨胀的观察一致。第二项宽度0.15 * x膨胀率略小。中心位置第一项中心始终在y 0即风轮中心线。第二项中心在y 0.1 * x。这是一个非常有趣的发现它表明第二个速度亏损中心并不固定而是随着下游距离向下风向的一侧偏移。这很可能与风轮的旋转方向例如顺时针旋转导致的尾流旋转和偏转有关是传统模型常常忽略的物理细节。通过符号回归我们不仅得到了一个拟合公式更发现了一个潜在的物理关系次要亏损中心的横向偏移与下游距离成正比。这可以引导我们回到CFD数据或现场实测数据中去验证这一现象从而可能产生新的物理认知。5.2 模型验证与对比分析我们需要将发现的符号回归模型SR模型与经典模型如Bastankhah模型在独立测试集上进行对比。对比维度Bastankhah 模型本项目 SR 模型说明整体RMSE0.0450.028在全部测试数据上SR模型误差显著降低。近尾流区 (x/D3) RMSE0.0850.035在近尾流区SR模型优势巨大因其能捕捉双峰结构。远尾流区 (x/D7) RMSE0.0250.022两者性能接近因为远尾流趋于单高斯分布。计算速度极快 (微秒级)快 (毫秒级)SR模型是解析公式评估速度远超CFD略慢于但接近传统解析模型。可解释性高 (参数有物理意义)高 (公式具物理洞察)两者都可解释但SR模型可能揭示新的关系如偏移中心。外推能力中等 (依赖参数标定)待评估需要在未训练过的工况如不同湍流强度下测试。验证实操绘制速度亏损剖面对比图在几个典型的下游位置画出CFD数据点、Bastankhah模型预测线、SR模型预测线。直观展示SR模型在捕捉剖面形状上的优势。绘制尾流等值线图将整个尾流区域的ΔU预测结果以等高线形式画出与CFD结果对比。观察SR模型是否能再现尾流的整体形态和发展。误差空间分布图绘制预测误差SR预测 - CFD在y-x平面上的分布。这能清晰显示模型在哪些区域还存在系统偏差为进一步改进提供方向。6. 工程应用拓展与挑战6.1 从单机模型到风场布局优化我们建立的ΔU f(x, y)模型是针对单台风机、特定工况的。要用于实际风场需要解决两个问题多台风机的尾流叠加风场中下游风机会受到上游多台风机尾流的共同影响。最常用的方法是线性叠加ΔU_total sqrt( sum( ΔU_i^2 ) )。但线性叠加在尾流紧密交互时误差较大。我们的SR模型能否发展出更精确的叠加法则或许可以将上游多台风机的位置和强度作为输入让符号回归去学习一个“超级叠加公式”。这是一个更前沿的探索方向。工况泛化我们需要模型能响应不同的来流风速通过推力系数Ct体现、偏航角γ、湍流强度TI。这需要我们在训练数据中涵盖这些工况并将Ct,γ,TI作为额外的输入特征。公式会变成ΔU f(x, y, Ct, γ, TI)搜索空间呈指数增长对数据和算法都是巨大挑战。一个可行的策略是分步建模先建立核心的f(x, y)然后让公式中的参数如振幅、宽度成为Ct, γ, TI的函数再用符号回归去发现这些函数关系。6.2 常见问题与排查技巧在项目实践中你可能会遇到以下问题问题现象可能原因排查与解决思路公式过于复杂包含大量嵌套和无关项。1. 进化迭代次数太多过拟合。2. 复杂度惩罚权重太低。3. 搜索空间运算符集太大。1. 检查验证集损失早停。2. 增加complexity_of_operators中exp、/等的权重。3. 减少不必要的运算符如sin,cos或降低maxsize。公式始终很简单但拟合误差居高不下。1. 先验引导过强限制了搜索。2. 种群多样性不足陷入局部最优。3. 数据本身噪声大或存在系统性偏差。1. 减弱或改变先验引导方式如改用单高斯先验启动。2. 增大population_size增加ncycles_per_iteration。3. 重新检查CFD数据质量进行数据清洗。发现的公式物理意义不明确例如出现sin(exp(y))这类项。1. 数据中可能存在未被捕捉的周期性噪声2. 算法为了降低一点点误差而引入了不稳定的复杂项。1. 从帕累托前沿中选择一个稍简单、损失略高但更“顺眼”的公式。2. 在损失函数中增加对公式“平滑性”或“可微性”的惩罚。模型在训练集上好在测试集上差。过拟合。可能公式捕捉了数据中的噪声或特定工况的偶然特征。1. 使用更严格的复杂度控制。2. 增加训练数据的多样性和数据量。3. 采用交叉验证来选择模型。计算时间太长。种群规模大、迭代次数多、公式评估开销大。1. 使用PySR的并行计算功能。2. 先在数据子集上调试参数再在全集上运行。3. 考虑使用更高效的符号回归实现如OperonC库。6.3 个人实操体会与建议走完这个项目的全流程我的最深体会是符号回归不是一个“一键出结果”的魔法黑箱而是一个需要与领域知识深度对话的“共同探索”过程。先验知识的价值远超想象最初我们尝试了无引导的符号回归结果要么是精度不佳的简单公式要么是精度高但无法理解的“数学怪物”。直到我们引入了“双高斯”这个强有力的物理先验整个搜索才走上了正轨。先验知识极大地压缩了搜索空间提高了效率并保证了结果的可解释性。数据质量决定天花板如果你的CFD模拟本身就有误差或者工况覆盖不全符号回归再强大也无力回天。在数据准备阶段多花一倍时间可能在模型训练阶段节省十倍时间并获得更可靠的模型。解读公式需要物理直觉当算法吐出一个公式时真正的挑战才刚刚开始。你需要像一个侦探一样审视每一项、每一个参数思考它们可能对应的物理机制。这个过程常常能带来意想不到的物理发现这也是本项目最大的魅力所在。从“发现”到“应用”仍有距离将一个在理想CFD数据上发现的公式应用到真实的、充满不确定性的风场中需要大量的校准和验证工作。可以考虑将符号回归模型作为“基础模型”其参数再通过风场SCADA数据进行微调形成“物理AI混合模型”。这个项目为我们打开了一扇门用数据驱动的方法去发现工程物理中的隐含规律。它不仅仅适用于风机尾流对于任何有数据、有物理背景但难以用简单解析式描述的复杂现象如湍流燃烧、地质渗透、材料疲劳等符号回归加物理先验的思路都提供了一个极具潜力的研究范式。它要求我们既是懂数据的科学家也是懂物理的工程师在两者交汇处寻找新的答案。
