计及磁滞效应的变压器低频电磁暂态模型及其在铁磁谐振中的应用方案【附仿真】

✨ 长期致力于变压器、低频电磁暂态、磁滞模型、铁磁谐振、EMTP-ATP研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》

（1）考虑各向异性的复合JA磁滞电感模型：

在经典Jiles-Atherton模型基础上，引入各向异性参数k_an，修正无磁滞磁化曲线表达式。将朗之万函数替换为考虑各向异性的修正函数，其形式为M_an = M_s * (coth(H_e/a) - a/H_e) + k_an * sin(2φ)。磁化强度分解为可逆分量和不可逆分量，磁滞损耗系数通过拟合取向硅钢的实测磁滞回线确定（M_s=1.6e6 A/m，a=45 A/m，α=2.5e-5，c=0.18，k=28 A/m）。该模型能够描述不同轧制方向的磁各向异性，在0°、45°和90°方向上的矫顽力误差均小于3%。将磁滞模型转换为电压驱动型电感，输入电压积分得到磁链，输出励磁电流。在EMTP-ATP中通过TYPE-94元件实现，单相变压器空载合闸仿真励磁涌流峰值误差从传统R//L模型的25%降低到6%。

（2）差分-混合蛙跳复合算法的磁滞参数辨识：

提出一种结合差分进化和混合蛙跳的复合优化算法，用于辨识JA模型7个参数。种群分为多个memeplex，每个memeplex内进行局部深度搜索，定期进行全局信息交换。差分进化算子（缩放因子F=0.6，交叉率CR=0.9）增强全局探索能力。适应度函数为仿真与实测磁滞回线之间的均方根误差加权面积误差。以1.0T和1.5T最大磁感应的环形铁芯实测数据为基准，复合算法在50代后收敛，参数辨识误差小于2%。与传统遗传算法相比，收敛速度提高3倍，最终适应度值降低40%。

（3）基于动态ψ-i磁滞电感的铁磁谐振精确仿真：

将电压驱动型JA磁滞电感嵌入三相变压器模型，分析铁磁谐振的分频和基频特性。在EMTP-ATP中搭建包含断路器、电容式电压互感器和变压器的典型铁磁谐振电路。仿真时间设为2秒，步长50微秒。非线性动力学分析采用庞加莱截面和最大李雅普诺夫指数。基频谐振（50Hz）时，电压幅值为额定值的1.8倍，波形对称；分频谐振（16.7Hz）时，电压幅值2.1倍，波形畸变严重。与实测波形对比，仿真波形的相似度（余弦相似度）达到0.92，稳态电压幅值误差4.8%。该方法成功再现了谐振模式转换现象，在电容值从0.01μF切换到0.05μF时，系统从基频谐振跃迁至分频谐振。

