Loss models for efficiency mapping

[KizzieNova]

Hi

anyone that got the loss models working on an SPM motor, in order to realize efficiency maps?

best regards
Kizzie

[SebGue]

Hi @KizzieNova ,

if you could post your actual simulation script, I will try to have a look.

Best regards,
Sebastian

[KizzieNova]

Hi Sebastian.

Below my code. I tried first with only copper losses, but I'm strugging with the structure of Pyleecan regarding the loss models.
Also iron losses don't work atm.

import numpy as np
import matplotlib.pyplot as plt

from pyleecan.Functions.load import load
from pyleecan.Classes.Simu1 import Simu1
from pyleecan.Classes.InputCurrent import InputCurrent
from pyleecan.Classes.OPdq import OPdq
from pyleecan.Classes.VarLoadCurrent import VarLoadCurrent
from pyleecan.Classes.DataKeeper import DataKeeper
from pyleecan.Classes.EEC_PMSM import EEC_PMSM
from pyleecan.Classes.Electrical import Electrical


from pyleecan.Classes.Loss import Loss
from pyleecan.Classes.LossModelJoule import LossModelJoule


machine = load("motor.json")

simu = Simu1(name="EfficiencyMap", machine=machine)
simu.input = InputCurrent()
simu.input.OP = OPdq(N0=1000, Id_ref=0, Iq_ref=0)

simu.elec = Electrical()
simu.elec.eec = EEC_PMSM()


loss = Loss()
loss.model_dict = {"Pcu": LossModelJoule(name="Pcu")}
print(loss)  # Dit zou nu moeten werken zonder error

simu.loss = loss


# --- 3. Rooster voor snelheid en Iq ---
Nspeed = 10
Niq = 10

speed_range = np.linspace(1000, 10000, Nspeed)
Iq_range = np.linspace(0, 50, Niq)

Speed_grid, Iq_grid = np.meshgrid(speed_range, Iq_range)

OP_matrix = np.column_stack((
    Speed_grid.flatten(),                # N0
    np.zeros_like(Speed_grid.flatten()),# Id_ref = 0
    Iq_grid.flatten()                   # Iq_ref
))

# --- 4. VarLoadCurrent instellen ---
varload = VarLoadCurrent()
varload.set_OP_array(OP_matrix, "N0", "Id_ref", "Iq_ref")
varload.is_reuse_femm_file = True

# --- 5. DataKeepers toevoegen ---
varload.datakeeper_list = [
    DataKeeper(name="Tem_av", symbol="Tem_av", unit="Nm",
               keeper=lambda out: getattr(out.elec, "Tem_av", 0)),
    DataKeeper(name="Pem_av", symbol="Pem_av", unit="W",
               keeper=lambda out: getattr(out.elec, "Pem_av", 0)),
    DataKeeper(name="Pcu", symbol="Pcu", unit="W",
               keeper=lambda out: getattr(out.elec, "losses", {}).get("Pcu", 0)),
    DataKeeper(name="I0", symbol="I0", unit="A",
               keeper=lambda out: out.elec.OP.get_I0_Phi0()["I0"]),
    DataKeeper(name="Phi0", symbol="Phi0", unit="rad",
               keeper=lambda out: out.elec.OP.get_I0_Phi0()["Phi0"])
]

# --- 6. Simulatie uitvoeren ---
simu.var_simu = varload



Xout = simu.run()
print(Xout.mag.meshsolution)
print(list(Xout.mag.meshsolution.values()))

# --- 7. Resultaten verwerken ---
Tem_av = np.array(Xout["Tem_av"].result).reshape((Niq, Nspeed))
Pem_av = np.array(Xout["Pem_av"].result).reshape((Niq, Nspeed))
Pcu = np.array(Xout["Pcu"].result).reshape((Niq, Nspeed))
I0 = np.array(Xout["I0"].result).reshape((Niq, Nspeed))

Ptot = Pcu
efficiency = np.where(Pem_av > 1e-3, Pem_av / (Pem_av + Ptot), 0)

# --- 8. Efficiëntiemap plotten ---
plt.figure(figsize=(8,6))
cp = plt.contourf(Speed_grid, Iq_grid, efficiency * 100, levels=20, cmap="viridis")
plt.colorbar(cp, label="Efficiëntie [%]")
plt.xlabel("Snelheid [rpm]")
plt.ylabel("Iq [A]")
plt.title("Efficiëntiemap (Id=0)")
plt.grid(True)
plt.tight_layout()
plt.show()

