Моделирование САР скорости асинхронного двигателя с переменными is – ψr в Matlab-Script в системе относительных единиц



Моделирование САР скорости асинхронного двигателя с переменными i_s – ψ_r в Matlab-Script в системе относительных единиц

Емельянов Александр Александрович, доцент;

Гусев Владимир Михайлович, магистрант;

Пестеров Дмитрий Ильич, студент;

Даниленко Дмитрий Сергеевич, студент

Российский государственный профессионально-педагогический университет (г. Екатеринбург)

Бесклеткин Виктор Викторович, магистрант.

Уральский федеральный университет имени первого Президента России Б. Н. Ельцина (г. Екатеринбург)

Иванин Александр Юрьевич, техник-метролог.

НПО «НТЭС» (Республика Татарстан, г. Бугульма)

В работе [1] приведена модель САР скорости асинхронного двигателя в Simulink. В этой статье покажем поэтапное преобразование всех элементов САР скорости в Matlab-Script. На рис. 1 приводим всю систему, в которой даны модель асинхронного двигателя (номер 7), в контурах тока по проекциям x и y соответствующие ПИ-регуляторы тока (номера 4 и 6), в контуре скорости П-регулятор скорости (номер 1).

Рис. 1. Математическая модель САР скорости асинхронного двигателя

Важным элементом является контур потока с ПИ-регулятором потока (номер 2). Для ориентации системы координат по потокосцеплению ротора вводится наблюдатель (номер 8). В модели учтена компенсация перекрестных связей (номер 5). Сигнал задания по скорости выполнен на задатчике интенсивности. В цепи задания скорости перед регулятором скорости предусмотрен фильтр.

Алгоритм перевода всех элементов САР скорости:

‒ приводится математическая формула той или иной переменной, выраженной в Simulink;

‒ приводится его структурная схема;

‒ переход от изображений к оригиналу (от s к d/dt) и решение с помощью простого метода Эйлера.

Математическая модель асинхронного двигателя с переменными i_s – ψ_r

А) Выражение для статорного тока i_sx по проекции x, подготовленное для структурной схемы, имеет следующий вид [1]:

(1)

где - электрическая скорость вращения ротора;

- механическая угловая скорость на валу двигателя.

Структурная схема (рис. 2).

Рис. 2. Структурная схема для определения тока i_sx в Simulink

Преобразуем уравнение (1) для программирования в Matlab-Script:

Обозначим , тогда:

Переходим к оригиналу :

Переходим к конечным разностям (простой метод Эйлера):

Отсюда ток i_sx в Matlab-Script определится следующим образом:

Б) Уравнение для определения тока i_sy в Simulink, полученное в работе [1], имеет следующий вид:

(2)

Структурная схема реализации уравнения (2) приведена на рис. 3.

Рис. 3. Структурная схема для определения тока i_sy в Simulink

Аналогично преобразуем выражение тока i_sy в форму, удобную для программирования в Matlab-Script:

Переходим к оригиналу:

Переходим к конечным разностям:

Ток i_sy в Matlab-Script определится следующим образом:

В) Уравнение для определения потокосцепления ψ_rxв Simulink имеет следующий вид:

(3)

Структурная схема представлена на рис. 4.

Рис. 4. Структурная схема для определения потокосцепления ψ_rx в Simulink

Преобразуем уравнение (3) для программирования в Matlab-Script:

Обозначим , тогда:

Переходим к оригиналу:

Переходим к конечным разностям:

Отсюда потокосцепление ψ_rx в Matlab-Script определится следующим образом:

Г) Уравнение для определения потокосцепления ψ_ry в Simulink имеет следующий вид:

(4)

Структурная схема реализации уравнения (4) приведена на рис. 5.

Рис. 5. Структурная схема для определения потокосцепления ψ_ry в Simulink

Преобразуем выражение потокосцепления ψ_ry в форму, удобную для программирования в Matlab-Script:

Д) На рис. 6 представлена структурная схема для реализации уравнения электромагнитного момента в Simulink:

Рис. 6. Математическая модель определения электромагнитного момента m в Simulink

