CN109086561B - 一种计及各向异性的ja磁滞模型参数提取方法 - Google Patents

Abstract

本发明公开了一种计及各向异性的JA磁滞模型参数提取方法，该方法包括以下步骤：S1.建立考虑各向异性的无磁滞磁化曲线模型；S2.求解JA磁滞电感；S3.采用差分‑混合蛙跳复合算法获取JA磁滞模型参数。本发明提出的计及各向异性的JA磁滞模型参数提取方法，从原理上解决了JA磁滞模型参数传统提取方法所存在的技术经济难题，且计算速度快、精度高，具有十分重要的应用价值。

Description

一种计及各向异性的JA磁滞模型参数提取方法

技术领域

本发明涉及电力技术领域，特别是一种计及各向异性的JA(Jiles-Atherton)磁滞模型参数提取方法。

背景技术

变压器饱和特性和铁芯损耗在低频电磁暂态过程中具有十分重要的作用。在磁化过程中存在磁感应强度变化滞后于磁场强度变化的现象，即磁滞效应。现有研究通常采用单值非线性电感和电阻的并联组合表征磁滞特性，不能准确模拟光滑的磁滞回线；现有JA磁滞模型采用朗之万函数描述无磁滞磁化曲线，其结构简单、计算方便，但在模拟各向异性特性方面还存在明显不足。此外，对经典JA磁滞模型进行修正后，模型复杂程度显著增加，而在JA模型参数的优化算法方面，现有研究表明GA和SA在对经典JA磁滞模型参数辨识时能够取得一定的求解精度但在求解速度和收敛性能方面还存在不足。因此，开展计及各向异性的JA磁滞模型参数提取方法研究对于保障变压器电磁暂态模型的应用具有极其重要的理论意义。

发明内容

有鉴于此，为了解决上述问题，本发明提供一种计及各向异性的JA磁滞模型参数提取方法，为变压器电磁暂态模型研究提供了基础。

为实现上述目的及其他目的，本发明提供一种计及各向异性的JA磁滞模型参数提取方法，该方法包括以下步骤：

S1.建立考虑各向异性的无磁滞磁化曲线模型；

其中，M_s为饱和磁化强度，E_m(θ)表示磁矩能量，k_B表示Boltzman常数，θ表示原子磁矩m_at和磁场方向之间的夹角；

S2.求解JA磁滞电感；

S3.采用差分-混合蛙跳复合算法获取JA磁滞模型参数。

可选地，所述求解JA磁滞电感，具体包括以下步骤：

①当t＝(t₀+Δt)时，该时刻电感两端电压u(t₀+Δt)已知，由式dB/dt＝u(t₀+Δt)/(N·S)对电压进行积分变换，依据B(t₀+Δt)＝B(t₀)+ΔB/Δt·Δt≈B(t₀)+dB/dt·Δt得到磁感应强度B(t₀+Δt)；

N为绕组匝数，S为铁芯横截面积，Δt为时间增量，ΔB为磁感应强度增量；

②将磁感应强度B(t₀)，磁场强度H(t₀)，磁化强度M(t₀)，dB/dt，以及dH/dt作为Matlab Function模块输入，计算得到dM/dt；

③利用表达式M(t₀+Δt)＝M(t₀)+dM/dt·Δt，计算(t₀+Δt)时刻磁化强度M(t₀+Δt)，由式H(t₀+Δt)＝B(t₀+Δt)/μ₀–M(t₀+Δt)可以计算得到磁场强度H(t₀+Δt)，根据安培环路定律N·i＝H·l，可以输出此时电流i(t₀+Δt)。

可选地，所述采用差分-混合蛙跳复合算法获取JA磁滞模型参数，具体包括以下步骤：

S31.获取试验磁场强度数据H_T；

S32.随机生成参数饱和磁化强度M_s，无磁滞磁化曲线形状参数a，平均场参数α，磁场方向和晶体特征轴方向夹角ψ，牵制系数k，磁化因数c，单轴各向异性的平均能量密度K_an；

S33.求解JA微分方程；

S34.根据试验磁感应强度B_T，计算每一青蛙适应度值E，并据此排序；

S35.划分m_p个青蛙子组族群，并在每一个子族群内进行局部深度搜索；

S36.混合整个族群并重新排序；

S37.生成新一代青蛙种群；

S38.判断是否达到最大迭代次数，如果是则输出最优青蛙个数；如果否，则重复步骤S35～步骤S37。

可选地，所述在每一个子族群内进行局部深度搜索，具体包括以下步骤；

S351.记录本子族群的局部最优青蛙U_b和局部最差青蛙U_w以及整个种群的全局最优青蛙U_g；

S352.进行最差青蛙第一次更新；

