script
 %     Расчет переходных процессов реостатного асинхронного электропривода. 
 %                 Реостатный пуск при питании от цеховой сети, 
%                 торможение динамическое или противовключением
                          %при U1=U1н, f1=f1н, R2доб=var,
 prog='read';       
 %                  Перед расчетами ПОВТОРИТЬ разделы:
 %          "Электромеханические свойства и характеристики АД"
 %                  "Переходные процессы электропривода 
%             с прямолинейной механической характеристикой" 

%                              РАССЧИТАТЬ:
%              - правильную пусковую диаграмму;  
%                добавочные сопротивления на пусковых ступенях; 
%                скорости двигателя для переключения сопротивлений;
%             -  добавочное сопротивление на ступени торможения.
%             
%                               ВНИМАНИЕ!
%      Предложенная программа является  образцом (примером) для Ваших расчетов.
%      Перед выполнением расчетов следует сохранить данную отладочную программу 
%                    под Вашим именем, а в этой программе
%                          НИЧЕГО НЕ ИЗМЕНЯТЬ!!!!!!!!!!!
%
%                          РАБОТА С ПРОГРАММОЙ
%      1. Внимательно прочитайте свою программу, уясните алгоритм расчета.
%          Обозначения физических величин, примененные в программе, часто не 
%          совпадают с общепринятыми. Следует разобраться с их обозначениями.
%       2. Дифференциальные и алгебраические уравнения программы описаны
%          в учебном пособии к курсовому проектированию по "Теории электропривода".
%      3. Ввести:
%      - каталожные и обмоточные данные двигателя; 
%      - добавочные сопротивления в цепи якоря при пуске и торможении;
%      - скорости переключения ступеней сопротивлений (при  меньшем числе 
%        ступеней на отсутствующей очередной ступени вводится Rd=0);
%      - времена работы на участках;
%      - расчетные данные рабочего органа РО. 
%      4. Ввести начальные условия для пуска и торможения:
%      - шаг  интегрирования h (можно не менять); 
%      - учет кривой намагничивания;
%      - предполагаемые времена пуска tрр с выходом на установившийся режим, 
%       торможения противовключение ttр или динамического торможения tdt с 
%       достижением остановки двигателя. 

%     5. ПРИ ПЕРВОМ ВКЛЮЧЕНИИ ПРОГРАММЫ СЛЕДУЕТ ВЫБРАТЬ 
%                 ПЕРЕХОДНЫЙ ПРОЦЕСС ПУСКА nu=0%                            
%     Первый пуск выполняется с контролем переключения ступеней сопротивлений
%                      по скорости per = 1.
%        Нажать знак "Save and Run".%     
%     Скорости переключения определены при расчете правильной пусковой диаграммы. 
%     После наладки программы при пуске под нагрузкой с контролем переключения
%     ступеней сопротивлений по скорости определяются времена работы
%     электропривода на ступенях. На эти времена будут настроены реле времени. 
%     Все остальные пуски следует выполнять с контролем переключения ступеней 
%     сопротивлений в функции времени  per = 2.
%     По окончании   пуска программа выводит нагрузочные диаграммы c интегральными
%     показателями пуска (figure 1) и динамические  характеристики (figure 2).
%      6. Не убирая с экрана figure1 и figure2, перейти к расчету торможения.     
%     установив (nu=1 или nu=2 в зависимости от типа торможения) и вновь нажать 
%                         знак "Save and Run".
%      По окончании торможения figure1 и figure2 будут дополнены данными 
%      переходного процесса торможения. 
%      7. Сохраните figure1 и figure2 для печати и дальнейшей обработки
%      результатов расчетов (проверок по нагреву, по перегрузочной способности, 
%      энергетических показателей электропривода и др.)
    
%                         ДВИГАТЕЛЬ

%Каталожные данные двигателя (при ПВ=40%):
tip='4MTH225M8';%Тип двигателя
ppo=4;           %Число пар полюсов
u1n=220;        %Номинальное фазное напряжение,В
f1n=50;         %Номинальная частота,Гц
Pn=30000;       %Номинальная мощность, Вт
nn=720;         %Номинальная скорость,об/мин
I1n=74.6;       %Номинальный ток статора,А 
kpdn=0.84;      %Номинальный кпд
cosfi1n=0.72;   %Номинальный cos fi
Mk=1030;        %Критический момент, Нм
I2n=71;         %Номинальный ток ротора, А
E20=252;        %Напряжение на кольцах неподвижного разомкнутого ротора, В
Jd=1.07;        %Момент инерции двигателя,кгм2