Уравнение электромагнитного момента для реализации в Matlab-Script:

Е) Механическая угловая скорость вращения вала двигателя в Simulink (рис. 7):

Рис. 7. Математическая модель определения механической угловой скорости вращения вала двигателя в Simulink

Отсюда механическая угловая скорость вращения вала двигателя в Matlab-Script:

Ж) Электрическая скорость вращения ротора в Simulink (рис. 8):

Рис. 8. Математическая модель определения электрической скорости вращения ротора в Simulink

Электрическая скорость вращения ротора в Matlab-Script:

Реализация математической модели асинхронного двигателя с короткозамкнутым ротором с переменными i_s – ψ_r в Matlab-Script в системе относительных единиц приведена в листинге 1.

Листинг 1

% Номинальные данные

PN=320000; UsN=380; IsN=324; fN=50; Omega0N=104.7;

OmegaN=102.83; nN=0.944; cos_phiN=0.92; zp=3;

% Параметры Т-образной схемы замещения при номинальной частоте

Rs=0.0178; Xs=0.118; Rr=0.0194; Xr=0.123; Xm=4.552; J=28;

% Базисные величины системы относительных единиц

Ub=sqrt(2)*UsN;

Ib=sqrt(2)*IsN;

OmegasN=2*pi*fN;

Omegab=OmegasN;

Omegarb=Omegab/zp;

Zb=Ub/Ib;

kd=1.0084;

Mb=kd*PN/OmegaN;

Pb=Mb*Omegarb;

% Расчет коэффициентов АД

rs=Rs/Zb;

lbs=Xs/Zb;

lbr=Xr/Zb;

lm=Xm/Zb;

Tj=J*Omegarb/Mb;

betaN=(Omega0N-OmegaN)/Omega0N;

SsN=3*UsN*IsN;

ZetaN=SsN/Pb;

kr=lm/(lm+lbr);

lbe=lbs+lbr+lbs*lbr*lm^(-1);

roN=0.9962;

rrk=roN*betaN;

Tr=lm/(rrk*kr);

Tr1=Tr/Omegab;

re=rs+rrk*kr^2;

Te=kr*lbe/re;

Te1=Te/Omegab;

% Расчет асинхронного двигателя (номер 7)

K=input('Длительность цикла k=');

for k=1:(K+1)

dt=0.000001;

usx(k+1)=0;usy(k+1)=1; wk(k)=1;

isx(1)=0; isy(1)=0; psirx(1)=0; psiry(1)=0;

wm(1)=0; w(1)=0; mc=0;

% Ток isx (А)

isx(k+1)=isx(k)+(-isx(k)+(1/re)*usx(k+1)+rrk*(kr^2)/(re*lm)*psirx(k)+ (kr/re)*w(k)*psiry(k)+(kr*lbe/re)*wk(k)*isy(k))*dt/Te1;

% Ток isy (Б)

isy(k+1)=isy(k)+(-isy(k)+(1/re)*usy(k+1)+rrk*(kr^2)/(re*lm)*psiry(k)-(kr/re)*w(k)*psirx(k)-(kr*lbe/re)*wk(k)*isx(k))*dt/Te1;

% Поток psirx (В)

psirx(k+1)=psirx(k)+(-psirx(k)+lm*isx(k)+(lm/(rrk*kr))*(wk(k)-w(k))*psiry(k))*dt/Tr1;

% Поток psiry (Г)

psiry(k+1)=psiry(k)+(-psiry(k)+lm*isy(k)-(lm/(rrk*kr))*(wk(k)-w(k))*psirx(k))*dt/Tr1;

% Электромагнитныймомент (Д)

m(k+1)=ZetaN*kr*(psirx(k+1)*isy(k+1)-psiry(k+1)*isx(k+1));

% Механическая скорость (Е)

wm(k+1)=wm(k)+(m(k+1)-mc)*dt/Tj;

% Электрическая скорость (Ж)

w(k+1)=wm(k+1)*zp;

end;

Математическое моделирование регуляторов тока

В работе [1] была получена передаточная функция для регуляторов тока по проекциям x и y:

