imrunsi1.m

关于电机的.m程序

%                <<<   INDUCTION MOTOR TRANSIENTS  >>>
%          Motor is initialized at a specified operating condition
%              Choice of reference frame: Stator, rotor, or synchronous
close all,clear,clc
disp(date)
info= ...
        ['                   INDUCTION MOTOR TRANSIENTS                                   '                                
         '                     [Sinusoidal excitation]                                    '  
         '                           (SI units)                                           '
         ' The motor is running in steady state at a light load (established by a speed   '  
         ' close to synchronous) when rated load torque is applied; when the new steady-  '  
         ' state conditions are reached, a light load is suddenly impressed followed by   '
         ' a 5/6 step reduction in both applied voltage and frequency .                   '
         ' The user has the option of selecting one of three common (dq)reference frames: '
         ' STATOR, ROTOR, or SYNCHRONOUS.                                                 ' 
         ' The torque, speed,and current responses are plotted.                           '
         ' First, the time evolution of speed (electric rad/s)is plotted as the simulation'
         ' proceeds and it can be stopped at any time.                                    '
         ' *Refer to script files <imrunsi1.m> and <imrunsif.m>.                          '];
      disp(info)
      disp('                                                         ')
% ---------------------------  DATA  ------------------------------------------
global Rs Rr L J Vr fs fk wo km Tl p 

% --------------------     MACHINE PARAMETERS ------------------------------ 
Vll=220    ; %Line-to-line rms voltage rating [V]
Rs=0.22    ; %Stator resistance (Ohm)
Rr=0.45    ; %Rotor resistance  (Ohm)
Xsl=2.2    ; %Stator leakage reactance (Ohm)
Xrl=2.2    ; %Rotor leakage reatance  (Ohm)
Xm=98.55   ; %Magnetizing reactance  (Ohm) 
fs=1.0     ; %Stator frequency (pu)
freq0=60   ; %Base frequency (Hz)
wo=2*pi*freq0 ; %Base frequency (rad/s)
p=4        ; %number of poles
J=0.05     ; %Moment of inertia (kg.m^2)
n=0.998    ; %per-unit operating speed (adjust to modify initial loading)
wm=n*wo    ; %initial rotor speed in electrical rad/s
Vr=Vll*sqrt(2/3); %peak per-phase voltage [V]
Vs=Vr      ;
% --------------------  REFERENCE FRAME SELECTION  ----------------------------
      km=menu('Select a common reference frame','STATOR','ROTOR','SYNCHRONOUS','QUIT');
                if km==1,  fk=0;                  
                   elseif km==2,  fk=n;           
                   elseif km==3, fk=fs;
                   elseif km<=0 | km>3, break
                end
% ---------------- STEADY-STATE PHASOR CALCULATIONS OF OPERATING POINT ------
     disp('          Parameters of the induction machine in SI units    ')
     disp('    Rs        Xsl       Rr        Xrl       Xm        Vs          J    ')
param=[Rs Xsl Rr Xrl Xm Vs J];
disp(param)
Lsl=Xsl/wo ;Lrl=Xrl/wo; Lm=Xm/wo;
Ls=Lsl+Lm ; Lr=Lrl+Lm;
L=[Lsl+Lm   0      Lm     0   
       0   Lsl+Lm    0     Lm
       Lm    0     Lrl+Lm  0
       0    Lm      0     Lrl+Lm];
   j=sqrt(-1);
    Zr=Rr/(fs-n)+j*Xrl; Zm=j*Xm*Zr/(j*Xm+Zr);  Zin=Rs/fs+j*Xsl+Zm;
    Is=Vs/(fs*Zin);  Ir=-Is*j*Xm/(j*Xm+Zr);
    Req=real(Zm);    Tl=0.75*p/wo*Is*conj(Is)*Req;Tl=real(Tl);
    phis=Ls*Is+Lm*Ir;  phir=Lm*Is+Lr*Ir;    phids=real(phis);  phiqs=imag(phis);
    phidr=real(phir);  phiqr=imag(phir);
    spd0=wm*60/(p*pi);
    fprintf(1,'  The operating point is established at: %6.0f rpm and %2.4f N.m\n',spd0,Tl)
% ---------------------------  ODE  -------------------------------------------
 fo=[phids phiqs phidr phiqr wm 0 0]';        %initial conditions 
 tf=1.4  ;                                    % final time
 tspan=[0 tf] ; 
 options=odeset('RelTol',1e-2,'maxstep',1e-3,'outputfcn','odeplot','outputsel',5) ;
 [t,f]=ode45('imrunsif',tspan,fo,options);
 % --------------------  OUTPUT VECTORS  --------------------------------------
 flux=f(:,1:4)'  ; cur=L\flux ;
 spd=f(:,5)*60/(p*pi);  %rpm
 phidr=f(:,3); phiqr=f(:,4);
 ids=cur(1,:)'  ;  iqs=cur(2,:)' ; is=sqrt(ids.^2+iqs.^2);
 Te=0.75*p*(iqs.*f(:,1)-ids.*f(:,2));
 ias=ids.*cos(f(:,6))-iqs.*sin(f(:,6));
 ibs=ids.*sin(f(:,6))+iqs.*cos(f(:,6));
  % --------------------  GRAPHICAL DISPLAYS  ----------------------------------
  h0=figure('Position',[100 60 900  560],'Name','Induction motor transients');
   subplot(2,2,1),plot(t,spd);title('Speed vs time');xlabel('Time [s]'),ylabel('Speed [rpm]'),grid
   subplot(2,2,2),plot(t,Te,'b');title('Torque versus time');xlabel('Time [s]'),ylabel('Torque [N.m]'),grid
   subplot(2,2,3),plot(spd,Te,'k');title('Torque versus speed');xlabel('Speed [rpm]'),ylabel('Torque [N.m]'),grid
   subplot(2,2,4),plot(t,phidr,'g',t,phiqr,'b');title('fdr&fqr vs time');xlabel('Time [s]'),ylabel('Flux [Wb]'),grid
      pause(4)
  h1=figure('Position',[100 60 900  560],'Name','Induction motor transients');
   subplot(2,2,1),plot(t,ids,'r',t,iqs,'b');title('ids&iqs vs time');xlabel('Time [s]'),ylabel('Stator current [A]]'),grid
   subplot(2,2,2),plot(t,ias,'r',t,is,'k');title('ias versus time');xlabel('Time [s]'),ylabel('Stator current [A]'),grid
   subplot(2,2,3),plot(ids,iqs,'b');title('Current space vector trajectory in dq frame');grid
   axis equal
   subplot(2,2,4),plot(ias,ibs,'r'),title('Current space vector trajectory in stator frame'),grid
  axis equal
   uicontrol('Parent',h1, ...
          'Units','normalized', ...
	      'BackgroundColor',[0.5 0.5 0.5], ...
	      'Callback','clc,clear,close all,break', ...
	      'FontSize',10, ...
	      'FontWeight','bold', ...
          'ForegroundColor',[0 0 0.5], ...
 	      'Position',[0.88 0.01 0.08 0.071], ...
	      'String','CLOSE', ...
	      'TooltipString','Close');  
% ------------------------------------------------------------------------------