%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fichier correspondant au TD 3 de Machines Electriques (ENSEM 2A NRJ S7) :
%  - modification du fichier initial fourni par N. Takorabet 
%  - fonctionne sous Matlab ou GNU/Octave
%
%     v3.1 : J. Fontcchastagner / nov. 2022
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear 
close all
clc


rs = 1.25;     % resistance d'une phase statorique a 20 degres C
% test pour savoir si on est sous Matlab ou Octave 
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if (~isOctave)
  j=1i; % pour eviter warning dans l'editeur de Matlab...
end 

%%%====================================
%%% Essai a vide separation des pertes
%%%=====================================  
Ess_vide = [  404    380     360    340    320   300     280    260    240    220   200      180     100      %U0 [V]
             7.26   5.81    4.92   4.42   3.93  3.58    3.26   2.95   2.68   2.47  2.16     2.0     1.1    %I0 [A]
             0.68   0.56    0.48   0.42   0.38  0.33    0.29   0.27   0.23   0.20  0.17     0.15   0.08 ];  %P0 [kW]   

U0 = Ess_vide(1,:);      % vecteur contenant les tensions entre-phases mesurees
V0 = U0/sqrt(3);         % vecteur contenant les tensions phase-neutre 
I0 = Ess_vide(2,:);      % vecteur contenant les courant statoriques mesures
P0 = Ess_vide(3,:)*1e3;  % vecteur contenant les Puissances actives mesurees


 %Calculer : 
   % Le vecteur S0 : Puissance apparante
   % Le vecteur Q0 : Puissance reactive
   % Le vecteur cosfi_0 : facteur de puissance a vide 
   % Le vecteur sinfi_0 : sinus phi a vide
   
S0 = sqrt(3)*U0.*I0;
Q0 = sqrt(S0.^2 - P0.^2);
cosphi_0 = P0./S0;
sinphi_0 = -sqrt(1.-cosphi_0.^2);


% Calculer le tableau courant a vide en complexe sous la forme : 
% I0_cmplx=I0_re-1i*I0_re

I0_cmplx = I0.*(cosphi_0 + j*sinphi_0);


% Calculer : 
   % Le tableau E0 : Fem
   % le tableau PfpPm=P0-3*rs*I0.*I0 qui represente les pertes fer + les
   % pertes mecaniques 

E0_cmplx = V0 - rs*I0_cmplx;
E0 = abs(E0_cmplx);
E02 = E0.^2;
PfpPm = P0-3*rs*I0.^2;


%% interpolation PfpPm au 1er ordre :PfpPm =alfa x V^2 + beta
%% avec alfa = 3/Rf   et beta = Pm
% Utiliser les fonctions polyfit et polyval de Matlab
% Extraire les Pertes mecaniques : P_meca = beta
% Extraire Rf= 3/alfa

cp = polyfit(E02,PfpPm,1);
alpha = cp(1);
P_meca = cp(2);
Rf = 3/alpha;


% Calculer la Reactance magnetisante  en utilisant le dernier point a vide 404V
Xm = 3*E0(1)^2/Q0(1);


% Tracer Les points de mesures Po-3*Rs*IO ainsi le polynome d'approximation : 
% alfa x V^2 + beta

figure(1)
plot(E02,PfpPm,'x','linewidth',5);
hold on
plot([E02 0],[alpha*E02+P_meca P_meca],'linewidth',3);
grid on
title('Essai a vide (separation des pertes)')
xlabel('E^2 (V^2)')
ylabel('Pf + Pm  (W)')


%%%==========================
%%%   Essai a rotor bloque
%%%==========================

Vcc=70/sqrt(3);
Icc=10;
Pcc=500;

%%% Calculer 
%% Scc Puissance apparante
%% Qcc Puissance reactive
%%% cosfi_cc et sinfi_cc 
% Icc_complx sous la forme Icc*(cosfi_cc-1*i*sinfi_cc)
% Ecc_cmplx 

