aci3_3_m.m

异步电机间接磁通量控制

% ACI3-3 system (Sensored Indirect Field Oriented Control of Induction Motor)
% 
% Required following function m-files:
% 1. "aci.m" - induction model
% 2. "cur_mod.m" - Current model
% 3. "pid_reg3.m" - PID controller
% 4. "park.m" - Park transformation
% 5. "inv_park.m" - Inverse park tranformation
% 6. "ramp_gen.m" - Ramp generator

clear all
close all

T = 5e-04;               % Sampling time for simulation (sec)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Incremental build %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phase1_inc_build = 0;   % Current_model test (CUR_MOD tests)
phase2_inc_build = 1;   % Closed-loop speed loop test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Incremental build %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%% Machine parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1-hp WEG Induction Motor
Rs = 1.723;                  % stator resistance (ohm) 
Rr = 2.011;                  % rotor resistance referred to stator (ohm)
Ls = (7.387+159.232)*1e-03;  % stator inductance (H)
Lr = (9.732+159.232)*1e-03;  % rotor inductance referred to stator (H)
Lm = 159.232*1e-03;          % magnetizing inductance (H)
P = 4;                       % number of poles
J = 0.001;                   % moment of inertia of the rotor mass (kg.m^2)
B = 0.0001;                  % damping coefficient (N.m.sec/rad) - always ignored
Tl = 1;                      % load torque (N.m)
np = P/2;                    % number of pole pairs
%%%%%%%%%%%%%%%%%%%%%%%%%% Machine parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Base quantities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fb = 60;                    % Base supplied frequency (Hz)
Wb = 2*pi*fb;               % Base electrically angular velocity (rad/sec)
Ib = 5;                     % Base phase current (amp)
Vb = 320/sqrt(3);           % Base phase voltage (volt)
%Lb = 220*sqrt(2/3)/(2*pi*60);     % Base flux linkage (volt.sec/rad)
Lb = Lm*Ib;                        % Base flux linkage (volt.sec/rad)
Tb = (3*Vb*Ib/2)*(np/(2*pi*60));   % Base torque (N.m)
SPb = 120*fb/P;                    % Base synchronous speed (rpm)  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Base quantities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%% Initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Induction motor
X = [0;0;0;0];       % X = [psi_rq; psi_rd; i_sq; i_sd] 
Wr = 0;              % Rotor electrically angular velocity (rad/sec)
Iq = 0;
Id = 0;
% Parameter vector of induction motor
p_im = [T; Rs; Rr; Ls; Lr; Lm; np; J; B; Tl/Tb; Wb; Ib; Vb; Lb; Tb];  

% Current model
h_cm = [0;0];              %  h = [im; theta_old]
p_cm = [T;Rr;Lr;fb];       %  p = [T; Rr; Lr; fb]
theta_psi_r = 0;

% PID - IQ control
h_iq = [0;0;0];            % h_pid  = [up_reg3; ui_reg3; ud_reg3]  
p_iq = [T; 2*Ib/Vb; 0.001; 0.0001; 0.1; 0.71; -0.71]; % [T; Kp; Ti; Td; Kc; Umax; Umin];  

% PID - ID control
h_id = [0;0;0];            % h_pid  = [up_reg3; ui_reg3; ud_reg3]  
p_id = [T; 2*Ib/Vb; 0.001; 0.0001; 0.1; 0.71; -0.71]; % [T; Kp; Ti; Td; Kc; Umax; Umin]; 

% PID - Speed control
h_sp = [0;0;0];            % h_pid  = [up_reg3; ui_reg3; ud_reg3]  
p_sp = [T; 0.02*SPb/Ib; 0.01; 0.0001; 0.9; 1; -1]; % [T; Kp; Ti; Td; Kc; Umax; Umin];  

%%%%%%%%%%%%%%%%%%%%%%%%%% Initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%% Simulation inputs  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Tt = 0.2;                 % Total time for simulation (sec) 
%Tt = 29*T;                 % Total time for simulation (sec) 
Tt = 1024*T;                 % Total time for simulation (sec) 

t = 0:T:Tt;               % Time for simulation (sec) 

speed_ref = 900/SPb;     % Reference speed (pu)
Id_ref = 2/Ib;            % Reference syn. rotating d-axis current (pu)   
Iq_ref = 0/Ib;            % Reference syn. rotating q-axis current (pu)  
%%%%%%%%%%%%%%%%%%%%%%%%%% Simulation inputs  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%% Simulation section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:length(t),

   if phase1_inc_build==1         
      % Phase 1 incremental build
     
      [rmp_out] = ramp_gen(2*pi/(2*pi),0,speed_ref*fb,t(i));
     [Uq_ref,Y,h_iq] = pid_reg3(Iq_ref,Iq,h_iq,p_iq);
     [Ud_ref,Y,h_id] = pid_reg3(Id_ref,Id,h_id,p_id);
     
     [uq_ref,ud_ref] = inv_park(Uq_ref,Ud_ref,rmp_out);
   
     [Te,Wr,X] = aci(Wr,X,[uq_ref; ud_ref],p_im);
   
     [Iq,Id] = park(X(3),X(4),rmp_out);

     [theta_psi_r,h_cm] = cur_mod(Id,Iq,Wr,h_cm,p_cm);     
         
     rmp(i) = rmp_out;
     
    elseif phase2_inc_build==1         
      % Phase 2 incremental build

     [Iq_ref,Y,h_sp] = pid_reg3(speed_ref,Wr,h_sp,p_sp);  % sensored
     [Uq_ref,Y,h_iq] = pid_reg3(Iq_ref,Iq,h_iq,p_iq);
     [Ud_ref,Y,h_id] = pid_reg3(Id_ref,Id,h_id,p_id);
  
     [uq_ref,ud_ref] = inv_park(Uq_ref,Ud_ref,theta_psi_r);
   
     [Te,Wr,X] = aci(Wr,X,[uq_ref; ud_ref],p_im);
  
     [Iq,Id] = park(X(3),X(4),theta_psi_r);

     [theta_psi_r,h_cm] = cur_mod(Id,Iq,Wr,h_cm,p_cm);     
      
 end
   
   i_qe(i) = Iq;
   i_de(i) = Id;
   v_qe(i) = Uq_ref;
   v_de(i) = Ud_ref;

   ibeta(i) = X(3);
   ialfa(i) = X(4);
   vbeta(i) = uq_ref;
   valfa(i) = ud_ref;
     
   psi_r_beta(i) = X(1);
   psi_r_alfa(i) = X(2);

   theta_r(i) = theta_psi_r;
   
   torque(i) = Te;
   wr(i) = Wr; 
   
   t(i)
     
 end   
