script
        
        %Базовый для расчета переходных процессов в системе ПЧ-АДКЗ 
                       %с задатчиком интенсивности 
                      %при U1=var, f1=var, xm=const                      
  prog='zipchad';
%                    Перед расчетами ПОВТОРИТЬ главы:
%    "Регулирование скорости в системе ПЧ-АД асинхронного двигателя"
%    "Переходные процессы в системе ПЧ-АД с задатчиком интенсивности"
              
 %                         РАССЧИТАТЬ:
%              wona - Начальное значение синхронной скорости,рад/с
%              woko - Конечное значение синхронной скорости,рад/с
%              Tzi - постоянную времени задатчика интенсивности,с

%                               ВНИМАНИЕ!
%        При выполнении расчетов сохранить данную отладочную программу 
%                    под другим именем, а в этой программе
%                             НИЧЕГО НЕ ИЗМЕНЯТЬ!!!!!!!!!!!! 

%                            РАБОТА С ПРОГРАММОЙ
%       1. Внимательно прочитайте свою программу, уясните алгоритм расчета.
%          Обозначения физических величин, примененные в программе, часто не 
%          совпадают с общепринятыми. Следует разобраться с их обозначениями.
%       2. Дифференциальные и алгебраические уравнения программы описаны
%          в учебном пособии к курсовому проектированию по "Теории электропривода".
%       3. Ввести перечисленные параметры двигателя, преобразователя
%       частоты ПЧ, рабочего органа РО;
%       4. Рассчитать и ввести
%         - постоянную времени задатчика интенсивности; 
%         - скачок синхронной скорости при ПИ-задатчике интенсивности - при необходимости; 
%         - напряжение форсировки,о.е. - при необходимости; 
%         - компенсация скольжения(0...1) - коэффициент обратной связи 
%            по активной составляющей тока статора  - при необходимости;    
%        - начальные wona и конечные woko синхронные скорости; 
%          - скорости wna вала  двигателя;
%          - шаг интегрирования h (можно не менять) ;  
%          - предполагаемые времена пуска tk с выходом на установившийся режим 
%          и торможения tt с достижением остановки двигателя;
%        5. При  включении программы ПЕРВЫМ  следует выбрать переходный процесс
%           пуска nu=0 и нажать знак "Save and Run".
%            По окончании  пуска программа  выводит
%            figure 1  -  нагрузочные диаграммы c интегральными  показателями пуска и
%            figure 2 - динамические механические характеристики.   
%        6. Не убирая с экрана figure1 и figure2, перейти к расчету торможения,     
%           установив nu=1 и вновь нажать знак "Save and Run".
%            По окончании торможения 
%           figure1 и figure2 будут дополнены данными переходного процесса торможения. 
%        7. Сохраните figure1 и figure2 для печати и дальнейшей обработки
%           результатов расчетов (проверок по нагреву, перегрузочной
%           способности, энергетических показателей электропривода и др.)
%        8. Повторите пп. 3...7 для второго режима работы.

%         При  ПЕРВОМ  включении программы следует выбрать 
%                                    переходный процесс пуска      nu=0                                                 
                                                    
%Каталожные данные двигателя
tip='4MTKF(H)200L6'; %Тип двигателя
pp=3;            %Число пар полюсов
u1n=220;         %Номинальное фазное напряжение,В
f1n=50;          %Номинальная частота,Гц
Pn=22000;        %Номинальная мощность, Вт
nn=935;          %Номинальная скорость,об/мин
I1n=51;          %Номинальный ток статора,А 
kpdn=0.83;       %Номинальный кпд
cosfi1n=0.79;    %Номинальный cos fi
Mk=760;          %Критический момент, Нм
Mp=706;          %Пусковой момент,Нм
I1p=275;         %Пусковой ток статора, A
Jd=0.57;         %Момент инерции двигателя,кгм^2
      
%Обмоточные данные двигателя
r1=0.29;         %Невыключаемое активное сопротивление статора,Ом
Rk=0.5;         %Активное сопротивление к.з.,Ом
Xk=0.63;         %Индуктивное сопротивление к.з.,Ом
r2p=Rk-r1;       %Невыключаемое активное сопротивление ротора,Ом
x1=Xk/2;         %Невыключаемое индуктивное сопротивление статора,Ом
x2p=Xk/2;        %Невыключаемое индуктивное сопротивление ротора,Ом
 %Imn=31.3;        %Ток холостого хода, А 