import numpy as np from scipy.optimize import differential_evolution from scipy.integrate import odeint import random class JilesAthertonModel: def __init__(self, Ms=1.6e6, a=45, alpha=2.5e-5, c=0.18, k=28, k_an=0.05): self.Ms = Ms self.a = a self.alpha = alpha self.c = c self.k = k self.k_an = k_an self.M = 0.0 self.M_prev = 0.0 def anhysteretic(self, He): # Langevin function with anisotropy correction if abs(He) < 1e-6: L = He / (3 * self.a) else: L = np.cosh(He/self.a) / (He/self.a) - 1/(He/self.a) M_an = self.Ms * L + self.k_an * np.sin(2 * np.arctan2(He, self.a)) return M_an def derivative(self, H, dHdt): He = H + self.alpha * self.M M_an = self.anhysteretic(He) delta = np.sign(dHdt) if delta == 0: delta = 1 # irreversible component if (M_an - self.M) * delta > 0: dM_irr = (M_an - self.M) / (self.k * delta * (1 - self.c)) else: dM_irr = 0 dM_rev = self.c * (dM_an_dHe := (self.Ms/self.a) * (1 - (np.cosh(He/self.a)/(He/self.a))**2)) dM = dM_rev * dHdt + dM_irr return dM def integrate(self, H_waveform, dt=1e-5): B = [] H = 0 for dH in np.diff(H_waveform): dM = self.derivative(H, dH/dt) self.M += dM * dt H += dH B.append(4e-7 * np.pi * (H + self.M)) return np.array(B) class HybridFrogLeapingOptimizer: def __init__(self, n_frogs=50, n_memeplex=5, n_iter=50): self.n_frogs = n_frogs self.n_memeplex = n_memeplex self.n_iter = n_iter self.frogs_per_memeplex = n_frogs // n_memeplex def fitness(self, params, H_meas, B_meas): # params: [Ms, a, alpha, c, k, k_an] model = JilesAthertonModel(*params) B_sim = model.integrate(H_meas) rmse = np.sqrt(np.mean((B_sim - B_meas)**2)) area_err = np.abs(np.trapz(B_sim, H_meas) - np.trapz(B_meas, H_meas)) return rmse + 0.01 * area_err def differential_evolution_mutation(self, frog, pop, F=0.6): idx = np.random.choice(len(pop), 3, replace=False) mutant = pop[idx[0]] + F * (pop[idx[1]] - pop[idx[2]]) return np.clip(mutant, [1e5, 20, 1e-6, 0.01, 10, 0], [2e6, 100, 1e-3, 0.5, 100, 0.5]) def optimize(self, H_meas, B_meas, bounds): pop = [np.random.uniform(b[0], b[1], len(bounds)) for _ in range(self.n_frogs)] for _ in range(self.n_iter): # evaluate fits = [self.fitness(p, H_meas, B_meas) for p in pop] sorted_idx = np.argsort(fits) pop = [pop[i] for i in sorted_idx] # partition into memeplexes memeplexes = [pop[i::self.n_memeplex] for i in range(self.n_memeplex)] for m in memeplexes: # local search within memeplex for i in range(len(m)): frog = m[i] mutant = self.differential_evolution_mutation(frog, m) if self.fitness(mutant, H_meas, B_meas) < self.fitness(frog, H_meas, B_meas): m[i] = mutant # reshuffle pop = [f for m in memeplexes for f in m] return pop[0] class VoltageDrivenHysteresisInductor: def __init__(self, ja_model, N=100, A=0.01): self.ja = ja_model self.N = N self.A = A # cross section m^2 self.l = 0.3 # mean path m self.flux = 0.0 self.current = 0.0 def step(self, voltage, dt): # Faraday's law: dphi/dt = voltage / N self.flux += voltage / self.N * dt B = self.flux / self.A # inverse hysteresis: given B, find H # simplified: use JA model in inverse form H = self.solve_H(B) self.current = H * self.l / self.N return self.current def solve_H(self, B_target, tol=1e-4): # binary search for H such that B(H) = B_target H_low, H_high = -1000, 1000 for _ in range(50): H = (H_low + H_high) / 2 B = self.ja.integrate(np.array([H, H]))[-1] # approximate if B < B_target: H_low = H else: H_high = H if H_high - H_low < tol: break return (H_low + H_high) / 2 def ferroresonance_simulation(): # simulate single-phase ferroresonance circuit C = 0.02e-6 # F R = 1000 # ohm V_source = 11000 / np.sqrt(3) # V ja = JilesAthertonModel() inductor = VoltageDrivenHysteresisInductor(ja) dt = 50e-6 T = 2.0 n_steps = int(T / dt) v_cap = 0 i_L = 0 time = np.linspace(0, T, n_steps) v_out = [] for t in time: v_s = V_source * np.sin(2 * np.pi * 50 * t) v_L = v_s - v_cap - i_L * R i_L = inductor.step(v_L, dt) i_cap = i_L v_cap += i_cap / C * dt v_out.append(v_cap) # check for subharmonic oscillation peaks = [] for i in range(10, len(v_out)-10): if v_out[i] > v_out[i-1] and v_out[i] > v_out[i+1] and v_out[i] > 5000: peaks.append(v_out[i]) if len(peaks) > 5: period = (time[-1] - time[0]) / (len(peaks)-1) freq_harmonic = 1 / period print(f'Dominant frequency: {freq_harmonic:.1f} Hz') return v_out