%Обмоточные данные двигателя
r1=0.137;        %Невыключаемое активное сопротивление статора,Ом
r2=0.08;         %Невыключаемое активное сопротивление ротора,Ом
x1=0.23;         %Невыключаемое индуктивное сопротивление статора,Ом
x2=0.254;        %Невыключаемое индуктивное сопротивление ротора,Ом    

 % Добавочные ступени сопротивления в цепи ротора
 % Если ступень отсутствует, Rd = 0           
   Rd1=0.41;    %  1 ступень
   Rd2=0.205;   %  2 ступень
   Rd3=0.1;     %  3 ступень
   Rd4=0.05;    %  4 ступень
   Rd5=2.05;    %  ступень при торможении противовключением
   Rd6=1.025;   %  ступень при динамическом торможении

   %Контроль переключения ступеней по скорости 
   w1=36;    %Участок 1 - скорость переключения на ступень 2 0.45*won
   w2=56;    %Участок 2 - скорость переключения на ступень 3 0.68*won
   w3=66;    %Участок 3 - скорость переключения на ступень 4 0.8*won
   w4=72;    %Участок 4 - скорость переключения  на ест.характеристику 0.9*won

   %Контроль переключения ступеней по времени 
   t1=0.19;       %Участок 1 - время разгона на ступени 1, с
   t2=0.28;       %Участок 2 - время разгона на ступенях 1-2, с
   t3=0.34;       %Участок 3 - время разгона на ступенях 1-2-3, с
   t4=0.37;       %Участок 4 - время разгона на ступенях 1-2-3-4, с
   %t5=0.5;        %Участок 5 - время дотягивания и начало установившегося
              % режима работы на естественной характеристике, с
              
   %Динамическое торможение   
    Ipost=1.5*I1n;  %Постоянный ток статора в заданной точке, А
    kc=0.816;       %Коэффициент схемы включения постоянного тока
    i1s=kc*Ipost    %Эквивалентный переменный ток фазы статора 
                    %в заданной точке, А 
    krn=1;      %При учете кривой намагничивания krn=1
  
   Mc1=300;   % Статический момент при пуске,Нм
   Mc2=270;   % Статический момент при торможении 
   Mca=0;     %Статический момент активный Мса=0,реактивный Мса=1 
   Jpr=0.2;   %Приведенный момент инерции РО,кгм^2
 
   tpp=0.4;   %конец расчета пуска, с65
   tpt=0.77;    %конец расчета торможения, с

   h=0.001;    %Шаг интегрирования,с

               %  НАЧАЛЬНЫЕ УСЛОВИЯ
 % Переходный процесс:
                  %nu=0  %пуск  
                  %nu=1 %торможение противовключением 
                  %nu=2 %динамическое торможение                  
                   nu=2
 %Контроль переключения ступеней сопротивлений 
 %per=1   % по скорости
 %per=2   % по времени  
                   per=1  
                   
                   %Расчетные данные двигателя
pi=3.1416;
won=2*pi*f1n/ppo; %Номинальная синхронная скорость,рад/с
wn=nn/9.55;      %Номинальная скорость,рад/с
sn=(won-wn)/won; %Номинальнгое скольжение
Mn=Pn/wn;        %Номинальный момент на валу,Нм
Mko=Mk/Mn;       %Критический момент, o.e.
sk=sn*(Mko+sqrt(Mko^2-1));
sinfi1n=sqrt(1-cosfi1n^2);
Imn=I1n*(sinfi1n-cosfi1n*sn/sk); 
xmn=u1n/Imn-x1;%Индуктивное сопротивление контура намагничивания,Ом
E1n=Imn*xmn;
ke=1.73*E1n/E20;
r2p=r2*ke^2;
x2p=x2*ke^2;
a1=r1/r2p;
Men=2*Mk*(1+a1*sk)/(sn/sk+sk/sn+2*a1*sk); %Номинальный электромагнитный момент,Нм
dMxn=Men-Mn;      %Потери холостого хода в номинальном режиме, Нм
R2n=E20/(1.73*I2n);%Номинальное активное сопротивление ротора, Ом
J=Jd+Jpr;

