% EE452 Lab 12
% AC Induction Machine Parameter Measurements
% Data Collection and Calculation script
% Revision 3/30/12
% Nicholas David
 
% Each group will be given a machine to measure the parameters of.
% Machines to be measured:
%   Marathon 5K33GN2A (SCIM used in Labs)
%   Lab-Volt Wound Rotor Machine (SCIM used in Labs)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is a script to measure parameters of an induction machine.  
%This procedure has been extracted from information contained in
%IEEE Standard 112-2004.  In particular, this script is used to implement
%test via "Method 4" described in section 5.9.5 of the standard.

%Circuit model contains:
% Stator Resistance
% Stator Leakage Reactance
% Rotor Resistance
% Rotor Leakage Inductance
% Magnetizing Inductance
% Core Resistance

%Calculation of model parameters is performed via measurements of
%rotor speed, line-to-line voltage phasor, and line current phasor, for
%No-Load, Locked-Rotor, and Load-Point operating conditions.  All
%measurements are performed with a three phase source derived from the utility
%grid (60Hz).

%These tests can be applied to a squirrel-cage machine, or a wound rotor
%machine.  For wound rotor machine, short circuit rotor windings for all
%tests.  For additional details, refer to IEEE Standard 112-2004.

%User Inputs are denoted by %**

clear all, clc;

%Machine quantities and test conditions
f = 60; %**Stator frequency applied for No Load Test (Hz)
P = 4; %**Number of machine poles
m = 3; %**Number of phases
ns = (120*f/P); %Synchronous speed for machine under test (rpm)
TF = 72; %**Room Temp (deg F)
NemaRatio = 0.67; %**Ratio of X1 to X2 (X1/X2 = 0.67 for NEMA Design B) 
XmRatio = 1; %**Define a ratio (X1/Xm) for use in calculations 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 1 - Stator Resistance, R1
%Directly measure stator winding resistance with Digital Multi-Meter
%Measure value of all three windings and take average
%This resistance includes all cabling and connections
%When Stator is connected in WYE, measurement between terminals is twice
%the resistance of each winding.


TC = (TF - 32)*5/9;
ta = TC; %**Enter temperature of coil at time of measurement; ambient temp (deg C)
tb = 40; %**Temperate at which resistance is reported; value on nameplate (deg C) 
k1 = 234.5; %factor for 100% IACS conductivity copper (inferred temperature for zero resistance)

Rab = ; %**Enter measured value (ohm)
Rbc = ; %**Enter measured value (ohm)
Rca = ; %**Enter measured value (ohm)
R1 = (Rab/2 + Rbc/2 + Rca/2)/3; %Average stator resistance (ohm)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 2: Do No-Load and Locked Rotor tests for determination of Xm, X1, X2

%NO NOAD TEST
%Run machine at rated voltage and rated frequency

%Measure stator terminal voltage, line current, phase, rotor speed
VabcLL0 = 2*[104 103 103]; %**[Vab Vbc Vca] %Measured line-to-line voltage magnitude (Vrms)
IabcL0 = [0.928 0.905 0.914]; %**[Ias Ibs Ics] %Measured line current magnitude (Arms)
tabc0 = (10^-3)*[5.1 5.1 5.1]; %**[ta tb tc] %time measured between line-to-line voltage and current peaks (s)
nt0 = 1795; %**Rotor speed during No Load Test (rpm)

%Condition Data
VabcLN0 = VabcLL0/sqrt(3); %Line-to-Neutral voltage magnitdue
VabcM0 = sum(VabcLN0)/3; %Average Line-to-Neutral voltage magnitude
I10 = sum(IabcL0)/3; %Average Line current magnitude
alphaabc0 = tabc0/(1/60)*2*pi; %Angle measured between line-to-line voltage and line current (rad)
thetaabc0 = alphaabc0 - pi/6;  %Angle between line-to-neutral voltage and line current (rad)
theta0 = sum(thetaabc0)/3; %Average angle between line-to-neutral voltage and line current (rad)
PF0 = cos(theta0); %Average power factor

%Define voltage phasor as reference
Va0 = VabcLN0(1).*exp(i*0);
Vb0 = VabcLN0(2).*exp(i*(-2*pi/3)); 
Vc0 = VabcLN0(3).*exp(i*(+2*pi/3)); 
Vabc0 = [Va0 Vb0 Vc0]; %Vabc line-to-neutral voltage phasors

