Моделирование САР скорости системы «АИН ШИМ – АД» с переменными ψr - is в Matlab-Script в системе относительных единиц



В работе [1] приведена модель САР скорости системы «АИН ШИМ – АД» в Simulink. В этой статье покажем поэтапное преобразование всех элементов САР скорости «АИН ШИМ – АД» в Matlab-Script. На рис. 1 приводим всю систему, в которой даны модель асинхронного двигателя (номер 13), автономный инвертор напряжения с широтно-импульсной модуляцией (АИН ШИМ) (номер 10), генератор пилообразного напряжения (ГПН) (номер 9), преобразователи координат под номерами 7, 8, 11, 12, 15 и 16. В контурах тока по проекциям x и y соответствующие ПИ-регуляторы тока (номера 4 и 6), в контуре скорости П-регулятор скорости (номер 1).

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

Алгоритм перевода всех элементов САР скорости системы «АИН ШИМ – АД»:

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

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

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

Рис. 1. Математическая модель САР скорости системы «АИН ШИМ – АД»

Математическая модель асинхронного двигателя с переменными 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;

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;

Математическое моделирование АИН ШИМ

Математические модели АИН ШИМ (номер 10) и генератора пилообразного напряжения ГПН (номер 9) в Simulink даны на рис. 9 и 10. Работа АИН ШИМ была рассмотрена нами в статьях за 2016 г.

Рис. 9. Генератор пилообразного напряжения (ГПН) в Simulink

Рис. 10. Математическая модель АИН ШИМ в Simulink

Математическая модель АИН ШИМ с генератором пилообразного напряжения в Matlab-Script представлена в листинге 2.

Листинг 2

% Моделирование ГПН (номер 9)

U0=1;

uop(1)=1;

tau(1)=0;

f_op=1000;

tau(k+1)=tau(k)+dt*f_op;

if tau(k+1)>=1

tau(k+1)=tau(k+1)-1;

end

if (tau(k+1)>=0) && (tau(k+1)<0.5)

f(k)=1-4*tau(k+1);

else

f(k)=4*tau(k+1)-3;

end

uop(k+1)=U0*f(k);

% Моделирование АИН ШИМ (номер 10)

% Выходные сигналы нуль-органов

if usa(k+1)>=uop(k+1)

fa(k+1)=0.9;

else

fa(k+1)=-0.9;

end

if usb(k+1)>=uop(k+1)

fb(k+1)=0.9;

else

fb(k+1)=-0.9;

end

if usc(k+1)>=uop(k+1)

fc(k+1)=0.9;

else

fc(k+1)=-0.9;

end

% Импульсные напряжения на выходе АИН ШИМ

up=2.2;

usa_pwm(k+1)=up*(1/6)*(2*fa(k+1)-fb(k+1)-fc(k+1));

usb_pwm(k+1)=up*(1/6)*(-fa(k+1)+2*fb(k+1)-fc(k+1));

usc_pwm(k+1)=up*(1/6)*(-fa(k+1)-fb(k+1)+2*fc(k+1));

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

Преобразователи координат u_sx, u_sy → u_sα, u_sβ (номер 7) и u_sα, u_sβ → u_sa, u_sb, u_sc (номер 8) в Simulink приведены на рис. 11 и 12 [2].

Рис. 11. Преобразователь координат u_sx, u_sy → u_sα, u_sβ

Рис. 12. Преобразователь координат u_sα, u_sβ → u_sa, u_sb, u_sc

Реализация преобразователей координат под номерами 7 и 8 в Matlab-Script представлена в листинге 3.

Листинг 3

% Преобразователь координат usx,usy -> us_alpha,us_beta (номер 7)

teta_psir(1)=0;

rox=cos(teta_psir(k));

roy=sin(teta_psir(k));

us_alpha(k+1)=rox*usx(k+1)-roy*usy(k+1);

us_beta(k+1)=roy*usx(k+1)+rox*usy(k+1);

% Преобразователькоординат us_alpha,us_beta -> usa,usb,usc (номер 8)