S353.判断更新前后青蛙性能是否提高；若是，则取代本子族群局部最差青蛙U_w；若否，则用全局最优青蛙U_g替代局部最优青蛙U_b进行第二次更新；

S354.再次比较更新前后青蛙性能是否提高，如果是，则取代本子族群局部最差青蛙U_w，如果否，则随机产生一个青蛙取代本子族群局部最差青蛙U_w；

S355.判断是否达到局部搜索次数，若是则局部搜索结速，若否则重复步骤S351～S354。

由于采用了上述技术方案，本发明具有如下的优点：

本发明提出的计及各向异性的JA磁滞模型参数提取方法，从原理上解决了JA磁滞模型参数传统提取方法所存在的技术经济难题，且计算速度快、精度高，具有十分重要的应用价值。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述：

图1为本发明所述方法的流程图；

图2为JA磁滞电感在MATLAB中求解；

图3为JA参数DE-SFLA优化流程图；

图4为取向和无取向硅钢磁滞回线测量的原理图；

图5为GO Bmax＝1.0T的试验与优化后仿真磁滞回线；

图6为GO Bmax＝1.5T的试验与优化后仿真磁滞回线；

图7为NGO Bmax＝1.0T的试验与优化后仿真磁滞回线；

图8为NGO Bmax＝1.5T的试验与优化后仿真磁滞回线。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。

需要说明的是，本实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。

如图1所示，本发明提供一种计及各向异性的JA磁滞模型参数提取方法，包括：

S1：建立考虑各向异性的无磁滞磁化曲线模型。

S2：求解JA磁滞电感。

S3：采用差分-混合蛙跳复合算法获取本文JA磁滞模型参数，并开展取向硅钢和无取向硅钢磁滞回线测量，分别比较不同算法在优化速度、精度以及收敛速度等方面的性能差异。

下面以一具体实施方式对本发明的计及各向异性的JA磁滞模型参数提取方法进行说明。具体步骤为：

S1：建立考虑各向异性的无磁滞磁化曲线模型。

当施加外磁场后，磁性材料磁化强度可能从易磁化轴方向转出，表现出磁晶各向异性，晶体磁晶各向异性能E_k可以用自发磁化与相互正交的晶体学主轴的方向余弦函数表示：

E_K＝K₁(α₁ ²α₂ ²+α₂ ²α₃ ²+α₃ ²α₁ ²)+K₂α₁ ²α₂ ²α₃ ² (1)

式中：K₁，K₂为立方晶体磁晶各向异性常数，α₁，α₂和α₃为方向余弦。

根据式(1)建立一种考虑各向异性的无磁滞磁化曲线模型，对于各向异性材料，由于各向异性能在不同方向不一致，此时无磁滞磁化曲线只能用磁矩在外磁场方向的总和进行计算，可采用Boltzman分布表示为：

式中：M_s为饱和磁化强度，θ表示原子磁矩m_at和磁场方向之间的夹角，k_B为Boltzman常数。磁矩能量E_m(θ)可以表示为：

E_m(θ)＝-μ₀·m_at·H·cosθ (3)

式中：H表示磁场强度；μ₀为真空磁导率。

对于各向异性材料，磁矩与磁场方向夹角不尽相同。在二维平面上，设磁矩m₁，m₂与磁场夹角为θ，简化后磁晶各向异性能为：

式中：K₀和K₁为各向异性常数，可通过试验测量得到，ψ为磁场方向和晶体特征轴方向夹角。

S2：求解JA磁滞电感。

比较各向异性函数与传统朗之万函数的性能差异，分析ode45，ode23以及四阶定步长Runge-Kutta算法对复合JA磁滞电感求解时间的区别；

在Matlab中求解上述JA磁滞电感的过程。如图2所示，具体求解步骤如下：

①当t＝(t₀+Δt)时，该时刻电感两端电压u(t₀+Δt)已知，由式dB/dt＝u(t₀+Δt)/(N·S)对电压进行积分变换，依据B(t₀+Δt)＝B(t₀)+ΔB/Δt·Δt≈B(t₀)+dB/dt·Δt得到磁感应强度B(t₀+Δt)；

N为绕组匝数，S为铁芯横截面积，Δt为时间增量，ΔB为磁感应强度增量；

②将磁感应强度B(t₀)，磁场强度H(t₀)，磁化强度M(t₀)，dB/dt，以及dH/dt作为Matlab Function模块输入，计算得到dM/dt；

③利用表达式M(t₀+Δt)＝M(t₀)+dM/dt·Δt，计算(t₀+Δt)时刻磁化强度M(t₀+Δt)，由式H(t₀+Δt)＝B(t₀+Δt)/μ₀–M(t₀+Δt)可以计算得到磁场强度H(t₀+Δt)，根据安培环路定律N·i＝H·l，可以输出此时电流i(t₀+Δt)。