гдеT_μ - некомпенсируемая постоянная времени (примем T_μ = 0,0025 с).

Обозначим:

Математические модели ПИ-регуляторов тока по проекциям x и y в Simulink приведены на рис. 9 и 10. Преобразуем их для программирования в Matlab-Script.

Рис. 9. ПИ-регулятор тока по проекции x в Simulink

Рис. 10. ПИ-регулятор тока по проекции y в Simulink

Пропорциональная часть регулятора тока по оси x в Simulink:

Выразим пропорциональную часть в Matlab-Script:

где

Интегральная часть регулятора тока по оси x:

Переходим от изображения к оригиналу:

Выразим интегральную часть через конечные разности:

Уравнение напряжения задания на выходе регулятора тока по оси x будет иметь следующий вид:

Аналогично преобразуем регулятор тока по оси y.

Пропорциональная часть:

где

Интегральная часть:

Уравнение на выходе регулятора тока по оси y:

Реализация математической модели регуляторов тока в Matlab-Script представлена в листинге 2.

Листинг 2

Tm=0.0025; dt=0.000001;

Ki=Te1/(2*Tm/re); Ti=2*Tm/re;

isx(1)=0; isy(1)=0;ux2(1)=0;uy2(1)=0;

ixsum(k+1)=ixzad(k+1)-isx(k);

iysum(k+1)=iyzad(k+1)-isy(k);

% Регулятор тока по оси x

%Пропорциональная часть задания usx

ux1(k+1)=ixsum(k+1)*Ki;

%Интегральнаячастьзадания usx

ux2(k+1)=ux2(k)+ixsum(k+1)*dt/Ti;

%Задание usx

uxzad(k+1)=ux1(k+1)+ux2(k+1);

% Регулятор тока по оси y

%Пропорциональная часть задания usy

uy1(k+1)=iysum(k+1)*Ki;

%Интегральная часть задания usy

uy2(k+1)=uy2(k)+iysum(k+1)*dt/Ti;

%Задание usy

uyzad(k+1)=uy1(k+1)+uy2(k+1);

Математическое моделирование наблюдателя потокосцепления ротора

Модель наблюдателя потокосцепления ротора в Simulink, полученная в работе [1], приведена на рис. 11. Преобразуем эту модель в Matlab-Script.

Рис. 11. Модель наблюдателя потокосцепления ротора в Simulink

Приведем уравнение модуля потокосцепления ротора к оригиналу:

Переходим к конечным разностям:

Уравнение скольжения для программирования в Matlab-Script будет иметь вид [1], [2], [3]:

Отсюда угловая скорость вращения системы координат :

Математическая модель наблюдателя в Matlab-Script приведена в листинге 3.

Листинг 3

dt=0.000001; psirx_oc(1)=0.001;

% Модуль потокосцепления ротора

psirx_oc(k+1)=psirx_oc(k)+(-psirx_oc(k)+lm*isx(k+1))*dt/Tr1;

% Скольжение

beta_psir(k+1)=isy(k+1)*rrk*kr/psirx_oc(k+1);

% Угловая скорость вращения системы координат

wk(k+1)=beta_psir(k+1)+w(k+1);

Математическое моделирование регулятора потока

Модель ПИ-регулятора потока в Simulink дана на рис. 12.

Рис. 12. ПИ-регулятор потока в Simulink

Номинальное потокосцепление ротора в соответствии с [3] определяется по следующей формуле и при векторном управлении поддерживается постоянным:

Передаточная функция регулятора потока из работы [1]:

где n = 2.

Выразим коэффициенты ПИ-регулятора потока:

Определим пропорциональную часть:

где

Интегральная часть регулятора потока:

Переходим от изображения к оригиналу:

Выразим интегральную часть через конечные разности:

Определим задание тока на выходе регулятора потока в Matlab-Script:

Реализация математической модели регулятора потока в Matlab-Script приведена в листинге 4.

Листинг 4

Tm=0.0025; psirN=0.942; n=2; dt=0.000001;

psirx_oc(1)=0.001; ixzad2(1)=0;

Kpsi=Tr1/(4*n*Tm*lm);

Tpsi=4*n*Tm*lm;

