В работе [1] приведена модель САР скорости асинхронного двигателя в Simulink в системе абсолютных единиц. В этой статье покажем поэтапное преобразование всех элементов САР скорости в Matlab-Script. На рис. 1 приводим всю систему, в которой даны модель асинхронного двигателя (номер 7), в контурах тока по проекциям x и y соответствующие ПИ-регуляторы тока (номера 4 и 6), в контуре скорости П-регулятор скорости (номер 1).
Рис. 1. Математическая модель САР скорости асинхронного двигателя в Simulink в системе абсолютных единиц
Важным элементом является контур потока с ПИ-регулятором потока (номер 2). Для ориентации системы координат по потокосцеплению ротора вводится наблюдатель (номер 8). В модели учтена компенсация перекрестных связей (номер 5). Сигнал задания по скорости выполнен на задатчике интенсивности. В цепи задания скорости перед регулятором скорости предусмотрен фильтр.
Алгоритм перевода всех элементов САР скорости:
‒ приводится математическая формула той или иной переменной, выраженной в Simulink;
‒ приводится его структурная схема;
‒ переход от изображений к оригиналу (от s к d/dt) и решение с помощью простого метода Эйлера.
Математическая модель асинхронного двигателя с переменными IS – ΨR
А) Выражение потокосцепления ΨRx по проекции x из работы [1] имеет вид:
|
|
(1) |
где
- электрическая скорость вращения ротора;
- механическая угловая скорость на валу двигателя.
Структурная схема для определения ΨRx в Simulink приведена на рис. 2.
Рис. 2. Структурная схема для определения потокосцепления ΨRx в Simulink
Преобразуем уравнение (1) для программирования в Matlab-Script:
Обозначим 
, тогда:
Переходим к оригиналу 
:
Переходим к конечным разностям (метод Эйлера):
Отсюда потокосцепление ΨRx в Matlab-Script определится следующим образом:
Б) Уравнение для определения тока ISx в Simulink, полученное в работе [1], имеет следующий вид:
|
|
(2) |
Структурная схема для определения ISx в Simulink приведена на рис. 3.
Рис. 3. Структурная схема для определения тока ISx в Simulink
Преобразуем выражение тока ISx для программирования в Matlab-Script:
Обозначим 
, тогда:
Переходим к оригиналу:
Переходим к конечным разностям:
Ток ISx в Matlab-Script определится следующим образом:
В) Уравнение для определения потокосцепления ΨRyв Simulink имеет вид:
|
|
(3) |
Структурная схема для определения ΨRy в Simulink представлена на рис. 4.
Рис. 4. Структурная схема для определения потокосцепления ΨRy в Simulink
Преобразуем уравнение (3) для программирования в Matlab-Script:
Г) Выражение тока ISy в Simulink имеет следующий вид [1]:
|
|
(4) |
Структурная схема реализации уравнения (4) дана на рис. 5.
Рис. 5. Структурная схема для определения тока ISy в Simulink
Отсюда ток ISy в Matlab-Script определится следующим образом:
Д) На рис. 6 представлена структурная схема для реализации уравнения электромагнитного момента в Simulink:
Рис. 6. Математическая модель определения электромагнитного момента M в Simulink
Уравнение электромагнитного момента для реализации в Matlab-Script:
Е) Механическая угловая скорость вращения вала двигателя в Simulink (рис. 7):
Рис. 7. Математическая модель определения механической угловой скорости вращения вала двигателя в Simulink
Отсюда механическая угловая скорость вращения вала двигателя в Matlab-Script:
Ж) Электрическая скорость вращения ротора в Simulink (рис. 8):
Рис. 8. Математическая модель определения электрической скорости вращения ротора в Simulink
Электрическая скорость вращения ротора в Matlab-Script:
Реализация математической модели асинхронного двигателя с короткозамкнутым ротором с переменными ΨR - IS в системе абсолютных единиц в Matlab-Script представлена в листинге 1. Параметры асинхронного двигателя рассмотрены в работах [2] и [3].
Листинг 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;
Psib=Ub/Omegab;
Lb=Psib/Ib;
kd=1.0084;
Mb=kd*PN/OmegaN;
Pb=Mb*Omegarb;
% Расчет коэффициентов АД
rs=Rs/Zb;
lbs=Xs/Zb;
lbr=Xr/Zb;
lm=Xm/Zb;
Lm=lm*Lb;
betaN=(Omega0N-OmegaN)/Omega0N;
kr=lm/(lm+lbr);
lbe=lbs+lbr+lbs*lbr*lm^(-1);
Lbe=lbe*Lb;
roN=0.9962;
rrk=roN*betaN;
Rrk=rrk*Zb;
Tr=lm/(rrk*kr);
Tr1=Tr/Omegab;
re=rs+rrk*kr^2;
Re=re*Zb;
Te=kr*lbe/re;
Te1=Te/Omegab;
% Расчет асинхронного двигателя (номер 7)
K=input('Длительность цикла k=');
for k=1:K
dt=0.00001;
Usx(k)=0; Usy(k)=Ub; Omegak=314;
Isx(1)=0; Isy(1)=0; Psirx(1)=0; Psiry(1)=0;
Omegam(1)=0; Omega(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)*Omega(k)*Psiry(k)+(kr*Lbe/Re)*Omegak(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)*Omega(k)*Psirx(k)-(kr*Lbe/Re)*Omegak(k)*Isx(k))*dt/Te1;
% Поток Psirx (В)
Psirx(k+1)=Psirx(k)+(-Psirx(k)+Lm*Isx(k)+(Lm/(Rrk*kr))*(Omegak(k)-Omega(k))*Psiry(k))*dt/Tr1;
% Поток Psiry (Г)
Psiry(k+1)=Psiry(k)+(-Psiry(k)+Lm*Isy(k)-(Lm/(Rrk*kr))*(Omegak(k)-Omega(k))*Psirx(k))*dt/Tr1;
% Электромагнитныймомент (Д)
M(k+1)=(3/2)*zp*kr*(Psirx(k+1)*Isy(k+1)-Psiry(k+1)*Isx(k+1));
% Механическая скорость (Е)
Omegam(k+1)=Omegam(k)+(M(k+1)-Mc)*dt/J;
% Электрическая скорость (Ж)
Omega(k+1)=Omegam(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.00001;
KI=Te1/(2*Tm/Re); TI=2*Tm/Re;
Isx(1)=0; Isy(1)=0;Ux2(1)=0;Uy2(1)=0;
% Моделирование регуляторов тока (номера 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);
Математическое моделирование наблюдателя потокосцепления ротора
Модель наблюдателя потокосцепления ротора в Simulink, полученная в работе [1], приведена на рис. 11. Преобразуем эту модель в Matlab-Script.
Рис. 11. Модель наблюдателя потокосцепления ротора в Simulink
Приведем уравнение модуля потокосцепления ротора к оригиналу:
Переходим к конечным разностям:
Уравнение скольжения 
для программирования в Matlab-Script будет иметь вид [1], [2], [3]:
Отсюда угловая скорость вращения системы координат 
:
Математическая модель наблюдателя в Matlab-Script приведена в листинге 3.
Листинг 3
dt=0.00001; 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);
% Угловая скорость вращения системы координат
Omegak(k+1)=Beta_Psir(k+1)+Omega(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.00001;
Psirx_oc(1)=0.001; Ixzad2(1)=0;
KPsi=Tr1/(4*n*Tm*Lm);
TPsi=4*n*Tm*Lm;
% Моделирование регулятора потока (номер 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);
Математическое моделирование регулятора скорости
Математическая модель П-регулятора скорости в Simulink [1] дана на рис. 13.
Рис. 13. Пропорциональный регулятор скорости в Simulink
Передаточная функция регулятора скорости:
Отсюда определим задание момента 
:
где
Математическая модель регулятора скорости в Matlab-Script представлена в листинге 5.
Листинг 5
Tm=0.0025; Omega(1)=0;
% Моделирование регулятора скорости (номер 1)
Omega_sum(k+1)=Omega_zad1(k+1)-Omega(k);
% Задание момента M
Mzad(k+1)=Omega_sum(k+1)*J/(4*Tm);
Математическое моделирование компенсации перекрестных связей
Математическая модель компенсации перекрестных связей в Simulink [1] дана на рис. 14.
Рис. 14. Компенсация внутренних перекрестных связей в Simulink
Компенсационные составляющие каналов управления определятся следующим образом:
Реализация математической модели компенсации перекрестных связей в Matlab-Script представлена в листинге 6.
Листинг 6
Isx(1)=0; Isy(1)=0; Psirx_oc(1)=0.001; Omegak(1)=0;
% Моделирование звена компенсации (номер 5)
% Звено компенсации x
Ukx(k+1)=-Omegak(k)*kr*Lbe*Isy(k);
% Звено компенсации y
Uky(k+1)=Omegak(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);
Математическое моделирование задатчика интенсивности
Задание на скорость Ω* в Simulink формируется в блоке Signal Builder (рис. 15).
Рис. 15. Сигнал задания на скорость Ω* в Simulink
Программирование сигнала задания на скорость в Matlab-Script представлено в листинге 7.
Листинг 7
tn=0.4;
tk=1.3;
dt=0.00001;
% Задание на скорость
if((k*dt>=0)&&(k*dt<=tn))
Omega_zad(k+1)=0;
end;
if((k*dt>=tn)&&(k*dt<=tk))
Omega_zad(k+1)=314*(k*dt-tn)/(tk-tn);
end;
if(k*dt>tk)
Omega_zad(k+1)=314;
end;
Математическое моделирование задания по скорости на выходе фильтра
Передаточная функция фильтра:
Определим задание скорости 
на выходе фильтра:
Перейдем от изображения к оригиналу:
Переходим к конечным разностям:
Математическая модель задания скорости на выходе фильтра в Matlab-Script дана в листинге 8.
Листинг 8
dt=0.00001;
Tm1=0.0075;
Omega_zad1(1)=0;
% Задание скорости на выходе фильтра
Omega_zad1(k+1)=Omega_zad1(k)+(Omega_zad(k+1)-Omega_zad1(k))*dt/Tm1;
Математическое моделирование задания статорного тока по проекции y
Математическая модель задания тока 
в Simulink дана на рис. 16.
Рис. 16. Реализация задания статорного тока 
в Simulink
Задание на статорный ток по проекции y:
Математическая модель задания 
в Matlab-Script представлена в листинге 9.
Листинг 9
Psirx_oc(1)=0.001;
% Задание Isy (номер 3)
Iyzad(k+1)=Mzad(k+1)/(Psirx_oc(k)*(3/2)*zp*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;
Psib=Ub/Omegab;
Lb=Psib/Ib;
kd=1.0084;
Mb=kd*PN/OmegaN;
Pb=Mb*Omegarb;
% Расчет коэффициентов АД
rs=Rs/Zb;
lbs=Xs/Zb;
lbr=Xr/Zb;
lm=Xm/Zb;
Lm=lm*Lb;
betaN=(Omega0N-OmegaN)/Omega0N;
kr=lm/(lm+lbr);
lbe=lbs+lbr+lbs*lbr*lm^(-1);
Lbe=lbe*Lb;
roN=0.9962;
rrk=roN*betaN;
Rrk=rrk*Zb;
Tr=lm/(rrk*kr);
Tr1=Tr/Omegab;
re=rs+rrk*kr^2;
Re=re*Zb;
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=1.612;
tn=0.4;
tk=1.3;
dt=0.00001;
% Расчет САР скорости АД
K=input('Длительность цикла k=');
for k=1:K
% Параметры САР скорости в начальный момент времени
Isx(1)=0; Isy(1)=0; Psirx(1)=0; Psiry(1)=0;
Omegam(1)=0; Omega(1)=0; Mc=0;
Psirx_oc(1)=0.001;
Ixzad2(1)=0;
Ux2(1)=0;
Uy2(1)=0;
Omegak(1)=0;
Omega_zad1(1)=0;
% Задание на скорость
if((k*dt>=0)&&(k*dt<=tn))
Omega_zad(k+1)=0;
end;
if((k*dt>=tn)&&(k*dt<=tk))
Omega_zad(k+1)=314*(k*dt-tn)/(tk-tn);
end;
if(k*dt>tk)
Omega_zad(k+1)=314;
end;
% Задание скорости на выходе фильтра
Omega_zad1(k+1)=Omega_zad1(k)+(Omega_zad(k+1)-Omega_zad1(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)
Omega_sum(k+1)=Omega_zad1(k+1)-Omega(k);
% Задание момента M
Mzad(k+1)=Omega_sum(k+1)*J/(4*Tm);
% Задание Isy (номер 3)
Iyzad(k+1)=Mzad(k+1)/(Psirx_oc(k)*(3/2)*zp*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)=-Omegak(k)*kr*Lbe*Isy(k);
% Звено компенсации y
Uky(k+1)=Omegak(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)*Omega(k)*Psiry(k)+(kr*Lbe/Re)*Omegak(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)*Omega(k)*Psirx(k)-(kr*Lbe/Re)*Omegak(k)*Isx(k))*dt/Te1;
% Поток Psirx (В)
Psirx(k+1)=Psirx(k)+(-Psirx(k)+Lm*Isx(k)+(Lm/(Rrk*kr))*(Omegak(k)-Omega(k))*Psiry(k))*dt/Tr1;
% Поток Psiry (Г)
Psiry(k+1)=Psiry(k)+(-Psiry(k)+Lm*Isy(k)-(Lm/(Rrk*kr))*(Omegak(k)-Omega(k))*Psirx(k))*dt/Tr1;
% Электромагнитный момент (Д)
M(k+1)=(3/2)*zp*kr*(Psirx(k+1)*Isy(k+1)-Psiry(k+1)*Isx(k+1));
% Механическая скорость (Е)
Omegam(k+1)=Omegam(k)+(M(k+1)-Mc)*dt/J;
% Электрическая скорость (Ж)
Omega(k+1)=Omegam(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);
% Угловая скорость вращения системы координат
Omegak(k+1)=Beta_Psir(k+1)+Omega(k+1);
% mass
mass_t(k)=k*dt;
mass_M(k)=M(k+1);
mass_Omega(k)=Omega(k+1);
mass_Psirx_oc(k)=Psirx_oc(k+1);
mass_Psiry(k)=Psiry(k+1);
end;
% Построениеграфиков
figure(1);
plot(mass_t,mass_Omega,'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. Результаты моделирования САР скорости асинхронного двигателя даны на рис. 19.
Рис. 18. Функциональная схема модели САР скорости асинхронного двигателя в системе абсолютных единиц в Matlab-Script
Рис. 19. Графики скорости, электромагнитного момента и потоков
Литература:
- Емельянов А.А., Гусев В.М., Пестеров Д.И., Даниленко Д.С., Бесклеткин В.В., Быстрых Д.А., Иванин А.Ю. Моделирование САР скорости асинхронного двигателя с переменными ΨR - IS с контуром потока в системе абсолютных единиц // Молодой ученый. - 2018. - №13. - С. 22-40.
- Шрейнер Р.Т. Математическое моделирование электроприводов переменного тока с полупроводниковыми преобразователями частоты. - Екатеринбург: УРО РАН, 2000. - 654 с.
- Шрейнер Р.Т. Электромеханические и тепловые режимы асинхронных двигателей в системах частотного управления: учеб. пособие / Р.Т. Шрейнер, А.В. Костылев, В.К. Кривовяз, С.И. Шилин. Под ред. проф. д.т.н. Р.Т. Шрейнера. - Екатеринбург: ГОУ ВПО «Рос. гос. проф.-пед. ун-т», 2008. - 361 с.