%%%%%%%%%%%%%%%%%%%%%%%%%% Simulation section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
subplot(221)
plot(t,valfa,t,vbeta);
legend('v_{s\alpha}','v_{s\beta}',0);
title('Induction Motor Responses');
ylabel('stator voltages (pu)');
xlabel('time (sec)');
grid
subplot(222)
plot(t,v_de,t,v_qe);
legend('v_{de}','v_{qe}',0);
title('Induction Motor Responses');
xlabel('time (sec)');
grid
subplot(223)
plot(t,ialfa,t,ibeta);
legend('i_{s\alpha}','i_{s\beta}',0);
ylabel('stator currents (pu)');
xlabel('time (sec)');
grid
subplot(224)
plot(t,i_de,t,i_qe);
legend('i_{de}','i_{qe}',0);
xlabel('time (sec)');
grid

figure(2)
subplot(211)
plot(t,psi_r_alfa,t,psi_r_beta);
legend('\psi_{r\alpha}','\psi_{r\beta}',0);
title('Induction Motor Responses');
ylabel('Rotor Flux Linkages (pu)');
xlabel('time (sec)');
grid
subplot(212)
plot(t,torque);
ylabel('T_e (pu)');
xlabel('time (sec)');
grid

if phase1_inc_build==1 
figure(3)
subplot(211);
plot(t,rmp);
ylabel('\theta_{ramp} (pu)');
xlabel('time (sec)');
grid
subplot(212);
plot(t,theta_r);
ylabel('\theta_{\Psi_r} (pu)');
xlabel('time (sec)');
grid
end

if phase2_inc_build==1 
figure(3)
plot(t,theta_r);
title('Estimated rotor flux');
ylabel('\theta_{\Psi_r} (pu)');
xlabel('time (sec)');
grid
end

figure(4)
plot(t,speed_ref*ones(1,length(t)),':',t,wr);
legend('\omega^*_r','\omega_r',0);
title('Induction Motor Responses');
ylabel('rotor electrically angular velocity (pu)');
xlabel('time (sec)');
grid

figure(5)
subplot(211)
plot(t(1:29),wr(1:29));
title('Induction Motor Responses');
ylabel('Measured rotor speed (pu)');
xlabel('time (sec)');
grid
subplot(212)
plot(t(1:29),torque(1:29));
ylabel('Measured torque (pu)');
xlabel('time (sec)');
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(6)
subplot(211)
plot(t,valfa);
axis([0 Tt -0.5 0.5]);
title('Alpha-axis stator voltage (pu)');
grid
subplot(212)
plot(t,vbeta);
axis([0 Tt -0.5 0.5]);
title('Beta-axis stator voltage (pu)');
xlabel('time (sec)');
grid

figure(7)
subplot(211)
plot(t,ialfa);
axis([0 Tt -1.5 1.5]);
title('Alpha-axis stator current (pu)');
grid
subplot(212)
plot(t,ibeta);
axis([0 Tt -1.5 1.5]);
title('Beta-axis stator voltage (pu)');
xlabel('time (sec)');
grid

figure(8)
subplot(211)
plot(t,i_de);
axis([0 Tt 0 1]);
title('D-axis stator current (pu)');
grid
subplot(212)
plot(t,i_qe);
axis([0 Tt -0.5 1.5]);
title('Q-axis stator current (pu)');
xlabel('time (sec)');
grid

figure(9)
subplot(211)
plot(t,psi_r_alfa);
axis([0 Tt -0.5 0.5]);
title('Alpha-axis rotor flux (pu)');
grid
subplot(212)
plot(t,psi_r_beta);
axis([0 Tt -0.5 0.5]);
title('Beta-axis rotor flux (pu)');
xlabel('time (sec)');
grid

figure(10)
subplot(211)
plot(t,torque);
axis([0 Tt -0.5 0.5]);
title('Electromagnetic torque (pu)');
grid
subplot(212)
plot(t,wr);
axis([0 Tt -0.2 0.8]);
title('Rotor speed (pu)');
xlabel('time (sec)');
grid

figure(11)
subplot(211)
plot(t,ialfa);
axis([0 Tt -2 2]);
title('MATLAB Version');
ylabel('Alpha-axis stator current (pu)');
xlabel('time (sec)');
grid
subplot(212)
plot(t,wr);
axis([0 Tt -1 1]);
ylabel('Measured speed (pu)');
xlabel('time (sec)');
grid

figure(12)
subplot(211)
plot(t,torque);
axis([0 Tt -0.5 0.5]);
title('MATLAB Version');
ylabel('Electromagnetic torque (pu)');
xlabel('time (sec)');
grid
subplot(212)
plot(t,theta_r);
axis([0 Tt -1 1]);
ylabel('Computed rotor flux angle (pu)');
xlabel('time (sec)');
grid