Forward and Backward Vector Preisach Model

Download this Jupyter File

Forward and Inverse Vector Preisach Model¶

In [1]:

import cempy as cp
from myPackage import plotEVF, getLines, plotSPTriag
import matplotlib.pyplot as plt
import numpy as np

%matplotlib widget

Scalar Preisch Model¶

With and Everett function E for a maximal inputValue $H_{max}$ = 1640 A/m and a maximal output value of $B_{max}$ = 1.5 T; distribution of the alpha and beta values are not uniform.

In [2]:

ev = cp.preisachScalar.Everett_Lorentzian(401, 1640, 1.5)
plotEVF(ev, only3D=True)

\begin{align} H_{max, in}(t) = \max\limits_{t'\in [0, t]} |H(t')|, \end{align}

and

\begin{align} B(t) = -\frac{1}{2}E(H_{max, in}, -H_{max, in}) + \sum_{k= 0}^{N-1}\Big[E(M_k, m_{k-1}) - E(M_k, m_k)\Big] \end{align}

In [3]:

p = cp.preisachScalar.preisach(ev)


Hin = getLines([0, 20, -20, 100, -100, 1000, -1000, 100, -100, 0], N = 1000)
Hin = getLines([0, 100, -100, 100, -500, 500, -500, 1000, -1000, 100, -100], N = 20)
print("len of input", len(Hin))

fig = plt.figure(1, figsize=[8, 4])
ax = fig.subplots(1, 2)

Bout = np.zeros(np.shape(Hin))

triag_plot = plotSPTriag(p, ax=ax[0], fig=fig, show=True)
line_bh, = ax[1].plot([])
ax[1].set_xlim([-1700, 1700])
ax[1].set_ylim([-1.6, 1.6])


def update(i):
    print("i = ", i, end="\r")
    Bout[i] = p.addH(Hin[i])

    if i % 4 == 0:
        line_bh.set_data(Hin[:i+1], Bout[:i+1])
        triag_plot.Update()      


import matplotlib.animation as animation
# open a new window



if False:
    ani.save("test.gif"   )

if False:
    ani = animation.FuncAnimation(fig, update, frames=len(Hin))
    from IPython.display import HTML
    HTML(ani.to_jshtml())

len of input 200

Smooth Exceeding of the Modellimits¶

with $\partial_H B$ = const. @ $\pm H_{max}$

In [4]:

Hin = getLines([0, ev.maxH*2, -ev.maxH*2, ev.maxH*2], N = 2000)


fig = plt.figure(2, figsize=[8, 4])
ax = fig.subplots(1, 2)

ax[0].set_ylabel("B in T")
ax[0].set_xlabel("H in A/m")

ax[1].set_ylabel("B in T")
ax[1].set_xlabel("H in A/m")

ax[1].set_xlim([ev.maxH*0.9, ev.maxH*1.1])
ax[1].set_ylim([ev.maxB*0.99, ev.maxB*1.01])

for val in [1, 1.05, 1.1, 1.2]:
    p.demagnetise()
    p.maxValueForOutside = val * ev.maxH
    Bout = [p.addH(H) for H in Hin]
    ax[0].plot(Hin, Bout)
    ax[1].plot(Hin, Bout)

plt.show(block=False)

Precomputing Gradient Matrix¶

With \begin{align} B(t) = -\frac{1}{2}E(H_{max, in}, -H_{max, in}) + \sum_{k= 0}^{N-1}\Big[E(M_k, m_{k-1}) - E(M_k, m_k)\Big] \end{align} and a differential permeability $\mu^\partial=\frac{\partial B}{\partial H}$. \begin{align} B(t) = -\frac{1}{2}E(H_{max, in}, -H_{max, in}) + \sum_{k= 0}^{N-1}\Big[E(M_k, m_{k-1}) - E(M_k, m_k)\Big] \end{align} therefore \begin{align} \mu^\partial = \begin{cases} \nabla E (M_{N-1}, m_{N-2}) \cdot \mathbf{e}_\alpha \quad\text{for }\partial_t H \ge 0\\ \nabla E (M_{N-2}, m_{N-1}) \cdot \mathbf{e}_\beta \quad\text{for }\partial_t H < 0 \end{cases}, \end{align} where \begin{align} \nabla E(\alpha, \beta) = \frac{E\Big(\alpha, \beta\Big) - E\Big(\alpha-h(\alpha), \beta\Big)}{h(\alpha)}\mathbf{e}_\alpha + \frac{E\Big(\alpha, \beta\Big) - E\Big(\alpha, \beta+h(\beta)\Big)}{h(\beta)}\mathbf{e}_\beta \end{align}

In [5]:

plotEVF(ev.GetDERising(), only3D=True)
plotEVF(ev.GetDEFalling(), only3D=True)

In [6]:

import time

Hin = getLines([0, 500, -500, 700, -700, 1000, -1000, 300, -300, 200, -200, 600, -600, 1100, -1100], 300)
print(len(Hin))
p.demagnetise()
dmuout_conventional = np.zeros(np.shape(Hin))
dmuout_conventional_pilot = np.zeros(np.shape(Hin))
dmuout_precomp = np.zeros(np.shape(Hin))
dmuout_precomp_pilot = np.zeros(np.shape(Hin))

p.demagnetise()
p.conventionalMatDiff=True

time_conv = time.time()
for i in range(len(Hin)):
    dmuout_conventional_pilot[i] = p.pilotH(Hin[i]).fourth
    p.addH(Hin[i])
    dmuout_conventional[i] = p.getMuDiff()
time_conv = time.time() - time_conv

p.demagnetise()
p.conventionalMatDiff=False
time_precompt = time.time()
for i in range(len(Hin)):
    dmuout_precomp_pilot[i] = p.pilotH(Hin[i]).fourth
    p.addH(Hin[i])
    dmuout_precomp[i] = p.getMuDiff()
time_precompt = time.time() - time_precompt

print("time conventional", time_conv)
print("time precomputed ", time_precompt)

plt.figure(3)
plt.clf()
plt.plot(dmuout_conventional)
plt.plot(dmuout_conventional_pilot)
plt.plot(dmuout_precomp)
plt.plot(dmuout_precomp_pilot)
plt.xlabel("time")
plt.ylabel("mu diff")
plt.show(block=False)

4200
time conventional 0.012959718704223633
time precomputed  0.012608051300048828

Inverse Model¶

add B instead of H

straigt forward: find H such that B is correct (inverse mode)

In [7]:

Bin = getLines([0, 0.1, -0.1, 0.3, -0.3, 1, -1, 0.3, -0.3], 1000)
p.demagnetise()


time_calc_inverse = time.time()
Hout = [p.addB(B) for B in Bin]
time_calc_inverse = time.time() - time_calc_inverse

print("time calc inverse = ", time_calc_inverse)

time calc inverse =  0.16762471199035645

new approach is to invert the Everett Function in a preprocessing step

In [8]:

iev = ev.copy(False)
iev = ev.Load("./save/EverettMatrix_inverseLorentzian_401_1D.bin")