Engee documentation
Notebook

Parameterization of SDM based on measurements

This example shows how to determine the parameters of a permanent magnet synchronous motor based on experimental measurements.

In order to determine the parameters of the SDPM engine, we will conduct three experiments. The desired values of the engine parameters are known to us in advance, and we will compare them with the values obtained from the experiments.

Data from the engine documentation:

In [ ]:
nPolespairs = 2;        # Количество пар полюсов
iRated = 2.5;           # Ток при номинальной скорости и нагрузке (А)
rpmRated = 7500;        # Номинальная скорость (об/мин)
torqueRated = 0.07;     # Крутящий момент при номинальной скорости (Н·м)

torqueStall = 0.28;     # Заданный крутящий момент (Н·м)
iStall = 8.5;           # Пиковый ток (А)

inertia = 5e-6;                        # Момент инерции (кг·м2)
viscousDamping = 1.13e-6;              # Коэффициент вязкого трения
staticFriction = 7e-4;                 # Коэффициент сухого трения
coulombFriction = 0.8*staticFriction;  # Коэффициент кулоновского трения

R = 3.43;               # Сопротивление (Ом)
L = 0.53e-3;            # Индуктивность (Гн)

# Параметры блока PMSM (СДПМ)
torqueConstant = 2/3*torqueStall/iStall; # Постоянная момента (Нм/А)
emfConstant = torqueConstant; # Постоянная противо-ЭДС (В∙с/рад)
PM = torqueConstant/nPolespairs; # Потокосцепление постоянного магнита (Вб)

# Параметры и флаги контроллера
rpmDemand = rpmRated; # Заданная скорость (об/мин)
enableDrive = 1; # Включение управления
dt = 1e-4; # Период дискретизации контроллера (с)

Test 1. Determination of R and L with a stationary rotor

The first experiment is carried out on an engine with a locked rotor. A voltage step is applied between the phases and the transient response process is studied. The resulting time constant τ is determined by the values of the resistance and inductance of the stator:

image_2.png

Let's set the initial rotation angle of the rotor:

In [ ]:
rotorangle = -90/(nPolespairs);

Launching the model:

In [ ]:
engee.addpath(@__DIR__)

if "PMSMFromMeasurement1" in [m.name for m in engee.get_all_models()]
    m = engee.open( "PMSMFromMeasurement1" ) # загрузка модели
else
    m = engee.load( "PMSMFromMeasurement1.engee" )
end
results1 = engee.run(m, verbose=true)
Building...
Progress 0%
Progress 12%
Progress 64%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    run_id => 91,
    "ia" => WorkspaceArray{Float64}("PMSMFromMeasurement1/ia")
,
    "uab" => WorkspaceArray{Float64}("PMSMFromMeasurement1/uab")

)

Reading data from a model:

In [ ]:
t = results1["uab"].time;       # Время
uab_1 = results1["uab"].value;  # Напряжение между фазами a и b
ia  = results1["ia"].value;     # Ток на фазе a

Calculate the resistance from the values of steady-state voltage and current.
The windings a and b are connected in series in the experiment, so we divide the resulting resistance by 2:

In [ ]:
R_e = uab_1[end] / ia[end] / 2;

Let's find the time constant of the circuit - this is the time it takes for the current to reach the value from the established:

In [ ]:
idx = findfirst(x -> x >= (0.63 * ia[end]), ia);
tau = t[idx];

From the formula let's express and calculate the inductance L:

In [ ]:
L_e = R_e * tau;

The expected transition process:

In [ ]:
iEstimated = uab_1[end] / (2 * R_e) * (1 .- exp.(-t .* R_e ./ L_e));

In the graph below, you can compare the measured current from the model and the current graph obtained from the calculations.

In [ ]:
using Plots
plotlyjs();
plot(t, ia, xlabel="Время, c", ylabel="Ток, А", w = 2, label="Измеренный ток", linecolor =:blue)
plot!(t, iEstimated, xlabel="Время, c", ylabel="Ток, А", w = 2, label="Ожидаемый ток", linecolor =:green)
plot!([tau, tau], [0, ia[end]], linecolor =:red, line =:dashdot, label="τ" )
xlims!(0, 10*tau)   # Установка пределов оси x
ylims!(0, ia[end])  # Установка пределов оси y