%Данные рабочего органа
Mcx=50;           %Статический момент на валу на холостом ходу,Нм
Mcn=220;        %Статический момент на валу под нагрузкой,Нм
Mca=0;             %момент статический активный Мса=0,реактивный Мса=1
Jpr=0.525;      %Приведенный к валу двигателя момент инерции РО, кгм^2 
J=1.2*Jd+Jpr;%Момент инерции электропривода,кгм^2

%Каталожные данные преобразователя частоты AT04-37
Upn=380;
Ipn=75;
dPpn=2351;  % 3*220*75*0.05*0.95=2351.25
f1min=5;
f1max=100;
  vu=0;       %Тип выпрямителя преобразователя частоты:    
                    %vu=0-управляемый; vu=1-неуправляемый
                    
%Расчетные данные двигателя
pi=3.14;
won=2*pi*f1n/pp; %Номинальная синхронная скорость,рад/с
wn=nn/9.55;      %Номинальная скорость,рад/с
Mn=Pn/wn;        %Номинальный момент на валу,Нм
sn=(won-wn)/won; %Номинальное скольжение
Mko=Mk/Mn;       %Критический момент, о.е.
Mpo=Mp/Mn;       %Пусковой момент, о.е.
I1po=I1p/I1n;    %Пусковой ток статора, о.е.
r2v=Mp*won/(3*I1p^2); % Активное сопротивление ротора c учетом вытеснения тока,Ом 
a1=r1/r2p;
sk=sn*(Mko+sqrt((Mko^2-1)+2*a1*sn*(Mko-1))/(1-2*sn*a1*(Mko-1)));
Memn=2*Mk*(1+a1*sk)/(sn/sk+sk/sn+2*a1*sk); %Номинальный электромагнитный момент,Нм
dMxn=Memn-Mn;    %Потери момента в номинальном режиме, Нм
Te=(x1+x2p)/(314*(r1+r2p)); %Электромагнитная постоянная времени,с
sinfi1n=sqrt(1-cosfi1n^2);
Imn=I1n*(sinfi1n-cosfi1n*sn/sk);%Ток намагничивания  в номинальном режиме, А
xmn=u1n/Imn-x1;  %Индуктивное сопротивление контура намагничивания,Ом
                         
  % Расчетные данные для системы ПЧ-АД   
Tzi=1;        %Постоянная времени задатчика интенсивности,с.    0.5
dwo=20;    %Скачок синхронной скорости при ПИ-задатчике,рад/с
u1for=0;    %Напряжение форсировки,о.е.
kot=0;       %Компенсация скольжения(0...1)-Коэффициент обратной связи
                  %по активной составляющей тока статора                   
 womin=won*f1min/f1n;
 Td=J*won/Memn;  %Механическая постоянная времени электропривода,с
 
%НАЧАЛЬНЫЕ УСЛОВИЯ 

 h=0.005                           %Шаг интегрирования,с
 tk=1.5                               %Время пуска,с  
 tt=0.8;                               %Время торможения,с
 
                                            nu=0   %      nu=0-пуск; nu=1-торможение             
               
if nu==0        %     пуск;
         
   wonap=0;           %Начальное значение синхронной скорости, рад/с
   wokop=90;         %Конечное  значение синхронной скорости,рад/с  
   wna=0;                %Начальное значение скорости двигателя,рад/с
   w=wna;  
  
   Pv=0;Qd=0;Pd=0;Pc=0;Qc=0;I1kvt=0;al=0;km1=0;dPp=0;dPR=0;
   a=0; M=0; I1=0; u1=0; xm=xmn; wo=0; km1=0;km2=0; I1a=0;
   ks=0;kpdd=0;kpd=0;tz=0; wozi=0;wonat=wo; 
   
   to=h;                  %Расчет цикла пуска   
   t2=tk;  
end
    
if nu==1 % торможение
            
    wonat=wo;              %Начальное значение синхронной скорости,рад/с
   wokot=0;                 %Конечное  значение синхронной скорости,рад/с
 if wokot<womin  wokot=womin; end 
   wna=w;                  %Начальное значение скорости двигателя,рад/с
   w=wna;  
              
   Pv=0;Qd=0;Pd=0;Pc=0;Qc=0;I1kvt=0;al=0;km1=0;dPp=0;dPR=0;ks=0;
   a=a; M=M; I1=I1; u1=u1; w=wk; km1=0; km2=0;
   tz=0; wozi=0;
   
  to=tk,                       %Расчет цикла торможения
  t2=tk+tt;                                                     