usa(k+1)=us_alpha(k+1);

usb(k+1)=-(1/2)*us_alpha(k+1)+us_beta(k+1)*sqrt(3)/2;

usc(k+1)=-(1/2)*us_alpha(k+1)-us_beta(k+1)*sqrt(3)/2;

Преобразователи координат u_{а шим}, u_{b шим}, u_{c шим} → u_sα, u_sβ (номер 11) и u_sα, u_sβ → u_sx, u_sy (номер 12) в Simulink даны на рис. 13 и 14.

Рис. 13. Преобразователь координат u_{а шим}, u_{b шим}, u_{c шим} → u_sα, u_sβ

Рис. 14. Преобразователь координат u_sα, u_sβ → u_sx, u_sy

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

Листинг 4

% Преобразователь координат usa,usb,usc -> us_alpha,us_beta (номер 11)

us_alpha1(k+1)=(1/3)*(2*usa_pwm(k+1)-usb_pwm(k+1)-usc_pwm(k+1));

us_beta1(k+1)=(1/sqrt(3))*(usb_pwm(k+1)-usc_pwm(k+1));

% Преобразователькоординат us_alpha,us_beta -> usx,usy (номер 12)

usx1(k+1)=rox*us_alpha1(k+1)+roy*us_beta1(k+1);

usy1(k+1)=-roy*us_alpha1(k+1)+rox*us_beta1(k+1);

Обратные преобразователи координат по статорным токам с номерами 15 и 16 в Simulink приведены на рис. 15 и 16 [2].

Рис. 15. Преобразователь координат i_sx, i_sy → i_sα, i_sβ

Рис. 16. Преобразователь координат i_sα, i_sβ → i_sa, i_sb, i_sc

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

Листинг 5

% Преобразователь координат isx,isy -> is_alpha,is_beta (номер 15)

is_alpha(k+1)=rox*isx(k+1)-roy*isy(k+1);

is_beta(k+1)=roy*isx(k+1)+rox*isy(k+1);

% Преобразователькоординат is_alpha,is_beta -> isa,isb,isc (номер 16)

isa(k+1)=is_alpha(k+1);

isb(k+1)=-(1/2)*is_alpha(k+1)+is_beta(k+1)*sqrt(3)/2;

isc(k+1)=-(1/2)*is_alpha(k+1)-is_beta(k+1)*sqrt(3)/2;

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

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

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

Обозначим:

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

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

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

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

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

где

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

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

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

Уравнение напряжения задания будет иметь следующий вид:

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

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

где

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

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

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

Листинг 6

Tm=0.0025;dt=0.000001;

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

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

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

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

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

%Пропорциональная часть задания 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);

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

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

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

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

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

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

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

Угол потока ротора (системы координат):

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

Листинг 7

dt=0.000001; psirx_oc(1)=0.001;

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

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

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

%Угол потока ротора (системы координат)

teta_psir(k+1)=teta_psir(k)+wk(k+1)*dt;

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

Модель ПИ-регулятора потока в Simulink (номер 2) дана на рис. 20.

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

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

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

гдеn = 2.

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

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

где

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

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

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

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

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

Листинг 8

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;

% Моделирование регулятора потока (номер 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) в Simulink [1] дана на рис. 21.

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

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

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

где

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

Листинг 9

Tm=0.0025;

w(1)=0;

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

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

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

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

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

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

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

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

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

Листинг 10

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

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

% Звено компенсации 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 (рис. 23).

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

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

Листинг 11

tn=0.2;

tk=0.4;

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 дана в листинге 12.

Листинг 12

dt=0.000001;

Tm1=0.0075;

wzad1(1)=0;

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

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

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

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

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

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

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

Листинг 13

psirx_oc(1)=0.001;

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

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

Моделирование САР скорости системы «АИН ШИМ – АД»

Полная математическая модель САР скорости системы «АИН ШИМ – АД» в Matlab-Script приведена в листинге 14.

Листинг 14

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

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.2; tk=0.4;

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;

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

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

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

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

% Моделирование регулятора потока (номер 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);

% Задание 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);