Scc = 3*Vcc*Icc;
Qcc = sqrt(Scc^2-Pcc^2);
cosphi_cc = Pcc/Scc;
sinphi_cc = Qcc/Scc;
Icc_cmplx = Icc*(cosphi_cc - j*sinphi_cc);
Ecc_cmplx = Vcc - rs*Icc_cmplx;


%% Pfcc et Qfcc les puissances actives et reactives
%  consommees dans la branche magnetisante lors de cet essai

Pfcc = 3*abs(Ecc_cmplx)^2/Rf;
Qfcc = 3*abs(Ecc_cmplx)^2/Xm;


% Calculer Rp2 et Np2_omega avec les deux 
% hypotheses (negliger ou pas la branche magnetisante)

Iprim2_cmplx = Icc_cmplx - Ecc_cmplx/Rf - Ecc_cmplx/(j*Xm);
Pscc = 3*rs*abs(Icc_cmplx)^2;
Rprim2 = (Pcc - Pfcc - Pscc)/3/abs(Iprim2_cmplx)^2;
Nprim2w = (Qcc - Qfcc)/3/abs(Iprim2_cmplx)^2;


disp(' ')
disp(' Parametres du schema equivalent :')
disp('-----------------------------------')
disp(' ')
disp([' R1 = ',num2str(rs),' ohms'])
disp([' Rf = ',num2str(Rf),' ohms'])
disp([' L1w = ',num2str(Xm),' ohms'])
disp([' N2''w = ',num2str(Nprim2w),' ohms'])
disp([' R2'' = ',num2str(Rprim2),' ohms'])
disp(' ')
disp([' Pertes mecaniques = ',num2str(P_meca),' W'])
disp(' ')


% En négligeant :
Rprim2app = (Pcc - Pscc)/3/Icc^2;
Nprim2wapp = Qcc/3/Icc^2;

Ecart_R2prim = (Rprim2-Rprim2app)/Rprim2*100;
Ecart_Nprim2w = (Nprim2w-Nprim2wapp)/Nprim2w*100;

disp(['-> Ecart sur R2'' en negligeant la branche magnetisante : ', ... 
    num2str(Ecart_R2prim), ' %' ])
disp(['-> Ecart sur N2''w en negligeant la branche magnetisante : ', ... 
    num2str(Ecart_Nprim2w), ' %' ])
disp(' ')


%%% A ce niveau, les parametres du schema equivalent sont connus.
%%----------------------------------------------------------------

%%% Question 2-d : TRacer la caracteristique couple-vitesse dans les deux
%%% hypotheses
p = 2;
Ns = 1500;          % vitesse synchrone (tr.min^{-1})
Omegas = Ns*pi/30;  % vitesse synchrone (rad.s^{-1})
ws = p*Omegas;      % pulsation courants statoriques

% Plutôt qu'une boucle, on travaille directement sur un vecteur contenant
% toutes les valeurs de glissement :
npoints = 1000;         % nombre de point pour le trace
glisdeb = 1e-12;        % on ne peut pas commencer à 0 (division par g)
glissement = linspace(glisdeb,1,npoints);
vit = (1-glissement)*Ns;  % vitesse (pour le trace)

Rprim2surg = Rprim2./glissement;

Zp = Rf*j*Xm/(Rf+j*Xm);
Zs = j*Nprim2w + Rprim2surg;
Zeq = rs + Zp.*Zs./(Zp+Zs);

V1 = 230;
I1 = V1./Zeq;
E1 = V1-rs.*I1;
Ip = E1./Zp;
I2 = I1-Ip;
Pt = 3*Rprim2surg.*abs(I2).^2;
Gamma = Pt/Omegas;

% caracteristique ideale :
Gamma_max = 3*p*V1^2/(ws*Nprim2w)/2;
gM = Rprim2/Nprim2w;
Gamma_ideal = 2*Gamma_max./(gM./glissement+glissement/gM);