%Formulate Line Current phasors
Ia0 = IabcL0(1).*exp(i*(-thetaabc0(1)));
Ib0 = IabcL0(2).*exp(i*(-thetaabc0(2) - 2*pi/3));
Ic0 = IabcL0(3).*exp(i*(-thetaabc0(3) + 2*pi/3));
Iabc0 = [Ia0 Ib0 Ic0]; %Iabc line current phasors

%Determine slip
SS0 = ns - nt0; %Slip Speed (rpm)
s0 = SS0/ns; %slip speed (pu)

%Calculate input Real and Reactive Power
S0NL = Vabc0 .* conj(Iabc0);  % per phase
P0 = real(sum(S0NL)); % three-phase 
Q0 = imag(sum(S0NL)); % three-phase


%LOCKED ROTOR TEST
%Run machine at rated frequency, but prevent rotor from turning
%Increase voltage until line current is at rated value

%Measure stator terminal voltage, line current, phase, rotor speed
VabcLLLR = [36.1 35.6 36.3]; %**[Vab Vbc Vca] %Measured line-to-line voltage magnitude (Vrms); range is low enough to use 1:20 attenuation on probe
IabcLLR = [1.60 1.69 1.65]; %**[Ias Ibs Ics] %Measured line current magnitude (Arms)
tabcLR = (10^-3)*[3.4 3.4 3.4]; %**[ta tb tc] %time measured between line-to-line voltage and current peaks
ntLR = 0; %**Rotor speed during Locked Rotor Test (rpm)

%Condition Data
VabcLNLR = VabcLLLR/sqrt(3); %Line-to-Neutral voltage magnitdue
V1L = sum(VabcLNLR)/3; %Average Line-to-Neutral voltage magnitude
I1L = sum(IabcLLR)/3; %Average line current magnitude
alphaabcLR = tabcLR/(1/60)*2*pi; %Angle measured between line-to-line voltage and line current (rad)
thetaabcLR = alphaabcLR - pi/6; %Angle between line-to-neutral voltage and line current (rad)
thetaLR = sum(thetaabcLR)/3; %Average angle between line-to-neutral voltage and line current (rad)
PFL = cos(thetaLR); %Average power factor

%Define voltage phasor as reference
VaL = VabcLNLR(1).*exp(i*0); 
VbL = VabcLNLR(2).*exp(i*(-2*pi/3));
VcL = VabcLNLR(3).*exp(i*(+2*pi/3)); 
VabcL = [VaL VbL VcL];

%Formulate Line Current phasors
IaL = IabcLLR(1).*exp(i*(-thetaabcLR(1)));
IbL = IabcLLR(2).*exp(i*(-thetaabcLR(2) - 2*pi/3));
IcL = IabcLLR(3).*exp(i*(-thetaabcLR(3) + 2*pi/3));
IabcL = [IaL IbL IcL];

%Calculate Slip
SSLR = ns - ntLR; %Slip Speed (rpm)
sLR = SSLR/ns; %slip speed (pu)

%Calculate input Real and Reactive Power
SL = VabcL .* conj(IabcL);  %per phase
PL = real(sum(SL)); %three-phase 
QL = imag(sum(SL)); %three-phase


%Do Iterative Calculations to determine Xm, X1, and X2
%(Calculations will be iterated until convergense is reached)
%Define initial values of X1, Xm and X2
X1 = 2*pi*f*(0.001); %**(2*pi*f*Lohm)
Xm = X1/XmRatio; %(ohm)
X2 = X1/NemaRatio; %(ohm)
X1prev = X1; %Define initial values used for determination of convergence
Xmprev = Xm; %Define initial values used for determination of convergence

%Start Solving for X1 and Xm
Xm = m*VabcM0^2/(Q0-(m*I10^2*X1)) * 1/((1+X1/Xm)^2); %(eqt 30 in std. 112)
X1L = QL/(m*I1L^2*(1+X1/X2+X1/Xm)) * (X1/X2+X1/Xm); %(eqt 31 in std. 112)
X1 = 60/f*X1L;
%Check if new value is more than 0.1% different from previous value
while or((or((X1>1.001*X1prev), (X1<0.999*X1prev))), (or((Xm>1.001*Xmprev), (Xm<0.999*Xmprev))));
     X1prev = X1; %Save old value
     Xmprev = Xm; %Save old value
     Xm = m*(VabcM0^2)/(Q0-(m*(I10^2)*X1)) * 1/((1+X1/Xm)^2); %Calculate new value
     X1L = QL/(m*I1L^2*(1+X1/X2+X1/Xm)) * (X1/X2+X1/Xm); %Calculate new value
     X1 = 60/f*X1L; %Calculate new value  