以典型的单轴晶体各向异性材料模型参数为例，即k₁＝0.011，k₂＝0.049，a＝1000A/m，k＝5000A/m，c＝0.1，M_s＝1.3×10⁶A/m，α＝0.001，K_an＝4×10⁴J/m³，ψ＝45^o为例，分别采用ode45，ode23，以及四阶定步长Runge-Kutta法通过Matlab编程求解JA磁滞电感，，计算机操作系统为：Windows 7Enterprise,系统类型：64位，CPU:Intel Core i5–46703.40GHz。不同求解算法下单次JA磁滞模型求解所需时间如表1所示，四阶定步长Runge-Kutta算法求解效率最高，其求解时间仅为ode23算法的26.8％，ode45求解时间介于四阶定步长Runge-Kutta算法和ode23算法之间。因此，在对JA磁滞模型进行求解时采用四阶定步长Runge-Kutta算法最合适。

表1求解算法对本文JA磁滞模型计算时间影响

S3：采用差分-混合蛙跳复合算法获取本文JA磁滞模型参数，并开展取向硅钢和无取向硅钢磁滞回线测量，分别比较不同算法在优化速度、精度以及收敛速度等方面的性能差异：

引入DE/target-to-best/1型差分变异算子后的青蛙更新步长为：

S₁＝rand·(U_b-U_w)+rand·(U_r1-U_r2) (5)

式中：U_r1和U_r2为本族群中与最差和最优青蛙均不同的两随机个体。

提出如下自适应步长系数：

式中：M₁和M₂分别为当前迭代次数和最大迭代次数。

引入自适应步长系数后的青蛙更新步长变为：

S₂＝k₁·(U_b-U_w)+k₂·(U_r1-U_r2) (7)

引入惯性权值α，可提高原始算法的局部搜索和全局搜索能力，从而得到新的更新公式：

U_w ^new＝α·U_w+S₂ (8)

式中：α为惯性权值，其取值范围为[0,1]。

依据上述原理，对经典SFLA中青蛙更新策略分别引入差分变异思想、自适应步长系数以及惯性权值后即可得到DE-SFLA的优化原理。采用DE-SFLA对JA模型参数优化的流程如图3所示。

S31.获取试验磁场强度数据H_T；

S32.随机生成参数饱和磁化强度M_s，无磁滞磁化曲线形状参数a，平均场参数α，磁场方向和晶体特征轴方向夹角ψ，牵制系数k，磁化因数c，单轴各向异性的平均能量密度K_an；

S33.求解JA微分方程；

式中：δ表示方向系数，M_an表示无磁滞磁化强度。

S34.根据试验磁感应强度B_T，计算每一青蛙适应度值E，并据此排序；

式中：N为试验数据长度，B_Tj为试验磁感应强度，B_Sj为计算所得磁感应强度，B_Tmax为试验磁感应强度最大值。

S35.划分m_p个青蛙子组族群，并在每一个子族群内进行局部深度搜索；

S36.混合整个族群并重新排序；

S37.生成新一代青蛙种群；

S38.判断是否达到最大迭代次数，如果是则输出最优青蛙个数；如果否，则重复步骤S35～步骤S37。

可选地，所述在每一个子族群内进行局部深度搜索，具体包括以下步骤；

S351.记录本子族群的局部最优青蛙U_b和局部最差青蛙U_w以及整个种群的全局最优青蛙U_g；

S352.按式(11)进行最差青蛙第一次更新；

α＝(α₂-α₁)·(N₂-N₁)/N₂+α₁ (11)

式中：α₁和α₂分别为惯性权值最小值和最大值；N₁和N₂分别为当前局部搜索迭代次数和局部搜索最大迭代次数。

S353.判断更新前后青蛙性能是否提高；若是，则取代本子族群局部最差青蛙U_w；若否，则用全局最优青蛙U_g替代式(7)局部最优青蛙U_b进行第二次更新；

S354.再次比较更新前后青蛙性能是否提高，如果是，则取代本子族群局部最差青蛙U_w，如果否，则随机产生一个青蛙取代本子族群局部最差青蛙U_w；

S355.判断是否达到局部搜索次数，若是则局部搜索结速，若否则重复步骤S351～S354。

本实施例采用湖南联众科技有限公司的MATS–2010M硅钢测量装置对取向(GO)和无取向(NGO)硅钢环形铁芯(20VA，10V)进行测量，试验原理如图4所示。1为电流探头，2为绕组，3为样品，4为计算机。相同尺寸的GO和NGO铁芯几何参数如表2所示。采用MATS–2010M硅钢测量装置固定最大磁感应强度B_m模式，即设定B_m值，励磁电源自动加压直到B＝B_m，在频率为50Hz下，分别对B_m＝1.0T和B_m＝1.5T的磁滞回线进行测量。