annotate!(0.0008, 4.5, text(string("R эталонное = ", R, " Ом"), :left, 10))
annotate!(0.0008, 4, text(string("R измеренное  = ", round(R_e, digits=2), " Ом"), :left, 10))
annotate!(0.0008, 3.5, text(string("L эталонная = ", L*1000, " мГн"), :left, 10))
annotate!(0.0008, 3, text(string("L измеренная  = ", round(L_e*1000, digits=2), " мГн"), :left, 10))
annotate!(tau, 1, text(string("Измеренная постоянная времени  = ", round(1000*tau, digits=3), " мc"), :left, 10))
Out[0]:

Test 2. Determination of the constant counter-EMF

In the second experiment, the motor rotates without electrical load. This allows us to estimate the constant counter-EMF.

Note that the counter-EMF constant, expressed in SI units, is equal to the torque constant of the motor, so only one constant is calculated in this example.

image_2.png

Launching the model:

In [ ]:
if "PMSMFromMeasurement2" in [m.name for m in engee.get_all_models()]
    m = engee.open( "PMSMFromMeasurement2" ) # загрузка модели
else
    m = engee.load( "PMSMFromMeasurement2.engee" )
end
results2 = engee.run(m, verbose=true)
Building...
Progress 0%
Progress 34%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    run_id => 92,
    "uab" => WorkspaceArray{Float64}("PMSMFromMeasurement2/uab")

)

Reading data from a model:

In [ ]:
t = results2["uab"].time;       # Время
uab_2  = results2["uab"].value; # Напряжение  a-b

Determination of the counter-EMF constant

The counter-EMF constant it can be calculated using the formula:

where – voltage drop in phase a of the motor,

– angular velocity of the motor shaft.

The torque constant is equal to the counter-EMF constant expressed in SI units.

In [ ]:
wRef = rpmRated*2*pi/60;        # Угловая скорость ротора (рад/с)
vPeakLL = maximum(abs.(uab_2)); # Пиковое линейное напряжение (В)
vPeakLN = vPeakLL/sqrt(3);      # Пиковое напряжение фаза-нейтраль (В)
emfConstantE = vPeakLN/wRef;    # wq константы противо-ЭДС

The graph below shows the counter-EMF that occurs when the experimental engine is moving at rated speed.

The number of pole pairs can be determined by calculating the number of periods of current change per rotation of the rotor.

In [ ]:
plot(t, uab_2, xlabel="Время, c", ylabel="Противо-ЭДС", label="Противо-ЭДС экспериментального двигателя", w = 2)
annotate!(0.0013, -20, text(string("Константа противо-ЭДС эталонная = ", round(emfConstant,digits=4), " В∙с/рад"), :left, 10))
annotate!(0.0013, -25, text(string("Константа противо-ЭДС измеренная = ", round(emfConstantE,digits=4), " В∙с/рад"), :left, 10))
Out[0]:

3. Determination of coefficients of friction and moment of inertia

In the third experiment, the unloaded motor is controlled by a controller.
In a previous experiment, a constant of the moment was obtained. With its help, by converting the stator currents, friction coefficients can be calculated.

image_2.png

Variables required for controller simulation:

In [ ]:
rpm2rad = 2*pi/60;
rad2rpm = 60/(2*pi);
powermax = torqueRated*rpmRated*2*pi/60;
Kp = 0.01;
Ki = 2.0;
ts = 0.0001;
t2eq_Iq = 2/(3*nPolespairs*PM);

Let's measure the mechanical torque required to maintain a constant rotational speed at four different speeds.
The torque required at lower speeds is primarily aimed at overcoming dry friction, at higher speeds - viscous friction.
The graph shows the measured moments for the four speeds. A straight line is drawn through these points. The zero-velocity intersection gives the dry friction coefficient, and the slope gives the viscous friction coefficient.

In [ ]:
function trapz(x, y)
    n = length(x)
    integral = 0.0
    for i in 1:n-1
        integral += (x[i+1] - x[i]) * (y[i+1] + y[i]) / 2
    end
    return integral
end
rpmVec = [0.25, 0.5, 0.75, 1.0]*rpmRated;
trqVec = zeros(size(rpmVec));
for i=1:length(rpmVec)
    if "PMSMFromMeasurement3" in [m.name for m in engee.get_all_models()]
        m = engee.open( "PMSMFromMeasurement3" ) # загрузка модели
    else
        m = engee.load( "PMSMFromMeasurement3.engee" )   
    end
    rpmDemand = rpmVec[i]; 
    engee.set_param!("PMSMFromMeasurement3", "StopTime" => 2*60/rpmDemand)
    results3 = engee.run(m, verbose=true)
    t = results3["ia"].time;        # Время из модели
    ia_3  = results3["ia"].value;   # ток фазы-a из модели
    idx = findfirst(x -> x >= (60/rpmDemand), t); 
    temp = trapz(t[idx:end],(ia_3[idx:end]).^2)
    iaRms = sqrt(temp/(t[end]-t[idx]));
    trqVec[i] = 3*torqueConstant*iaRms;