end    
% end the iterations and keep final value

X1 = X1; %(ohm)
Xm = Xm; %(ohm)
X2 = X1/NemaRatio; %(ohm)

%Convert Reactance to Inductance
L1 = X1/(2*pi*f); %(H)
Lm = Xm/(2*pi*f); %(H)
L2 = X2/(2*pi*f); %(H)
 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 3: %Do No-Load test at range of voltages to determine friction and windage
%Range of voltage is from 125%rated down to point at which current magnitude increases

%Record voltage and current at voltage 
VFWLL = 2*[105 105 105; 100 99.7 100; 94.0 93.7 94.2; 88.2 87.8 88.2; 82.7 82.2 82.7; 72.6 72.1 72.6; 63.4 62.8 63.4; 50.8 50.4 50.9; 43.5 43.0 43.5; 35.5 35.0 35.5; 28.5 28.3 28.5; 23.3 23.0 23.3; ]; %**each row consists of [Vab, Vbc, Vca], columns are different voltage levels
IFWL = [0.928 0.922 0.928; 0.854 0.852 0.861; 0.786 0.790 0.786; 0.720 0.718 0.718; 0.668 0.667 0.661; 0.580 0.560 0.567; 0.503 0.501 0.491; 0.410 0.405 0.400; 0.360 0.355 0.348; 0.305 0.301 0.296; 0.285 0.283 0.275; 0.283 0.279 0.276; ]; %**each row consists of [Ia, Ib, Ic], columns are different voltage levels
tFW = (10^-3)*[5.1 5.1 5.1; 5.1 5.1 5.1; 5.1 5.1 5.1; 5.1 5.1 5.1; 5.1 5.1 5.1; 5.1 5.1 5.1; 5.0 5.0 5.0; 4.8 4.8 4.8; 4.6 4.6 4.6; 4.6 4.6 4.6; 3.9 3.9 3.9; 3.6 3.6 3.6; ]; %**time between [(Vab, Ia) (Vbc, Ib) (Vca, Ic)]

VFWLN = VFWLL/sqrt(3); %Line-to-Neutral voltage (V)
VFW = sum(VFWLN,2)/3; %Average Line-to-Neutral voltage magnitude (V)
IFW = sum(IFWL,2)/3; %Average line current magnitude (V)
alphaFW = tFW/(1/60)*2*pi; %Angle measured between line-to-line voltage and line current (rad)
thetaFW = alphaFW - pi/6; % Angle between line-to-neutral voltage and line current (rad)
thetaFW = sum(thetaFW,2)/3; %Average angle between line-to-neutral voltage and line current (rad)
PFFW = cos(thetaFW); %Average power factor

%Calculate Input Power and Stator I^2R loss,
%Define Line-to-Neutral voltage phasor as reference (V)
VFW = VFW.*exp(i*0);

%Formulate Line Current phasor (A)
IFW = IFW.*exp(i*(-thetaFW));

%Formulate Input Power
SinFW = 3*(VFW .* conj(IFW)); %Three-phase (VA)
PinFW = real(SinFW); %Three-phase input power(W)


%Formulate Stator Loss
PsIR = 3*(abs(IFW).^2 .* R1); %Three-phase stator loss(W)

%Subtract stator loss from input power and plot resulting power vs. voltage
PM = PinFW - PsIR; %Total Power Tranferred across the air-gap

%Plot remaining power vs. voltage, and extend curve to 0 voltage intercept (this is friction and windage loss)
figure(1), clf;
plot(abs(VFW), PM,'x'); %Average line-to-neutral voltage vs. Air-gap power
hold on;

% p = polyfit(VFW,PM,2);
% x = 0:120;
% plot(0:120,polyval(p,x));
%axis([0 130 0 200]);