end 
   
 %      РАСЧЕТ ПЕРЕХОДНОГО ПРОЦЕССА
   
for t = to:h:t2  
    
    t1=t 
          
      if and (nu==0, wokop>0), 
       s2=sign(wokop-wonap);;tz=t1;
       wozi=abs(won)*s2*tz/Tzi;     
    if wonap< womin   wonap=womin; end
    if wonap<s2*dwo     wonap=0;
    else wonap = wonap-s2*dwo;  end   
       wo1=wonap+s2*dwo+wozi;
   if wo1>(wokop+s2*dwo), wo1=wokop;end
      end
                       
      if and (nu==0, wokop<0), 
       s2=sign(wokop-wonap);;tz=t1; 
       wozi=abs(won)*s2*tz/Tzi;     
    if wonap>womin   wonap=-womin; end
    if wonap>s2*dwo     wonap=0;
    else wonap = wonap-s2*dwo;  end   
       wo1=wonap+s2*dwo+wozi;
   if wo1<(wokop+s2*dwo), wo1=wokop;end
   end 
         
     if and (nu==1, wonat<0), 
       s2=sign(wokot-wonat);    tz=t1-tk; 
       wozi=abs(won)*s2*tz/Tzi;    
       wo1=wonat+s2*dwo+wozi;
    if wokot >-womin   wokot=-womin; end
    if wo1>(wokot+s2*dwo), wo1=wokot;end
     end
                  
      if and (nu==1, wonat>0), 
       s2=sign(wokot-wonat);    tz=t1-tk; 
       wozi=abs(won)*s2*tz/Tzi;    
       wo1=wonat+s2*dwo+wozi;
    if wokot <womin   wokot=womin; end
    if wo1<(wokot+s2*dwo), wo1=wokot;end
   end
       
     %Приращение частоты от обратной связи 
%                       по активной составляющей тока статора
    dfi=f1n*I1a/(I1n*cosfi1n);
    ks=kot*sn*dfi;
    if w>wo, ks=0; end 
    
     wo=wo1+ks;    
    a=abs(wo/won);
    u1=(u1for+a-u1for*a)*u1n;
    if u1>u1n u1=u1n;end

%Расчет сопротивлений, токов, напряжений
%             Т-образной схемы замещения АД

    s=(wo-w)/wo;
    z1=r1+j*x1*a;
    if (r2v-r2p)>0 
        r2s=(r2p+(r2v-r2p))/s; %Учет вытеснения тока ротора
    else 
        r2s=r2p/s; 
    end
    z2p=r2s+j*x2p*a;
    zm=j*xm*a;
    zmr=zm*z2p/(zm+z2p);
    zc=z1+zmr;
    fi1=angle(zc);
    i1=u1/zc;
    I1u=abs(i1);
    I1=I1u*(1-exp(-t/Te));
    I1a=I1*cos(fi1);
    I1r=I1*sin(fi1);
    e=i1*zmr;
    i2p=e/z2p;
    I2pu=abs(i2p);
    I2p=I2pu*(1-exp(-t/Te));

%Уравнение движения методом Эйлера
   dMx=abs(dMxn)*sign(w);
   M=3*I2p^2*r2s/wo;
   if Mca==0,
      Mc=Mcn+dMx;             %Активный Мс,Нм
   else 
      Mc=abs(Mc1)*sign(w)+dMx; %Реактивный Мс,Нм
   end  
   w=w+h*(M-Mc)/J ;                 %Скорость двигателя,рад/с
   al=al+w*h;                        %Угол поворота вала двигателя,рад

% Энергетические показатели  
   Pv1=(M-dMx)*w;                          %Мощность на валу,Вт
   Pv=Pv+Pv1*h;                             %Активнаяэнергия на валу,Вт
   Pd1=3*u1*I1a;                          %Активная мощность двигателя,Вт
   Pd=Pd+Pd1*h;                           %Активная энергия двигателя,Вт*с
   Qd=Qd+3*u1*I1r*h;                  %Реактивная энергия двигателя,Вар*с
   I1kvt=I1kvt+I1^2*h;                    %Квадрат тока статора*время,А^2*с
   dPp1=dPpn*0.5*(1+(I1a/Ipn)^2);    %Потери мощности в ПЧ,Вт
   dPp=dPp+dPp1*h;                      %Потери энергии в ПЧ,Вт*с
   Pc1=Pd1+dPp1;                        %Активная мощность сети,Вт
   Pc=Pd+dPp;                             %Активная энергия сети,Вт*с
   