if nu==0, t0=h; w=0; tk=tpp; end
if nu==1, t0=tpp; w=wk; tk=tpt; end
if nu==2, t0=tpp; tk=tpt;w=wk; end

if nu==0 
   w=0;Pv=0;Pc=0;Qc=0;I1kvt=0;I2kvt=0;al=0;;dP2=0;dPd=0;dP1=0;dP=0;
   M=0;I1a=0;r2s=0;Im=Imn;I2p=0;I1q=0;zc=0;I1=0;I2=I2n; xm=xmn;t=0; 
else 
    w=wk;Pv=0;Pc=0;Qc=0;I1kvt=0;al=0;;dP2=0;dPd=0;dP1=0;dP=0;
    I2kvt=0;M=0;I1a=0;r2s=0;Im=Imn;I2p=0;I1q=0;zc=0;I1=0;I2=I2n; xm=xmn;t=tpp; 
end

for t=t0:h:tk
    t
    if and (nu==0, per==1), Rd=Rd1+Rd2+Rd3+Rd4; 
       if w>w1  Rd=Rd2+Rd3+Rd4; end 
       if w>w2  Rd=Rd3+Rd4; end
       if w>w3  Rd=Rd4; end
       if w>w4  Rd=0;   end 
    end  
    if and (nu==0, per==2), Rd=Rd1+Rd2+Rd3+Rd4;      
        if t>t1  Rd=Rd2+Rd3+Rd4; end 
        if t>t2  Rd=Rd3+Rd4; end
        if t>t3  Rd=Rd4; end
        if t>t4  Rd=0;   end 
    end    
    
   if nu==1, won=-abs(won); Rd=Rd5; end    %Торможение противовключением   
   if nu==2  s=-w/won; Rd=Rd6;             %Динамическое торможение 
   else s=(won-w)/won;
   end

   if abs(s)<0.001, s=0.001;end 

   if krn==1,
      imo=Im/Imn; 
      if imo<0.37 em=1.35*imo; end                  %kk=1.35;
      if imo>0.37 em=0.5+1.09*(imo-0.37); end       %kk=1.095; 
      if imo>0.6 em=0.75+0.625*(imo-0.6); end       %kk=0.625;
      if imo>1 em=1+0.5*(imo-1); end                %kk=0.5;
      if imo>1.2 em=1.1+0.167*(imo-1.2);end         %kk=0.167; 
      if imo>1.8 em=1.2+0.1*(imo-1.8); end          %kk=0.1;
      if em>1.3 em=1.3;end 
      xm=em*E1n/(imo*Imn);
   else 
      xm=xmn; 
   end
   Rdp=Rd*ke^2;
   z1=r1+j*x1;
   r2s=(r2p+Rdp)/s;
   x2d=x2p;
   z2p=r2s+j*x2d;
   zm=j*xm;
   zmr=zm*z2p/(zm+z2p);

   if nu==2 
      im=i1s*z2p/(z2p+zm);
      Im=abs(im);
      Te=(x2p+xm)/(314*(Rdp+r2p));
      i2p=i1s*zm/(z2p+zm);
      I2pu=abs(i2p);
      I2p=I2pu*(1-exp(-t/Te));
      i2=i2p*ke;
      I2=I2p*ke;     
   else  
      zc=z1+zmr;
      i1=u1n/zc; 
      fi1=angle(zc);
      I1a=real(i1);
      I1=abs(i1);  
      I1a=I1*cos(fi1);
      I1q=I1*sin(fi1);
      e=i1*zmr;
      i2p=e/z2p;
      i2=i2p*ke;
      I2p=abs(i2p);
      I2=I2p*ke;
      im=e/zm; 
      Im=abs(im);  
    end  
    if I2==0, I2p=0; M=0; I1=Im; end    
    M=3*I2p^2*r2s/won; 
    dMx=abs(dMxn)*sign(w);     %Реактивный момент потерь хх двигателя,Нм
      if or (nu==1, nu==2), Mc=Mc2;
    else        Mc=Mc1;   end    
    if Mca==0, 
      Mc=Mc+dMx;                      %Активный Мс,Нм
   else 
      Mc=abs(Mc)*sign(w)+dMx;         %Реактивный Мс,Нм
  end       
    w=w+h*(M-Mc)/J;                   %Скорость двигателя,рад/с 
    al=al+w*h;                        %Угол поворота вала двигателя,рад
    Pv1=(M-dMx)*w;                    %Механическая мощность на валу,Вт*с
    Pv=Pv+(M-dMx)*w*h;                %Механическая энергия на валу,Вт*с
    if nu==3 Pc=0; I1=Ipost;dP1=0;    %Потери энергии в статоре при ДТ не учитываются
    else
       Pc1=3*u1n*I1a;                 %Активная мощность из сети,Вт
       Pc=Pc+Pc1*h;                   %Активная энергия из сети,Вт*с
       Qc=Qc+3*u1n*I1q*h;             %Реактивная энергия из сети,Вар*с               
       cosfi=abs(cos(angle(zc))); 
       dP1=dP1+3*I1*I1*r1*h;          %Потери энергии в обмотках статора,Вт*с 
    end   
    I1kvt=I1kvt+I1^2*h;               %Квадрат тока статора*время,А^2*с     
    I2kvt=I2kvt+I2^2*h;               %Квадрат тока ротора*время,А^2*с
    if nu==0 I1kv=sqrt(I1kvt/tpp);    %Среднеквадратичный ток статора,А
       I2kv=sqrt(I2kvt/tpp);          %Среднеквадратичный ток ротора,А
    else 
       I1kv=sqrt(I1kvt/(tpt-tpp));
       I2kv=sqrt(I2kvt/(tpt-tpp));
    end