# --- 9. Pcu verliezen plotten ---
plt.figure(figsize=(8,6))
cp2 = plt.contourf(Speed_grid, Iq_grid, Pcu, levels=20, cmap="Reds")
plt.colorbar(cp2, label="Pcu [W]")
plt.xlabel("Snelheid [rpm]")
plt.ylabel("Iq [A]")
plt.title("Wikkelingsverliezen (Pcu)")
plt.grid(True)
plt.tight_layout()
plt.show()

Edit: Format code for readability.

[SebGue]

Hi @KizzieNova ,

you're right, there is some issue, i.e. I get an error in Meshsolution.get_solution with your script. But atm I can't tell if it's an issue with pyleecan or with your script. Unfortunately I need some time to have a closer look.

Nevertheless, have you already tried 'test_FEMM_Loss_SPMSM_18s16p.py'. Its a loss calculation on a single OP, but it should work. Maybe you could also adapt this to a simulation with ElecLUTdq to compute a LUT. If you save the LUT, you can load it into your electrical simulation and avoid long MagFEMM calculations every time you run your script.

Edit: You could set simu.input.Nt_tot = 16 * 4 in the test to get some quick results.

Best regards,
Sebastian

[KizzieNova]

I Made a new attempt, unfortunately still struggling with the saving of the torque and efficiency data, which are not correctly saved.

from os.path import join, split, exists

import numpy as np
import pytest

from SciDataTool.Functions.Plot.plot_3D import plot_3D

from pyleecan.Classes.ElecLUTdq import ElecLUTdq
from pyleecan.Classes.OPdq import OPdq
from pyleecan.Classes.Simu1 import Simu1
from pyleecan.Classes.InputCurrent import InputCurrent
from pyleecan.Classes.VarLoadCurrent import VarLoadCurrent
from pyleecan.Classes.MagFEMM import MagFEMM
from pyleecan.Classes.PostLUT import PostLUT
from pyleecan.Classes.Loss import Loss
from pyleecan.Classes.LossModelSteinmetz import LossModelSteinmetz
from pyleecan.Classes.LossModelJoule import LossModelJoule
from pyleecan.Classes.LossModelProximity import LossModelProximity
from pyleecan.Classes.LossModelMagnet import LossModelMagnet
from pyleecan.Classes.DataKeeper import DataKeeper

from pyleecan.Functions.load import load
from pyleecan.Functions.Load.load_json import LoadMissingFileError

@pytest.mark.long_5s
@pytest.mark.long_1m
@pytest.mark.MagFEMM
@pytest.mark.EEC_PMSM
@pytest.mark.SPM
@pytest.mark.periodicity
@pytest.mark.skip(reason="Work in progress")
def test_efficiency_map_spm():
"""Efficiency map for SPM motor (Id=0, sweep over speed and Iq)."""

DATA_DIR = r'c:\motor'
RESULT_DIR = r'c:\motor'

motor= load(join(DATA_DIR, "motor.json"))
LUT_file_name = f"LUT_eff_motor_SPM.h5"
path_to_LUT = join(RESULT_DIR, LUT_file_name)

if not exists(split(path_to_LUT)[0]):
    raise Exception("The path to LUT is not valid.")

    
Nspeed = 50
Nload = 2
speed_vect = np.linspace(500, 10000, Nspeed)
load_vect = np.linspace(0, 1, Nload)
Iq_max = 10

T = np.zeros((Nspeed, Nload))
eff = np.zeros((Nspeed, Nload))
J = np.zeros((Nspeed, Nload))

for i, N0 in enumerate(speed_vect):
    simu = Simu1(name=f"eff_map_{int(N0)}", machine=motor)
    simu.input = InputCurrent(
        OP=OPdq(N0=N0),
        Nt_tot=4 * 20,
        Na_tot=200 * 8,
        is_periodicity_a=True,
        is_periodicity_t=True,
    )
    simu.elec = ElecLUTdq(
        Urms_max=320,
        Jrms_max=27e6,
        n_interp=100,
        n_Id=1,
        n_Iq=Nload,
        Id_max=0,
        Id_min=0,
        Iq_min=0,
        Iq_max=Iq_max,
        LUT_enforced=None,
        is_grid_dq=True,
        Tsta=120,
        type_skin_effect=1,
        LUT_simu=Simu1(
            input=InputCurrent(
                OP=OPdq(),
                Nt_tot=4 * 10,
                Na_tot=200 * 8,
                is_periodicity_a=True,
                is_periodicity_t=True,
            ),
            var_simu=VarLoadCurrent(
                postproc_list=[PostLUT(is_save_LUT=True, file_name=LUT_file_name)],
                is_keep_all_output=True,
            ),
            mag=MagFEMM(
                is_periodicity_a=True,
                is_periodicity_t=True,
                nb_worker=4,
                is_get_meshsolution=True,
            ),
            loss=Loss(
                is_get_meshsolution=False,
                Tsta=100,
                model_dict={
                    "stator core": LossModelSteinmetz(group="stator core"),
                    "rotor core": LossModelSteinmetz(group="rotor core"),
                    "joule": LossModelJoule(group="stator winding"),
                    "proximity": LossModelProximity(group="stator winding"),
                    "magnets": LossModelMagnet(group="rotor magnets"),
                },
            ),
        ),
    )
    out = simu.run()
    T[i, :] = np.array(out.elec.Tem_av)
    print(T[i, :])
    eff[i, :] = np.array(out.elec.OP.efficiency)
    print(eff[i, :])
    J[i, :] = np.array(out.elec.get_Jrms() * 1e-6)
    print(J[i, :])