% Энергетические показатели: КПД kpd и коэффициент мощности km

  % Энергетические показатели: КПД kpd и коэффициент мощности km
             % nu=0-пуск,
             % nu=1-торможение
             % vu==0 - управляемый выпрямитель                                                                                    
             % vu=1-неуправляемый выпрямитель,  
            
       if and (Pc<0, vu==1)        %nu==1
          Pc=0; kpd=0; dPR=-(dPp+Pd);
       else
       % Цикловый КПД системы ПЧ-АД 
      if and (Pc<0,Pv<0)   
          kpd=Pc/Pv; 
      else
       if and (Pc>0,Pv>0)  
           kpd=Pv/Pc;  
       else
            if and (kpd<0, kpd>1)     kpd=0;
             end
       end
      end
       end
      
      %Мгновенный КПД системы ПЧ-АД
        if and (Pc1>0,Pv1>0)  kpdd=Pv1/Pc1; 
        else    
      if and (Pc1<0,Pv1<0) kpdd=Pc1/Pv1;
      else
          if and (kpdd<0, kpdd>1), kpdd=0;      
      end 
      end
       end     
       
      %Мгновенный коэффициент мощности системы ПЧ-АД 
     if    vu==0    % , vu=0 - управляемый выпрямитель    
          km=u1/(0.608*Upn); %km=cos(alfa)=Ed/Edo=u1*1.73*cosfi pchn / Upn=u1/(0.608*Upn)
     end
    if  vu==1  % vu=1-неуправляемый выпрямитель
        km=0.95; %cosfi=1, коэффициент искажений kи=0,95
        end   
         km1=km1+km*h; %Цикловый коэффициент мощности системы ПЧ-АД  
         Qc=Pc*sqrt(1-km^2)/km;      %Реактивная энергия ,Вар*с
               
     if nu==0
      I1kv=sqrt(I1kvt/tk); %Среднеквадратичный ток статора
       km2=km1/tk;  
   else 
      I1kv=sqrt(I1kvt/(tt));
      km2=km1/tt;   
   end

%Печать результатов расчета

%Вывод текста на на figure 1
        tm=0.55;  %начало вывода интегральных  показателей пуска 
        %ty=Tzi;  %начало вывода показателей уст. режима на графике
        %ttt=tk+tt; %начало вывода интегральных  показателей торможения 
   b0=num2str(t);
   b1=num2str(Pv);
   b2=num2str(Pc);
   b3=num2str(Qc);
   b4=num2str(I1kv);
   b5=num2str(al);
   b6=num2str(dPp);
   b7=num2str(Pd);
   b8=num2str(a);
   b9=num2str(w);
   b10=num2str(M);
   b11=num2str(I1);
   b12=num2str(u1);
   b13=num2str(Td);
   b14=num2str(dwo);
   b15=num2str(u1for);
   b16=num2str(kot);
   b18=num2str(Pn*10^(-3));
   b19=num2str(won);
   b20=num2str(Memn);
   b21=num2str(I1n);
   b22=num2str(t-tk);
   b23=num2str(kpdd);
   b24=num2str(km);

   % Если какую-то из характеристик не нужно выводить на график,
   % поставьте знак % перед plot:   %plot()
   
   figure(1)         % - нагрузочные диаграммы w,M,I1,wo=f(t)
 
   xlabel('t')
   ylabel('w,M,I1,wo,kpdd,km'),   
   plot(t,wo1/abs(won),'r')
  plot(t,I1/I1n,'m')
  plot(t,w/abs(won),'b')
  plot(t,M/Memn,'k')
  plot(t,kpdd,'c')
   plot(t,km,'g')
   axis([0 4.5 -2.5 3.2])      %рамка
   hold on 
   grid on
   %end
    
 
      figure(2)            % Динамические механические характеристики

   xlabel('M,I1')
   ylabel('w')
   plot(M/Memn,w/abs(won),'k');      
   plot(I1/I1n,w/abs(won),'m');     
   axis([-4 6 -1.5 1.5])       %рамка 
   hold on 
   grid on