figure(2)
plot(vit,Gamma,'linewidth',3);
hold on
plot(vit,Gamma_ideal,'linewidth',3);
grid on
xlim([0 1.05*vit(1)])
title('Couple - vitesse')
xlabel('Vitesse (tr.min^{-1})')
ylabel('Couple  (N.m)')
legend('Caracteristique Reelle','Caracteristique Ideale',...
        'Location','NorthWest')


% question 4-b :
%-----------------

Radd = 0.5;
Rprim2surg = (Rprim2+Radd)./glissement;

Zp = Rf*j*Xm/(Rf+j*Xm);
Zs = j*Nprim2w + Rprim2surg;
Zeq = rs + Zp.*Zs./(Zp+Zs);

V1 = 230;
I1 = V1./Zeq;
E1 = V1-rs.*I1;
Ip = E1./Zp;
I2 = I1-Ip;
Pt = 3*Rprim2surg.*abs(I2).^2;

Gamma2 = p*Pt/ws;

figure(3)
plot(vit,Gamma,'linewidth',3);
hold on
plot(vit,Gamma2,'linewidth',3)
grid on
xlim([0 1.05*vit(1)])
title('Couple - vitesse')
xlabel('Vitesse (tr.min^{-1})')
ylabel('Couple  (N.m)')
legend('R2','R2 + Radd','Location','northwest')


% question 4-c : 
%-----------------
Gammar = 40;   % couple de charge 

% on rajoute Gamma_r sur les plots precedents
plot([0 1.05*vit(1)],[Gammar Gammar],'g','linewidth',2);
legend('R2','R2 + Radd','\Gamma_r','Location','northwest')
figure(2)
plot([0 1.05*vit(1)],[Gammar Gammar],'g','linewidth',2)
legend('Caracteristique Reelle','Caracteristique Ideale',...
        '\Gamma_r','Location','NorthWest')

% Cherchons le glissement correspondant (celui de faible valeur, 
%                                            i.e. de la zone stable) :
fonction_a_annuler = @(x) interp1(glissement,Gamma,x)- Gammar;
gdepart = 1e-3;  % point de depart
if (~isOctave) 
  % On est sous Matlab
  options = optimoptions('fsolve');
  options.Display = 'off';
  g = fsolve(fonction_a_annuler,gdepart,options);
else
  % On est sous Octave
  pkg load optim; % attention il faut avoir le pkg optim. A installer par :
                  % pkg install -forge struct statistics optim
  g = fsolve(fonction_a_annuler,gdepart);
end 

Rprim2surg = Rprim2/g;
Zp = Rf*j*Xm/(Rf+j*Xm);
Zs = j*Nprim2w + Rprim2surg;
Zeq = rs + Zp*Zs/(Zp+Zs);
V1 = 400/sqrt(3);
I1 = V1/Zeq;
E1 = V1-rs*I1;
Ip = E1/Zp;
I2 = I1-Ip;
Pt = 3*Rprim2surg*abs(I2)^2;
Gamma = Pt/Omegas; % pour verification, normalement : Gamma = Gammar
Omega = (1-g)*Omegas;
Pu = Gamma*Omega;
Ps = real(3*V1*conj(I1));
eta = Pu/Ps;

disp(' ')
disp(' Point de fonctionnement :')
disp('---------------------------')
disp(' ')
disp([' Couple = ',num2str(Gamma),' N.m'])
disp([' Vitesse = ',num2str(Omega*30/pi),' rpm'])
disp([' Glissement = ',num2str(g*100),' %'])
disp(' ')
disp([' Tension (entre-phases) = ',num2str(V1*sqrt(3)),' V'])
disp([' Courant = ',num2str(abs(I1)),' A'])
disp([' Puissance stator = ',num2str(Ps),' W'])
disp(' ')
disp([' Rendement = ',num2str(eta*100),' %'])
disp(' ')