if nu==2 I1=Ipost;I1kv=Ipost;dP1=dP1+0.67*Ipost^2*r1*h; end
    dP2=dP2+3*I2p*I2p*r2p*h;      %Потери энергии в обмотках ротора,Вт*с
    dP=dP1+dP2;                   %Потери энергии в двигателе,Вт*с
    dPd=dPd+3*I2p*I2p*Rdp*h;      %Потери энергии в добавочном сопротивлении,Вт*с
    
    if nu==3 Pc=0; kpd=0; end
    if and(Pv1>0,Pc1>0), kpd=Pv1/Pc1; end
    if and(Pv1<0,Pc1<0), kpd=Pc1/Pv1;end    
    if and(M==0,w==0)  kpd=0;end
    if kpd <=0 kpd=0; end 
    if kpd >1, kpd=0; end
    wk=w;
    
    %Печать результатов расчета(с.266-415    )
    tm=0.2;     % развертка показателей на рисунке, с

    b0=num2str(t);
    b1=num2str(Pv);
    b2=num2str(Pc);
    b3=num2str(Qc);
    b4=num2str(I1kv);
    b5=num2str(al);
    b6=num2str(dP);
    b7=num2str(dPd);
    b8=num2str(I2);
    b9=num2str(w);
    b10=num2str(M);
    b11=num2str(I1);
    b14=num2str(Mc);
    b15=num2str(I2kv);
    b18=num2str(I2n);
    b19=num2str(won);
    b20=num2str(Mn);
    b21=num2str(I1n);
    b22=num2str(t-tpp);
    b23=num2str(kpd);
    b24=num2str(cosfi);

    % Если какую-то из характеристик не нужно выводить на график,
    % поставьте знак % перед plot:   %plot()
    
    figure(1)         % - нагрузочные диаграммы w,M,I1,wo=f(t)
 
    xlabel('t')
    ylabel('w,M,I1,I2'), 
    plot(t,M/Mn,' k')
    plot(t,I1/I1n,' b')
    plot(t,I2p*ke/I2n,' r')
    plot(t,w/abs(won),' m')
    axis([0 tk+0.2 -3 5])      %рамка
    hold on 
    grid on

    figure(2)            % Динамические механические характеристики

    xlabel('M,I1,I2')
    ylabel('w')
    plot(M/Mn,w/abs(won),' k')
    plot(I1/I1n,w/abs(won),' b')
    plot(I2/I2n,w/abs(won),' r')
    %plot(kpd,w/abs(won),'g')
    %plot(cosfi,w/abs(won),'m')   
    axis([-3 5.5 -0.5 1.5])      %рамка 
    hold on 
    grid on
    end

   %динамические характеристики w=f(M,I1,I2) 