end
%end
for i=-3:0.01:5.5, 
    plot(i,0,'k'); %ось х  
end

for i=-1.4:0.01:1.4, 
    plot(0,i,'k'); %ось y 
end 
wk=w;   %Конечная скорость при пуске

if and(figure(1),nu==0),%nu=0-пуск 
         
   d=date;
   text(2.55,-2.2,d,'color','k');
   text(1.55,-2.2,prog,'color','k');
   text(0.55,-2.2,tip,'color','k');
   text(tm-0.25,-0.15,'tk =','color','k');
   text(tm,-0.15,b0,'color','k');    
   text(tm-0.3,-0.4,' Pv =','color','k');
   text(tm,-0.4,b1,'color','k');
   text(tm-0.3, -0.65,'Pd =','color','k');
   text(tm,-0.65,b7,'color','k');
   text(tm-0.35,-0.9,'dPp =','color','k');
   text(tm,-0.9,b6,'color','k');
   text(tm-0.3,-1.15,'Pc =','color','k');
   text(tm,-1.15,b2,'color','k');
   text(tm-0.3,-1.4,'Qc =','color','k');
   text(tm,-1.4,b3,'color','k');
   text(tm-0.35,-1.65,'I1kv =','color','k');
   text(tm,-1.65,b4,'color','k');
   text(tm-0.25,-1.9,'al =','color','k');
   text(tm,-1.9,b5,'color','k');

   text((Tzi+0.2),2.85,'wo/won =','color','r');
   text((Tzi+0.65),2.85,b8,'color','r');
   text((Tzi+0.2),2.6,'w =','color','b');
   text((Tzi+0.55),2.6,b9,'color','b');
   text((Tzi+0.2),2.35,'M =','color','k');
   text((Tzi+0.55),2.35,b10,'color','k');
   text((Tzi+0.2),2.1,'I1 =','color','m');
   text((Tzi+0.55),2.1,b11,'color','m');
   text((Tzi+0.2),1.85,'kpdd =','color','c');
   text((Tzi+0.55),1.85,b23,'color','c');
   text((Tzi+0.2),1.6,'km =','color','g');
   text((Tzi+0.55),1.6,b24,'color','g');

   text((tk+tt+0.45),2.85,' Td =','color','k');
   text((tk+tt+0.85),2.85,b13,'color','k');
   text((tk+tt+0.45), 2.6,'dwo =','color','k');
   text((tk+tt+0.85),2.6,b14,'color','k');
   text((tk+tt+0.35),2.35,' u1for =','color','k');
   text((tk+tt+0.85),2.35,b15,'color','k');
   text((tk+tt+0.45),2.1,'kot =','color','k');
   text((tk+tt+0.85),2.1,b16,'color','k');
end 

if and(figure(1),nu==1), %nu=1-торможение
 
   text(tk+tt+0.1,-0.15,'tt =','color','k');
   text(tk+tt+0.45,-0.15,b22,'color','k');    
   text((tk+tt+0.45),-0.4,b1,'color','k');
   text((tk+tt+0.45),-0.65,b7,'color','k');
   text((tk+tt+0.45),-0.9,b6,'color','k');
   text((tk+tt+0.45),-1.15,b2,'color','k');
   text((tk+tt+0.45),-1.4,b3,'color','k');
   text((tk+tt+0.45),-1.65,b4,'color','k');
   text((tk+tt+0.45),-1.9,b5,'color','k');
end

if and(figure(2),nu==0),

   d=date;
   text(3.8,-0.3,d,'color','k');
   text(3.8,-0.2,prog,'color','k');
   text(-3.8,1.3,tip,'color','k');
   text(-1.4,1.3,'Pn =','color','k');
   text(-1,1.3,b18,'color','k'); 
   text(0.2,1.3,'won =','color','k');
   text(1,1.3,b19,'color','k');
   text(2.7,1.2,' I1n =','color','m');
   text(3.6,1.2,b21,'color','m');
   text(2.4,1.3,'Memn =','color','k');
   text(3.6,1.3,b20,'color','k');

end
for i=0:0.01:5.5, 
    plot(i,0,'k'); %ось х  
end 
 %end
%  y-Жё   m-Фи   c-Го   r-Кр   g-Зе   b-Си   w-Бе   k-Чё