表2铁芯几何尺寸

分别采用上节所述优化算法，在Matlab 2013a中依据实测B-H磁滞回线对JA模型参数进行辨识(计算机操作系统：64位Windows 7 Enterprise,CPU:Intel Core i5–46703.40GHz)。

优化后的磁滞回线结果如图5～8所示，从图中可以看出，DE、GA和DE-SFLA三种算法在GO和NGO两种材料不同磁感应强度下均能达到全局最优解附近，而SA在磁感应强度未饱和时(B_max＝1.0T)，其结果能够收敛于全局最优解，而当磁感应强度饱和时(B_max＝1.5T)，其优化结果为局部最优。在对SA的温度函数以及参数邻域矩阵进行反复调整后，其结果依旧陷入局部最优。因此，SA在磁滞回线饱和情况下对JA模型参数进行优化时容易得到局部最优解。

本发明提出的计及各向异性的JA磁滞模型参数提取方法，从原理上解决了JA磁滞模型参数传统提取方法所存在的技术经济难题，且计算速度快、精度高，具有十分重要的应用价值。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的保护范围当中。

Claims (3)

1.一种计及各向异性的JA磁滞模型参数提取方法，其特征在于，该方法包括以下步骤：

S1.建立考虑各向异性的无磁滞磁化曲线模型；

其中，M_s为饱和磁化强度，E_m(θ)表示磁矩能量，k_B表示Boltzman常数，θ表示原子磁矩m_at和磁场方向之间的夹角；

S2.求解JA磁滞电感；

S3.采用差分-混合蛙跳复合算法获取JA磁滞模型参数；

所述采用差分-混合蛙跳复合算法获取JA磁滞模型参数，具体包括以下步骤：

S31.获取试验磁场强度数据H_T；

S32.随机生成参数饱和磁化强度M_s，无磁滞磁化曲线形状参数a，平均场参数α，磁场方向和晶体特征轴方向夹角ψ，牵制系数k，磁化因数c，单轴各向异性的平均能量密度K_an；

S33.求解JA微分方程；

S34.根据试验磁感应强度B_T，计算每一青蛙适应度值E，并据此排序；

S35.划分m_p个青蛙子组族群，并在每一个子族群内进行局部深度搜索；

S36.混合整个族群并重新排序；

S37.生成新一代青蛙种群；

S38.判断是否达到最大迭代次数，如果是则输出最优青蛙个体，如果否，则重复步骤S35～步骤S37，青蛙个体即JA磁滞模型参数。

2.根据权利要求1所述的一种计及各向异性的JA磁滞模型参数提取方法，其特征在于，所述求解JA磁滞电感，具体包括以下步骤：

①当t＝(t₀+Δt)时，该时刻电感两端电压u(t₀+Δt)已知，由式dB/dt＝u(t₀+Δt)/(N·S)对电压进行积分变换，依据B(t₀+Δt)＝B(t₀)+ΔB/Δt·Δt≈B(t₀)+dB/dt·Δt得到磁感应强度B(t₀+Δt)；

N为绕组匝数，S为铁芯横截面积，Δt为时间增量，ΔB为磁感应强度增量；

②将磁感应强度B(t₀)，磁场强度H(t₀)，磁化强度M(t₀)，dB/dt，以及dH/dt作为MatlabFunction模块输入，计算得到dM/dt；

③利用表达式M(t₀+Δt)＝M(t₀)+dM/dt·Δt，计算(t₀+Δt)时刻磁化强度M(t₀+Δt)，由式H(t₀+Δt)＝B(t₀+Δt)/μ₀–M(t₀+Δt)可以计算得到磁场强度H(t₀+Δt)，根据安培环路定律N·i＝H·l，可以输出此时电流i(t₀+Δt)。

3.根据权利要求1所述的一种计及各向异性的JA磁滞模型参数提取方法，其特征在于，所述在每一个子族群内进行局部深度搜索，具体包括以下步骤；

S351.记录本子族群的局部最优青蛙U_b和局部最差青蛙U_w以及整个种群的全局最优青蛙U_g；

S352.进行最差青蛙第一次更新；

S353.判断更新前后青蛙性能是否提高；若是，则取代本子族群局部最差青蛙U_w；若否，则用全局最优青蛙U_g替代局部最优青蛙U_b进行第二次更新；

S354.再次比较更新前后青蛙性能是否提高，如果是，则取代本子族群局部最差青蛙U_w，如果否，则随机产生一个青蛙取代本子族群局部最差青蛙U_w；

S355.判断是否达到局部搜索次数，若是则局部搜索结速，若否则重复步骤S351～S354。