Документация Engee
Notebook

Параметризация СДПМ на основе измерений

В этом примере показано, как определить параметры синхронного двигателя с постоянными магнитами (СДПМ) на основе экспериментальных измерений.

Для того чтобы определить параметры двигателя СДПМ, мы проведем три эксперимента. Искомые значения параметров двигателя известны нам заранее, и мы сравним их со значениями, полученными из экспериментов.

Данные из документации на двигатель:

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; # Период дискретизации контроллера (с)

Тест 1. Определение R и L при неподвижном роторе

Первый эксперимент проводится на двигателе c заблокированным ротором. Между фазами подается ступенька напряжения и изучается переходный процесс отклика. Результирующая постоянная времени τ определяется значениями сопротивления и индуктивности статора:

$$τ = \frac{R}{L}$$

image_2.png

Зададим начальный угол поворота ротора:

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

Запуск модели:

In [ ]: 
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 19%
Progress 100%
Out[0]:
SimulationResult(
    "ia" => WorkspaceArray{Float64}("PMSMFromMeasurement1/ia")
,
    "uab" => WorkspaceArray{Float64}("PMSMFromMeasurement1/uab")

)

Чтение данных из модели:

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

Рассчитаем сопротивление из значений, установившихся напряжения и тока. Обмотки a и b в эксперименте соединены последовательно, поэтому разделим полученное сопротивление на 2:

$$R_e = \frac{1}{2}\frac{u_{ab}}{i_a}$$

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

Найдем постоянную времени цепи - это время, за которое ток достигает значения $1-\frac{1}{e}=0.63$ от установившегося:

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

Из формулы $τ = \frac{R}{L}$ выразим и рассчитаем индуктивность L:

In [ ]: 
L_e = R_e * tau;

Ожидаемый переходный процесс:

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

На графике ниже можно сравнить измеренный ток из модели и график тока, полученный из расчетов.

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]:

Тест 2. Определение постоянной противо-ЭДС

Во втором эксперименте двигатель вращается без электрической нагрузки. Это позволяет оценить постоянную противо-ЭДС.

Обратите внимание, что постоянная противо-ЭДС, выраженная в единицах СИ, равна постоянной момента двигателя, поэтому в этом примере рассчитывается только одна константа.

image_2.png

Запуск модели:

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 18%
Progress 100%
Out[0]:
SimulationResult(
    "uab" => WorkspaceArray{Float64}("PMSMFromMeasurement2/uab")

)

Чтение данных из модели:

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

Определение константы противо-ЭДС

Константу противо-ЭДС $K_e$ можно рассчитать по формуле:

$$K_e= \frac{V_{a, pk}}{\omega},$$

где $V_{a,pk}$ – падение напряжения на фазе a двигателя,

$\omega$ – угловая скорость вала двигателя.

Постоянная момента равна постоянной противо-ЭДС, выраженной в единицах СИ.

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

На графике ниже показана противо-ЭДС, возникающая при движении экспериментального двигателя на номинальной скорости.

Количество пар полюсов можно определить, подсчитав количество периодов изменения тока за один оборот ротора.

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. Определение коэффициентов трения и момента инерции

В третьем эксперименте ненагруженный двигатель управляется контроллером. В прошлом эксперименте была получена постоянная момента. C её помощью, путем преобразования токов статора, могут быть рассчитаны коэффициенты трения.

image_2.png

Переменные, требуемые для моделирования контроллера:

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);

Измерим механический крутящий момент, необходимый для поддержания постоянной частоты вращения на четырех различных скоростях. Момент, необходимый на более низких скоростях, в первую очередь направлен на преодоление сухого трения, на более высоких скоростях - вязкого трения. На графике показаны измеренные моменты для четырех скоростей. Через эти точки проведена прямая. Пересечение нулевой скорости дает коэффициент сухого трения, а наклон - коэффициент вязкого трения.

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 5%
Progress 11%
Progress 17%
Progress 23%
Progress 28%
Progress 34%
Progress 40%
Progress 49%
Progress 54%
Progress 59%
Progress 66%
Progress 72%
Progress 78%
Progress 99%
Progress 100%
Building...
Progress 0%
Progress 6%
Progress 17%
Progress 28%
Progress 38%
Progress 48%
Progress 58%
Progress 67%
Progress 78%
Progress 88%
Progress 97%
Progress 100%
Building...
Progress 0%
Progress 9%
Progress 28%
Progress 50%
Progress 66%
Progress 88%
Progress 100%
Building...
Progress 0%
Progress 11%
Progress 33%
Progress 52%
Progress 71%
Progress 94%
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]:

На втором графике показан тест замедления двигателя. Для этого теста требуемый крутящий момент устанавливается равным нулю или двигатель отключается. Используя измеренное замедление, при заданных значениях моментов трения и демпфирования, можно определить момент инерции ротора двигателя.

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 37%
Progress 71%
Progress 89%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    "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]:

Вывод

В этом примере мы определили параметры синхронного двигателя с постоянными магнитами на основе экспериментальных измерений. Для этого мы провели три теста. В первом тесте мы нашли сопротивление и индуктивность обмотки двигателя, во втором определили константу противо-ЭДС, а в третьем – коэффициенты трения и момент инерции. Найденные экспериментально значения сопоставимы с параметрами двигателя из документации.

Блоки, использованные в примере