speed_grid = speed_vect[:, None] * np.ones((1, Nload))
Iq_grid = np.ones((Nspeed, 1)) * (load_vect * Iq_max)

plot_3D(
    Xdata=speed_grid,
    Ydata=T,
    Zdata=eff,
    xlabel="Speed [rpm]",
    ylabel="Torque [Nm]",
    zlabel="Efficiency",
    title="SPM Efficiency map",
    type_plot="pcolormesh",
    is_contour=True,
    levels=[0.7, 0.85, 0.9, 0.92, 0.93, 0.94, 0.95],
    gamma=5,
)

return out

To run it without pytest

if name == "main":
out = test_efficiency_map_spm()
input("Press Enter to continue...")


[SebGue]

Hello @KizzieNova ,

here is some modification of your code that runs on my side. Give it a try and feel free to ask if there is any question ...

from os.path import join, split, exists
import numpy as np

from SciDataTool.Functions.Plot.plot_3D import plot_3D


from pyleecan.Classes.ElecLUTdq import ElecLUTdq
from pyleecan.Classes.OPdq import OPdq
from pyleecan.Classes.Simu1 import Simu1
from pyleecan.Classes.InputCurrent import InputCurrent
from pyleecan.Classes.VarLoadCurrent import VarLoadCurrent
from pyleecan.Classes.MagFEMM import MagFEMM
from pyleecan.Classes.PostLUT import PostLUT
from pyleecan.Classes.Loss import Loss
from pyleecan.Classes.LossModelSteinmetz import LossModelSteinmetz
from pyleecan.Classes.LossModelJoule import LossModelJoule
from pyleecan.Classes.LossModelProximity import LossModelProximity
from pyleecan.Classes.LossModelMagnet import LossModelMagnet
from pyleecan.Classes.DataKeeper import DataKeeper


from pyleecan.Functions.load import load
from pyleecan.Functions.Load.load_json import LoadMissingFileError

from pyleecan.definitions import DATA_DIR, RESULT_DIR