% Преобразователькоординат usx,usy -> us_alpha,us_beta (номер 7)

teta_psir(1)=0;

rox=cos(teta_psir(k));

roy=sin(teta_psir(k));

us_alpha(k+1)=rox*usx(k+1)-roy*usy(k+1);

us_beta(k+1)=roy*usx(k+1)+rox*usy(k+1);

% Преобразователькоординат us_alpha,us_beta -> usa,usb,usc (номер 8)

usa(k+1)=us_alpha(k+1);

usb(k+1)=-(1/2)*us_alpha(k+1)+us_beta(k+1)*sqrt(3)/2;

usc(k+1)=-(1/2)*us_alpha(k+1)-us_beta(k+1)*sqrt(3)/2;

% Моделирование ГПН (номер 9)

U0=1; uop(1)=1; tau(1)=0; f_op=1000;

tau(k+1)=tau(k)+dt*f_op;

if tau(k+1)>=1

tau(k+1)=tau(k+1)-1;

end

if (tau(k+1)>=0) && (tau(k+1)<0.5)

f(k)=1-4*tau(k+1);

else

f(k)=4*tau(k+1)-3;

end

uop(k+1)=U0*f(k);

% Моделирование АИН ШИМ (номер 10)

% Выходные сигналы нуль-органов

if usa(k+1)>=uop(k+1)

fa(k+1)=0.9;

else

fa(k+1)=-0.9;

end

if usb(k+1)>=uop(k+1)

fb(k+1)=0.9;

else

fb(k+1)=-0.9;

end

if usc(k+1)>=uop(k+1)

fc(k+1)=0.9;

else

fc(k+1)=-0.9;

end

% Импульсные напряжения на выходе АИН ШИМ

up=2.2;

usa_pwm(k+1)=up*(1/6)*(2*fa(k+1)-fb(k+1)-fc(k+1));

usb_pwm(k+1)=up*(1/6)*(-fa(k+1)+2*fb(k+1)-fc(k+1));

usc_pwm(k+1)=up*(1/6)*(-fa(k+1)-fb(k+1)+2*fc(k+1));

% Преобразователькоординат usa,usb,usc -> us_alpha,us_beta (номер 11)

us_alpha1(k+1)=(1/3)*(2*usa_pwm(k+1)-usb_pwm(k+1)-usc_pwm(k+1));

us_beta1(k+1)=(1/sqrt(3))*(usb_pwm(k+1)-usc_pwm(k+1));

% Преобразователькоординат us_alpha,us_beta -> usx,usy (номер 12)

usx1(k+1)=rox*us_alpha1(k+1)+roy*us_beta1(k+1);

usy1(k+1)=-roy*us_alpha1(k+1)+rox*us_beta1(k+1);

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

% Ток isx (А)

isx(k+1)=isx(k)+(-isx(k)+(1/re)*usx1(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)*usy1(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;

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

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

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

% Угол потока ротора (системы координат)

teta_psir(k+1)=teta_psir(k)+wk(k+1)*dt;

% Преобразователь координат isx,isy -> is_alpha,is_beta (номер 15)

is_alpha(k+1)=rox*isx(k+1)-roy*isy(k+1);

is_beta(k+1)=roy*isx(k+1)+rox*isy(k+1);

% Преобразователькоординат is_alpha,is_beta -> isa,isb,isc (номер 16)

isa(k+1)=is_alpha(k+1);

isb(k+1)=-(1/2)*is_alpha(k+1)+is_beta(k+1)*sqrt(3)/2;

isc(k+1)=-(1/2)*is_alpha(k+1)-is_beta(k+1)*sqrt(3)/2;

% 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 (рис. 25).

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

Функциональная схема модели САР скорости системы «АИН ШИМ – АД» в Matlab-Script приведена на рис. 26.

Рис. 26. Функциональная схема модели САР скорости системы «АИН ШИМ – АД» в Matlab-Script

Результаты моделирования САР скорости системы «АИН ШИМ – АД» в Matlab-Script даны на рис. 27.

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

Литература:

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