psirxsum(k+1)=psirN-psirx_oc(k);

% Пропорциональная часть задания isx

ixzad1(k+1)=psirxsum(k+1)*Kpsi;

% Интегральная часть задания isx

ixzad2(k+1)=ixzad2(k)+psirxsum(k+1)*dt/Tpsi;

% Задание isx

ixzad(k+1)=ixzad1(k+1)+ixzad2(k+1);

Математическое моделирование регулятора скорости

Математическая модель П-регулятора скорости в Simulink [1] дана на рис. 13.

Рис. 13. Пропорциональный регулятор скорости в Simulink

Передаточная функция регулятора скорости:

Отсюда определим задание момента :

где

Математическая модель регулятора скорости в Matlab-Script представлена в листинге 5.

Листинг 5

Tm=0.0025; w(1)=0;

wsum(k+1)=wzad1(k+1)-w(k);

% Задание момента m

mzad(k+1)=wsum(k+1)*Tj/(4*Tm);

Математическое моделирование компенсации перекрестных связей

Математическая модель компенсации перекрестных связей в Simulink [1] дана на рис. 14.

Рис. 14. Компенсация внутренних перекрестных связей в Simulink

Компенсационные составляющие каналов управления определятся следующим образом:

Реализация математической модели компенсации перекрестных связей в Matlab-Script представлена в листинге 6.

Листинг 6

isx(1)=0; isy(1)=0; psirx_oc(1)=0.001; wk(1)=0;

% Звенокомпенсации x

ukx(k+1)=-wk(k)*kr*lbe*isy(k);

% Звенокомпенсации y

uky(k+1)=wk(k)*kr*(lbe*isx(k)+psirx_oc(k));

Математическое моделирование задатчика интенсивности

Задание на скорость ω* в Simulink формируется в блоке Signal Builder (рис. 15).

Рис. 15. Сигнал задания на скорость ω* в Simulink

Программирование сигнала задания на скорость в Matlab-Script представлено в листинге 7.

Листинг 7

tn=0.8;

tk=1.1;

dt=0.000001;

% Задание на скорость

if((k*dt>=0)&&(k*dt<=tn))

wzad(k+1)=0;

end;

if((k*dt>=tn)&&(k*dt<=tk))

wzad(k+1)=(k*dt-tn)/(tk-tn);

end;

if(k*dt>tk)

wzad(k+1)=1;

end;

Математическое моделирование задания по скорости на выходе фильтра

Передаточная функция фильтра:

Определим задание скорости на выходе фильтра:

Перейдем от изображения к оригиналу:

Переходим к конечным разностям:

Математическая модель задания скорости на выходе фильтра в Matlab-Script дана в листинге 8.

Листинг 8

dt=0.000001;

Tm1=0.0075;

wzad1(1)=0;

% Задание скорости на выходе фильтра

wzad1(k+1)=wzad1(k)+(wzad(k+1)-wzad1(k))*dt/Tm1;

Математическое моделирование задания статорного тока по проекции y

Математическая модель задания тока в Simulink дана на рис. 16.

Рис. 16. Реализация задания статорного тока в Simulink

Задание на статорный ток по проекции y:

Математическая модель задания в Matlab-Script представлена в листинге 9.

Листинг 9

psirx_oc(1)=0.001;

% Задание isy

iyzad(k+1)=mzad(k+1)/(psirx_oc(k)*kr);

Моделирование САР скорости асинхронного двигателя

Полная математическая модель САР скорости асинхронного двигателя в Matlab-Script приведена в листинге 10.

Листинг 10

% Номинальные данные АД

PN=320000; UsN=380; IsN=324; fN=50; Omega0N=104.7;

OmegaN=102.83; nN=0.944; cos_phiN=0.92; zp=3;

% Параметры Т-образной схемы замещения при номинальной частоте

Rs=0.0178; Xs=0.118; Rr=0.0194; Xr=0.123; Xm=4.552; J=28;

% Базисные величины системы относительных единиц

Ub=sqrt(2)*UsN;

Ib=sqrt(2)*IsN;

OmegasN=2*pi*fN;