end
rpmDemand = rpmRated;
Building...
Progress 0%
Progress 6%
Progress 14%
Progress 20%
Progress 25%
Progress 30%
Progress 36%
Progress 42%
Progress 47%
Progress 53%
Progress 59%
Progress 64%
Progress 70%
Progress 76%
Progress 83%
Progress 91%
Progress 97%
Progress 100%
Progress 100%
Building...
Progress 0%
Progress 7%
Progress 26%
Progress 37%
Progress 49%
Progress 59%
Progress 69%
Progress 79%
Progress 89%
Progress 99%
Progress 100%
Progress 100%
Building...
Progress 0%
Progress 5%
Progress 27%
Progress 47%
Progress 61%
Progress 76%
Progress 89%
Progress 100%
Progress 100%
Building...
Progress 0%
Progress 10%
Progress 47%
Progress 69%
Progress 91%
Progress 100%
Progress 100%
In [ ]:
scatter(rpmVec, 1000*trqVec, title="Механические потери крутящего момента", xlabel="Скорость, м/c", ylabel="Момент, мН*м", w = 2, label="Измеренные значения", legend=:bottomright )
using Polynomials
coef = coeffs(fit(2*pi/60*rpmVec, trqVec,1));
estimatedStaticFriction = coef[1];
estimatedViscousDamping = coef[2];
rpmVec0 = [0; rpmVec];
plot!(rpmVec0, 1000*(estimatedStaticFriction.+rpmVec0*2*pi/60*estimatedViscousDamping), label="Линейная аппроксимация")
annotate!(trqVec[end], 2.1, text(string("Коэффициент сухого трения эталонный = ", staticFriction), :left, 10))
annotate!(trqVec[end], 2.0, text(string("Коэффициент сухого трения измеренный = ", round(estimatedStaticFriction,digits=6)), :left, 10))
annotate!(trqVec[end], 1.9, text(string("Коэффициент вязкого трения эталонный = ", viscousDamping), :left, 10))
annotate!(trqVec[end], 1.8, text(string("Коэффициент вязкого трения измеренный = ", round(estimatedViscousDamping, digits=9)), :left, 10))
Out[0]:

The second graph shows an engine deceleration test. For this test, the required torque is set to zero or the motor is switched off.
Using the measured deceleration, at the specified values of the friction and damping moments, it is possible to determine the moment of inertia of the motor rotor.

In [ ]:
enableDrive = 0; 
engee.set_param!("PMSMFromMeasurement3", "StopTime" => 2*60/rpmDemand)
if "PMSMFromMeasurement3" in [m.name for m in engee.get_all_models()]
    m = engee.open( "PMSMFromMeasurement3" ) # загрузка модели
else
    m = engee.load( "PMSMFromMeasurement3.engee" )
end 
results4 = engee.run(m, verbose=true) 
Building...
Progress 0%
Progress 6%
Progress 50%
Progress 69%
Progress 87%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    run_id => 97,
    "ia" => WorkspaceArray{Float64}("PMSMFromMeasurement3/Controller/Sensing_iabc/ia")
,
    "rpm" => WorkspaceArray{Float64}("PMSMFromMeasurement3/Controller/Encoder/rpm")

)
In [ ]:
enableDrive = 1;
t = results4["rpm"].time;       # Время из модели
rpm  = results4["rpm"].value;   # скорость ротора
plot(t, rpm, xlabel="Время, c", ylabel="Скорость (об/мин)", label="Тест замедления", w = 2)
initialTorque = -staticFriction - rpmRated*2*pi/60*viscousDamping;
initialAccel = 2*pi/60*(rpm[end]-rpm[1])/(2*60/rpmDemand);
estimatedInertia = initialTorque/initialAccel;
annotate!(0, 7465, text(string("Момент инерции эталонный  = ", inertia, " кг∙м^2"), :left, 10))
annotate!(0, 7460, text(string("Момент инерции измеренный = ", round(estimatedInertia,digits=6), " кг∙м^2"), :left, 10))
Out[0]:

Conclusion

In this example, we have determined the parameters of a permanent magnet synchronous motor based on experimental measurements.
To do this, we conducted three tests. In the first test, we found the resistance and inductance of the motor winding, in the second we determined the counter-EMF constant, and in the third we determined the coefficients of friction and moment of inertia.
The values found experimentally are comparable to the engine parameters from the documentation.