%**Analyze plot and estimate power at V=0
PlossFW = input('Analyze plot and estimate power at V=0. PFW = '); %Ask user to analyze plot
if isempty(PlossFW);
    PlossFW = 15; %**Enter zero-voltage intercept (this is friction and windage loss)
end;
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate Total Core Loss, Ph
%Friction and Windage loss has a value only during No-Load test
%Calculate Core Loss for No-Load test measurements
%Ph0 equation assumes rotor current is zero (s=0, Zr=infinite)
Ph0 = PinFW - PsIR - PlossFW; %Core loss for each voltage applied at 60Hz No-Load

%Calculate Core Conductance for each No-Load Friction Windage test
Gfe0 = Ph0./(m*VFW.^2) * (1+X1/Xm)^2; %Core Conductance in from No Load Test

Gfe0 = sum(Gfe0)/length(Gfe0); %Average core conductance from Locked Rotor Tests

%Calculate Core Resistance, Rfe
Rfe = 1/Gfe0; %Core Resistance


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate Rotor Resistance, R2
%Do Full Load Slip Test (Reduced Voltage Slip Test)
%Run motor at rated frequency, and rated voltage, and apply a load

%Measure stator voltage, line current, phase angle, slip, and follow other
%calculations from section 5.9.2.2.  
%Perform at several loads and take an average.
VSTLL = 2*[105 104 105]; %**Line-to-Line voltage during slip test (V)
ISTL = [1.10 1.13 1.14]; %**Line current magnitude during slip test (A)
tST = (10^-3)*[3.8 3.8 3.8]; %**Time between Line-to-Line voltge and line current peak (Vab to Ia, etc.)
ntST = [1728]; %**Rotor speed during Slip Test (rpm)

%Condition Data
VSTLN = VSTLL/sqrt(3); %Line-to-Neutral voltage during slip test (V)
VST = sum(VSTLN)/3; %Average Line-to-Neutral voltage during slip test (V)
IST = sum(ISTL)/3; %Average Line current during slip test (A)
alphaST = tST/(1/60)*2*pi; %Angle measured between line-to-line voltage and line current (rad)
thetaST = alphaST - pi/6; %Angle between line-to-neutral voltage and line current (rad)
thetaST = sum(thetaST)/3; %Average angle between line-to-neutral voltage and line current (rad)
PFST = cos(thetaST); %Average power factor

%Calculate Slip
SSST = ns - ntST; %Slip Speed (rpm)
sST = SSST/ns; %slip speed (pu)

%Correct R1 to specspecified tempeature before proceding
R1 =  R1*(tb+k1)/(ta+k1); %Correct value for temperature (ohm)% Stator resistance

%Calculate V2 (Rotor Winding Voltage reflected to stator) using eqt. 44
V1 = VST; %Sub Value for equation taken from std. 112
I1 = IST; %Sub Value for equation taken from std. 112
theta1 = thetaST; %Sub Value for equation taken from std. 112
V2 = sqrt((V1-I1*(R1*cos(theta1) - X1*sin(theta1)))^2 + (I1*(R1*sin(theta1) + X1*cos(theta1)))^2);

%Calculate Theta2 using eqt. 45
theta2 = atan((-I1*(R1*sin(theta1) - X1*cos(theta1)))/(V1 - I1*(R1*cos(theta1) - X1*sin(theta1))));

%Calculate Ife and Ie using Eqation 49 and 46
Im = V2/Xm; %Equations gives Ie, but figure 2 shows V2/Xm = Im
Ife = V2/Rfe ;

%Calculate I2 using Equation 50
I2 = sqrt((I1*cos(theta1) - Im*sin(theta2) - Ife*cos(theta2))^2 + (-I1*sin(theta1) + Im*cos(theta2) + Ife*sin(theta2))^2);

%Calculate Z2 using Equation 53
Z2 = V2/I2;

%Calculate R2 using equation 54
R2 = sST*sqrt(Z2^2 - X2^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Final Parameter Values (Corrected for Temperature, eqt. 3)
display('Machine Parameters:')

L1 % Stator Leakage Inductance (H)
R1 % Stator Resistance (ohm)
R2 % Rotor Resistance (ohm)
L2 % Rotor Leakage Inductance (H)
Lm % Magnetizing Inductance (H)
Rfe %Core Resistance (ohm)
PlossFW %Friction and Windage Loss (W)

%END SCRIPT