Omegab=OmegasN;

Omegarb=Omegab/zp;

Zb=Ub/Ib;

kd=1.0084;

Mb=kd*PN/OmegaN;

Pb=Mb*Omegarb;

% Расчет коэффициентов АД

rs=Rs/Zb;

lbs=Xs/Zb;

lbr=Xr/Zb;

lm=Xm/Zb;

Tj=J*Omegarb/Mb;

betaN=(Omega0N-OmegaN)/Omega0N;

SsN=3*UsN*IsN;

ZetaN=SsN/Pb;

kr=lm/(lm+lbr);

lbe=lbs+lbr+lbs*lbr*lm^(-1);

roN=0.9962;

rrk=roN*betaN;

Tr=lm/(rrk*kr);

Tr1=Tr/Omegab;

re=rs+rrk*kr^2;

Te=kr*lbe/re;

Te1=Te/Omegab;

% Параметры САР скорости

Tm=0.0025; Tm1=0.0075;

Ki=Te1/(2*Tm/re);

Ti=2*Tm/re;

n=2;

Kpsi=Tr1/(4*n*Tm*lm);

Tpsi=4*n*Tm*lm;

psirN=0.942;

tn=0.8; tk=1.1;

dt=0.000001;

% Расчет САР скорости АД

K=input('Длительность цикла k=');

% Параметры САР скорости в начальный момент времени

isx(1)=0; isy(1)=0; psirx(1)=0; psiry(1)=0; wm(1)=0; mc=0;

psirx_oc(1)=0.001; ixzad2(1)=0; ux2(1)=0; uy2(1)=0;

wk(1)=0; w(1)=0; wzad1(1)=0;

% Задание на скорость

for k=1:(K+1)

if((k*dt>=0)&&(k*dt<=tn))

wzad(k+1)=0;

end;

if((k*dt>=tn)&&(k*dt<=tk))

wzad(k+1)=(k*dt-tn)/(tk-tn);

end;

if(k*dt>tk)

wzad(k+1)=1;

end;

% Задание скорости на выходе фильтра

wzad1(k+1)=wzad1(k)+(wzad(k+1)-wzad1(k))*dt/Tm1;

% Моделирование регулятора потока (номер 2)

psirxsum(k+1)=psirN-psirx_oc(k);

% Пропорциональная часть задания isx

ixzad1(k+1)=psirxsum(k+1)*Kpsi;

% Интегральная часть задания isx

ixzad2(k+1)=ixzad2(k)+psirxsum(k+1)*dt/Tpsi;

% Задание isx

ixzad(k+1)=ixzad1(k+1)+ixzad2(k+1);

% Моделирование регулятора скорости (номер 1)

wsum(k+1)=wzad1(k+1)-w(k);

% Задание момента m

mzad(k+1)=wsum(k+1)*Tj/(4*Tm);

% Задание isy (номер 3)

iyzad(k+1)=mzad(k+1)/(psirx_oc(k)*kr);

% Моделирование регуляторов тока (номера 4 и 6)

ixsum(k+1)=ixzad(k+1)-isx(k);

iysum(k+1)=iyzad(k+1)-isy(k);

% Регулятор тока по оси x (номер 4)

%Пропорциональная часть задания usx

ux1(k+1)=ixsum(k+1)*Ki;

%Интегральнаячастьзадания usx

ux2(k+1)=ux2(k)+ixsum(k+1)*dt/Ti;

%Задание usx

uxzad(k+1)=ux1(k+1)+ux2(k+1);

% Регулятор тока по оси y (номер 6)

%Пропорциональная часть задания usy

uy1(k+1)=iysum(k+1)*Ki;

%Интегральная часть задания usy

uy2(k+1)=uy2(k)+iysum(k+1)*dt/Ti;

%Задание usy

uyzad(k+1)=uy1(k+1)+uy2(k+1);

% Моделирование звена компенсации (номер 5)

% Звено компенсации x

ukx(k+1)=-wk(k)*kr*lbe*isy(k);

% Звено компенсации y

uky(k+1)=wk(k)*kr*(lbe*isx(k)+psirx_oc(k));