def test_efficiency_map_spm():
    """Efficiency map for SPM motor (Id=0, sweep over speed and Iq)."""


    # DATA_DIR = r'c:\motor'
    # RESULT_DIR = r'c:\motor'

    #motor= load(join(DATA_DIR, "motor.json"))

    motor= load(join(DATA_DIR, "Machine", "SPMSM_001.json"))
    LUT_file_name = f"LUT_eff_motor_SPM.h5"
    path_to_LUT = join(RESULT_DIR, LUT_file_name)

    if not exists(split(path_to_LUT)[0]):
        raise Exception("The path to LUT is not valid.")

    Nspeed = 5
    Nload = 5
    speed_vect = np.linspace(500, 10000, Nspeed)
    load_vect = np.linspace(0, 1, Nload)
    Iq_max = 10

    T = np.zeros((Nspeed, Nload))
    eff = np.zeros((Nspeed, Nload))
    J = np.zeros((Nspeed, Nload))

    # Steinmetz loss model parameter for lamination
    f_ref = 50 
    B_ref = 1.5
    kwarg = dict(k_hy=1.5/f_ref/B_ref**2, k_ed=1/f_ref**2/B_ref**2, alpha_f=1, alpha_B=2)

    # simulation definition
    simu = Simu1(name=f"eff_map", machine=motor)
    simu.input = InputCurrent(
        OP=OPdq(N0=speed_vect[0]),
        Nt_tot=4 * 20,
        Na_tot=200 * 8,
        is_periodicity_a=True,
        is_periodicity_t=True,
    )
    simu.elec = ElecLUTdq(
        Urms_max=320,
        Jrms_max=27e6,
        n_interp=100,
        n_Id=3,
        n_Iq=3,
        Id_max=0,
        Id_min=-Iq_max,
        Iq_min=0,
        Iq_max=Iq_max,
        LUT_enforced=None,
        is_grid_dq=True,
        Tsta=120,
        type_skin_effect=1,
        LUT_simu=Simu1(
            input=InputCurrent(
                OP=OPdq(),
                Nt_tot=4 * 10,
                Na_tot=200 * 8,
                is_periodicity_a=True,
                is_periodicity_t=True,
            ),
            var_simu=VarLoadCurrent(
                postproc_list=[PostLUT(is_save_LUT=True, file_name=LUT_file_name)],
                is_keep_all_output=True,
            ),
            mag=MagFEMM(
                is_periodicity_a=True,
                is_periodicity_t=True,
                nb_worker=4,
                is_get_meshsolution=True,
            ),
            loss=Loss(
                is_get_meshsolution=False,
                Tsta=100,
                model_dict={
                    "stator core": LossModelSteinmetz(group="stator core", **kwarg),
                    "rotor core": LossModelSteinmetz(group="rotor core", **kwarg),
                    "joule": LossModelJoule(group="stator winding"),
                    "proximity": LossModelProximity(group="stator winding"),
                    "magnets": LossModelMagnet(group="rotor magnets"),
                },
            ),
        ),
    )

    # try to load the LUT
    try:
        LUT_enforced = load(path_to_LUT)
        simu.elec.LUT_enforced = LUT_enforced
    except (FileNotFoundError, LoadMissingFileError):
        print("The LUT could not be loaded, so it will be computed.")
        LUT_enforced = None

    # run electrical simulation
    for jj, load_rate in enumerate(load_vect):
        for ii, N0 in enumerate(speed_vect):
            if ii > 0 and LUT_enforced is None:
                simu.elec.LUT_enforced = load(path_to_LUT)

            simu.elec.load_rate = load_rate
            simu.input.OP.N0 = N0
            simu.elec.Irms_max = None

            out = simu.run()
            T[ii, jj] = np.array(out.elec.Tem_av)
            
            eff[ii, jj] = np.array(out.elec.OP.efficiency)
            J[ii, jj] = np.array(out.elec.get_Jrms() * 1e-6)
            print(f"load = {load_rate:.2f} - speed = {N0:.1f} rpm - {T[ii, jj]:.1f} {eff[ii, jj]:.2f} {J[ii, jj]:.1f}")

    speed_grid = speed_vect[:, None] * np.ones((1, Nload))
    Iq_grid = np.ones((Nspeed, 1)) * (load_vect * Iq_max)

    plot_3D(
        Xdata=speed_grid,
        Ydata=T,
        Zdata=eff,
        xlabel="Speed [rpm]",
        ylabel="Torque [Nm]",
        zlabel="Efficiency",
        title="SPM Efficiency map",
        type_plot="pcolormesh",
        is_contour=True,
        # levels=[0.7, 0.85, 0.9, 0.92, 0.93, 0.94, 0.95],
        # gamma=5,
        is_show_fig=True,
    )

    return out

# To run it without pytest
if __name__ == "__main__":
    out = test_efficiency_map_spm()
    input("Press Enter to continue...")

[KizzieNova]

Thanks for the helpful suggestions.

This code runs without errors, although there seems to be a bug in where the LUT results are stored, the defined path is not used, but the native path in the Pyleecan folder.

Nevertheless, the efficiency results look unfeasible to me.
When plotting the OP after each simulation, I see still d current is applied, while I want a pure q current excitation for this SPM motor.
out = simu.run()
print("operation point:")
print(out.elec.OP)

Forcing the Id_ref to 0 seems not to work
Id_ref = 0*(I0_rmsexp(1jPhi0)).real # d axis current [Arms]
Iq_ref = (I0_rmsexp(1jPhi0)).imag # q axis current [Arms]
simu.input.OP = OPdq(N0=N0, Id_ref=Id_ref, Iq_ref=Iq_ref)

[SebGue]

Hello @KizzieNova ,

in your modified code I used the Pyleecan folder as a result directory. So that makes sense.

I guess ElectLUTdq electrical simulation has a MTPA solver, so you may get a d-current regardless of you Id_ref values.
If you want to do a simulation with reference current values for d and q axis, you have to use Electrical instead of ElecLUTdq.
I.e. you calcuate your LUT with ElecLUTdq as before and setup an Electrical simulation for your main calculation loop.

Best regards, Sebastian