if and(figure(1),nu==0),%nu=0-пуск 
   d=date;
   text(0.25,4.6,tip,'color','k');
   text(0.05,4.2,prog,'color','k');
   text(0.05,4.6,d,'color','k'); 
   text(0.2,4.2,' Mc =','color','k');
   text(0.3,4.2,b14,'color','k');   
   text(tm+0.01,-0.2,'tp =','color','k');
   text(0.41,-0.2,b0,'color','k');    
   text(tm+0.01,-0.5,'Pv=','color','k');
   text(0.41,-0.5,b1,'color','k');
   text(tm, -0.8,'Pc =','color','k');
   text(0.41,-0.8,b2,'color','k');
   text(tm,-1.2,'Qc=','color','k');
   text(0.41,-1.2,b3,'color','k');
   text(tm-0.03,-1.5,'  I1kv =','color','k');
   text(0.41,-1.5,b4,'color','k');
   text(tm+0.01,-1.8,'al =','color','k');
   text(0.41,-1.8,b5,'color','k');
   text(tm,-2.2,'dP =','color','k');
   text(0.41,-2.2,b6,'color','k');
   text(tm-0.02,-2.5,'dPd=','color','k');
   text(0.41,-2.5,b7,'color','k');
   text(tm-0.03,-2.8,'  I2kv =','color','k');
   text(0.41,-2.8,b15,'color','k');

   text((tpp-0.1),3.4,'w =','color','m');
   text(tpp,3.4,b9,'color','m');
   text((tpp-0.1),3.0,'M =','color','k');
   text(tpp,3.0,b10,'color','k');
   text((tpp-0.1),2.6,'I1 =','color','b');
   text(tpp,2.6,b11,'color','b');
   text((tpp-0.1),2.2,'I2 =','color','r');
   text(tpp,2.2,b8,'color','r');
   text((tpp-0.13),3.8,'cosfi =','color','c');
   text(tpp,3.8,b24,'color','c');
   text((tpp-0.11),4.2,'kpd=','color','g');
   text(tpp,4.2,b23,'color','g');

  %нагрузочные диаграммы w,M,I1,I2=f(t)
   if and(figure(2),nu==0), 
      d=date;
      text(-2.8,1.3,d,'color','k');
      text(-1.2,1.3,tip,'color','k');
      text(1.6,1.3,'Mc1=','color','k');
      text(2.25,1.3,b14,'color','k');
   else
      text(1.6,1.1,'Mct=','color','k');
      text(2.25,1.1,b14,'color','k');
   end

      text(-2.8,1.05,'won=','color','k');
      text(-2.1,1.05,b19,'color','k');
      text(-2.8,0.9,'I1n =','color','b');
      text(-2.1,0.9,b21,'color','k');
      text(-2.8,0.75,'I2n =','color','r');
      text(-2.1,0.75,b18,'color','k');
      text(-2.8,0.6,'Mn =','color','b');
      text(-2.1,0.6,b20,'color','k');
      %text(-2.8,0.45,'kpd','color','g');
      %text(-2.8,0.3,'cos fi','color','m');
      hold on
      grid on
end

if or (nu==1, nu==2)  %nu=1 - торможение противовключением 
                      %nu=2 - динамическое торможение
   text(0.2,3.8,' Mct=','color','k');
   text(0.3,3.8,b14,'color','k');    
   text((tpt-0.05),-0.2,'tt =','color','k');
   text(tpt+0.01,-0.2,b22,'color','k');    
   %text(tm+0.2,-0.5,'A =','color','k');
   text(tpt+0.01,-0.5,b1,'color','k');
   %text(tm+0.15, -0.8,'Pc =','color','k');
   text(tpt+0.01,-0.8,b2,'color','k');
   %text(tm+0.15,-1.2,' Qc=','color','k');
   text(tpt+0.01,-1.2,b3,'color','k');
   %text(tm+0.1,-1.5,'I1kvt =','color','k');
   text(tpt+0.01,-1.5,b4,'color','k');
   %text(tm+0.2,-1.8,'al =','color','k');
   text(tpt+0.01,-1.8,b5,'color','k');
   %text(tm+0.1,-2.2,'  dP =','color','k');
   text(tpt+0.01,-2.2,b6,'color','k');
   %text(tm+0.1,-2.5,' dPd =','color','k');
   text(tpt+0.01,-2.5,b7,'color','k');
   %text(tm+0.1,-2.8,'I2kvt =','color','k');
   text(tpt+0.01,-2.8,b15,'color','k');   
   hold on
   grid on
end  
for i=-3:0.001:5,
    plot(i,0,'k'); 
end;
for i=-3:0.001:5,
    plot(0,i,'k'); 
end;
%  y-Жё   m-Фи   c-Го   r-Кр   g-Зе   b-Си   w-Бе   k-Чё