% Моделирование напряжений usx и usy

usx(k+1)=uxzad(k+1)-ukx(k+1);

usy(k+1)=uyzad(k+1)+uky(k+1);

% Моделирование асинхронного двигателя (номер 7)

% Ток isx (А)

isx(k+1)=isx(k)+(-isx(k)+(1/re)*usx(k+1)+(rrk*(kr^2)/(re*lm))*psirx(k)+(kr/re)*w(k)*psiry(k)+(kr*lbe/re)*wk(k)*isy(k))*dt/Te1;

% Ток isy (Б)

isy(k+1)=isy(k)+(-isy(k)+(1/re)*usy(k+1)+(rrk*(kr^2)/(re*lm))*psiry(k)-(kr/re)*w(k)*psirx(k)-(kr*lbe/re)*wk(k)*isx(k))*dt/Te1;

% Поток psirx (В)

psirx(k+1)=psirx(k)+(-psirx(k)+lm*isx(k)+(lm/(rrk*kr))*(wk(k)-w(k))*psiry(k))*dt/Tr1;

% Поток psiry (Г)

psiry(k+1)=psiry(k)+(-psiry(k)+lm*isy(k)-(lm/(rrk*kr))*(wk(k)-w(k))*psirx(k))*dt/Tr1;

% Электромагнитный момент (Д)

m(k+1)=ZetaN*kr*(psirx(k+1)*isy(k+1)-psiry(k+1)*isx(k+1));

% Механическая скорость (Е)

wm(k+1)=wm(k)+(m(k+1)-mc)*dt/Tj;

% Электрическая скорость (Ж)

w(k+1)=wm(k+1)*zp;

% Моделирование наблюдателя (номер 8)

% Модуль потокосцепления ротора

psirx_oc(k+1)=psirx_oc(k)+(-psirx_oc(k)+lm*isx(k+1))*dt/Tr1;

% Скольжение

beta_psir(k+1)=isy(k+1)*rrk*kr/psirx_oc(k+1);

% Угловая скорость вращения системы координат

wk(k+1)=beta_psir(k+1)+w(k+1);

% mass

mass_t(k)=k*dt;

mass_psirx_oc(k)=psirx_oc(k+1);

mass_psiry(k)=psiry(k+1);

mass_m(k)=m(k+1);

mass_w(k)=w(k+1);

end;

% Построениеграфиков

figure(1);

plot(mass_t,mass_w,'b');

grid on;

figure(2);

plot(mass_t,mass_m,'b');

grid on;

figure(3);

plot(mass_t,mass_psirx_oc,'b',mass_t,mass_psiry,'r');

grid on;

Числовые значения параметров выводятся в окне Workspace (рис. 17).

Рис. 17. Числовые значения параметров в окне Workspace

Функциональная схема модели САР скорости асинхронного двигателя в Matlab-Script приведена на рис. 18.

Рис. 18. Функциональная схема модели САР скорости асинхронного двигателя в Matlab-Script

Результаты моделирования САР скорости асинхронного двигателя в Matlab-Script даны на рис. 19.

Рис. 19. Графики скорости, электромагнитного момента и потоков

Литература:

1. Емельянов А.А., Гусев В.М., Пестеров Д.И., Даниленко Д.С., Бесклеткин В.В., Быстрых Д.А., Иванин А.Ю. Моделирование САР скорости асинхронного двигателя с переменными ψr - is с контуром потока в системе относительных единиц // Молодой ученый. - 2018. - №11. - С. 14-32.
2. Шрейнер Р.Т. Математическое моделирование электроприводов переменного тока с полупроводниковыми преобразователями частоты. - Екатеринбург: УРО РАН, 2000. - 654 с.
3. Шрейнер Р.Т. Электромеханические и тепловые режимы асинхронных двигателей в системах частотного управления: учеб. пособие / Р.Т. Шрейнер, А.В. Костылев, В.К. Кривовяз, С.И. Шилин. Под ред. проф. д.т.н. Р.Т. Шрейнера. - Екатеринбург: ГОУ ВПО «Рос. гос. проф.-пед. ун-т», 2008. - 361 с.