%Ejercicio Capítulo 6 del texto de Carson W. Taylor
%
%Se obtienen las curvas Qinyectada - Voltaje en la barra REC del sistema de la Fig. 6.1 (Pag 140)


%Programa Principal
function simula
clear all
close all
clc


%Datos Nominales de la carga:
Pm = 0.6;            % Potencia activa del motor en pu.
Qm = 0.4;            % Potencia reactiva del motor en pu.
Pt = 0.6;            % Potencia activa de la carga termostatica en pu.
Bc = 0.324;          % Susceptancia del banco de condensadores

%Datos del transformador
Xt = 0.0833;         % Reactancia del transformador en pu.

% CASO 1: Sistema operando en condiciones normales : Prefalla

 %Datos del sistema, equivalente de thevenin visto desde el nudo REC
  Vth=1.0604;        
  Xth=0.2674;

 %Modelo Cargas:
 %     Motor              :  Pm = cte   Qm = cte
 %     Carga Termostática :  Pt = cte
 %Transformador           :  Las tomas actuan para mantener Vl =  1 pu.
 
  %Caso 1a: Curva Qin - Vrec en el caso que la toma alcanza su límite inferior toma=0.9
  k = 0;
  toma = 0.9;                            % Límite inferior de la toma
  for Vl=0.8:0.001:1                     % Variación de tensión en la carga
      Pmot = Modelo_P(0,Pm,Vl);          % Modelo motor a Pm=cte
      Pter = Modelo_P(0,Pt,Vl);          % Modelo carga termostatica a Pt=cte
      Ps   = Pmot + Pter;                % Potencia activa total
      Qmot = Q_motor(0,Qm,0,0,0,Vl);     % Modelo motor a Qm=cte
      Qcon = Q_banco(Bc,Vl);             % Reactiva del banco en la barra de carga
      Icar = I_carga(Ps,Qmot,Qcon,Vl);   % Corriente consumida por la carga
      Qtra = Q_trafo(Xt,Icar);           % Potencia reactiva absorvida por el transformador
      E2tr = Esec_trafo(Xt,Icar,Vl);     % Tensión secundaria en los enrrollados del transformador
      k=k+1;
      Vrec(k) = toma*abs(E2tr);          % Tensión en la barra REC (Primario del transformador)
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);  %Reactiva total inyectada en  REC
   end

  %Caso 1b: Curva Qin - Vrec en el caso que las tomas varian entre 0.9 y 1.1 para mantener Vl = 1 pu.
  %Note que en el caso 1a la tension final es 1.0 pu. que corresponde al valor que las tomas variando 
  %entre 0.9 y 1.1 deben mantener constante. Por lo tanto las potencias consumidas en la carga
  %son las calculadas en la ultima iteración del ciclo anterior pues en este caso Vl=1 pu. Por lo tanto
  %no es necesario re-calcularlas.
  for toma=0.9:0.001:1.1;
     k = k + 1;
     Vrec(k)=toma*abs(E2tr);
     Qin(k) = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
   
  %Caso 1c: Curva Qin - Vrec en el caso que la toma alcanza su límite superior toma=1.1
  toma=1.1;
  for Vl=1.0:0.001:1.2
      Pmot = Modelo_P(0,Pm,Vl);
      Pter = Modelo_P(0,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(0,Qm,0,0,0,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
  
   %Grafica resultados Caso 1
   plot(Vrec,Qin*500)
   grid on
   axis([ 0.4 1.2 -100 400])
   hold on

% CASO 2: Sistema pierde una linea, se analiza estado inmediatamente despues de la falla

   %Datos del sistema en condiciones de falla
    Vth=1.0604;
    Xth=2*0.2674;

    %Modelo Cargas:
    %     Motor              :  Zm = cte   
    %     Carga Termostática :  Zt = cte
    %Trafo                   :  No tiene tiempo de actuar, la toma queda en el valor de pre-falla
      E1tr = 1.05;         % Tension de pre-falla en el primario del trafo
      Vl   =1;             % Tension de pre-falla en la carga
      % Cálculo de la toma en condiciones de pre-falla
        Pmot = Modelo_P(1,Pm,Vl);
        Pter = Modelo_P(1,Pt,Vl);
        Ps   = Pmot + Pter;
        Qmot = Q_motor(1,Qm,0,0,0,Vl);
        Qcon = Q_banco(Bc,Vl);
        Icar = I_carga(Ps,Qmot,Qcon,Vl);
        Qtra = Q_trafo(Xt,Icar);
        E2tr = Esec_trafo(Xt,Icar,Vl);
      toma = abs(E1tr/E2tr);          % Esta toma permanece constante
   k=0;
   for Vl=0.7:0.001:1.2
      Pmot = Modelo_P(1,Pm,Vl);
      Pter = Modelo_P(1,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(1,Qm,0,0,0,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
    plot(Vrec,Qin*500)
 
  % CASO 3: Sistema pierde una linea, se analiza estado instantes despues de la falla

   %Datos del sistema en condiciones de falla
    Vth=1.0604;
    Xth=2*0.2674;

    %Modelo Cargas:
    %     Motor              :  Pm = cte   Qm = Qm(Vl) Van Cutsem
    %     Carga Termostática :  Zt = cte
    %Trafo                   :  la toma aun no alcanza actuar y queda en el valor del caso 2
    
    %Parametros del motor:
    xm = 3.06;
    x  = 0.2;
  k=0;
   for Vl=0.55:0.001:1.2
      Pmot = Modelo_P(0,Pm,Vl);
      Pter = Modelo_P(1,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(2,Qm,xm,x,Pmot,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
    plot(Vrec,Qin*500)

% CASO 4: Sistema pierde una linea, se analiza estado segundos después de la falla

%Modelo Cargas:
%     Motor              :  Pm = cte   Qm = Qm(Vl) Van Cutsem
%     Carga Termostática :  Zt = cte
%Trafo                   :  Reaccionan las tomas las que tratan de mantener Vl= 1pu.
 %Caso 4a: Curva Qin - Vrec en el caso que la toma alcanza su límite inferior toma=0.9
 k=0;
 toma=0.9;
 for Vl=0.6:0.001:1;
      Pmot = Modelo_P(0,Pm,Vl);
      Pter = Modelo_P(1,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(2,Qm,xm,x,Pmot,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end

  %Caso 4b: Curva Qin - Vrec en el caso que las tomas varian entre 0.9 y 1.1 para mantener Vl = 1 pu.
  %Note que en el caso 1a la tension final es 1.0 pu. que corresponde al valor que las tomas variando 
  %entre 0.9 y 1.1 deben mantener constante. Por lo tanto las potencias consumidas en la carga
  %son las calculadas en la ultima iteración del ciclo anterior pues en este caso Vl=1 pu. Por lo tanto
  %no es necesario re-calcularlas.
  for toma=0.9:0.001:1.1;
     k = k + 1;
     Vrec(k)=toma*abs(E2tr);
     Qin(k) = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
   
  %Caso 4c: Curva Qin - Vrec en el caso que la toma alcanza su límite superior toma=1.1
  toma=1.1;
  for Vl=1.0:0.001:1.2;
      Pmot = Modelo_P(0,Pm,Vl);
      Pter = Modelo_P(1,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(2,Qm,xm,x,Pmot,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
  
   %Grafica resultados Caso 4
   plot(Vrec,Qin*500,'r')
   
% CASO 5: Sistema pierde una linea, se analiza estado en regimen permanente de falla

%Modelo Cargas:
%     Motor              :  Pm = cte   Qm = Qm(Vl) Van Cutsem
%     Carga Termostática :  Pt = cte
%Trafo                   :  Las tomas  tratan de mantener Vl= 1pu.
 %Caso 5a: Curva Qin - Vrec en el caso que la toma alcanza su límite inferior toma=0.9
 k=0;
 toma=0.9;
 for Vl=0.7:0.001:1;
      Pmot = Modelo_P(0,Pm,Vl);
      Pter = Modelo_P(0,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(2,Qm,xm,x,Pmot,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end

  %Caso 5b: Curva Qin - Vrec en el caso que las tomas varian entre 0.9 y 1.1 para mantener Vl = 1 pu.
  %Note que en el caso 1a la tension final es 1.0 pu. que corresponde al valor que las tomas variando 
  %entre 0.9 y 1.1 deben mantener constante. Por lo tanto las potencias consumidas en la carga
  %son las calculadas en la ultima iteración del ciclo anterior pues en este caso Vl=1 pu. Por lo tanto
  %no es necesario re-calcularlas.
  for toma=0.9:0.001:1.1;
     k = k + 1;
     Vrec(k)=toma*abs(E2tr);
     Qin(k) = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
   
  %Caso 5c: Curva Qin - Vrec en el caso que la toma alcanza su límite superior toma=1.1
  toma=1.1;
  for Vl=1.0:0.001:1.2;
      Pmot = Modelo_P(0,Pm,Vl);
      Pter = Modelo_P(0,Pt,Vl);
      Ps   = Pmot + Pter;
      Qmot = Q_motor(2,Qm,xm,x,Pmot,Vl);
      Qcon = Q_banco(Bc,Vl);
      Icar = I_carga(Ps,Qmot,Qcon,Vl);
      Qtra = Q_trafo(Xt,Icar);
      E2tr = Esec_trafo(Xt,Icar,Vl);
      k=k+1;
      Vrec(k) = toma*abs(E2tr);
      Qin(k)  = Qtra + Qmot - Qcon - Q_sistema(Vth,Xth,Vrec(k),Ps);
   end
  
   %Grafica resultados Caso 5
   plot(Vrec,Qin*500,'g')
  
  
  %Curvas Q-Vrec de los bancos de condensadores en la barra REC
  B1=(153 / 500 )     ;   %Banco de condensadores de 153 MVAR en la Barra REC
  B2=(153 +  75 )/500 ;   %Se agregan 75 MVAR al banco
  B3=(153 + 125 )/500 ;   %Se agregan 125 MVAR al banco
  Vrec=0:0.01:1.2;
  Qbanco1=B1.*Vrec.^2;
  Qbanco2=B2.*Vrec.^2;
  Qbanco3=B3.*Vrec.^2;
  plot(Vrec,Qbanco1*500,'k',Vrec,Qbanco2*500,'k',Vrec,Qbanco3*500,'k')
 
%Fin Programa principal



%**********************   DEFINICION DE SUBRUTINAS   ******************************

%Función para obtener la potencia activa consumida por una carga considerando diferentes modelos

function Pactiva = Modelo_P(tipo,cte1,voltaje)
%tipo: modelo de potencia activa de la carga   0: Potencia cte  1: Z=cte  2: otros
%cte1 : Constante de entrada para definir : Po si tipo=0, Go si tipo=1
%Voltaje: Tensión en bornes del elemento
switch tipo
case 0
   Pactiva=cte1;
case 1
   Pactiva=cte1*voltaje^2;
case 2
   %otro modelo puede agregarse aqui
otherwise
   disp('Error tipo de modelo de carga activa no existe'); pause; quit
end

%Función Para obtener la potencia reactiva suministrada por un banco de condensadores

function Qc= Q_banco(Bc,voltaje)
%Bc : Susceptancia del banco de condensadores
%Voltaje: Tensión en bornes del banco
if Bc < 0
   disp('Error Bc debe ser mayor que cero'); pause; quit
end
Qc=Bc*voltaje^2;

%función para evaluar la potencia reactiva consumida por un motor 

function Qm = Q_motor(tipo,cte1,xm,x,pm,voltaje)
%tipo: modelo de potencia reactiva del motor  0: Potencia cte  1: Z=cte  2: Van Cutsen 3:otros
%cte1 : Constante de entrada para definir : Qo si tipo=0 , Bo si tipo=1
%Voltaje: Tensión en bornes del motor
switch tipo
case 0
   Qm = cte1;
case 1
   Qm = cte1*voltaje^2;
case 2
   Qm = voltaje^2/xm + voltaje^2/(2*x) - sqrt( ( voltaje^2/(2*x) )^2 - pm^2 );
otherwise
   disp('Error tipo de modelo de potencia reactiva del motor no existe'); pause; quit
end

%función para evaluar la potencia reactiva consumida por el transformador

function Iload= I_carga(Psis,Qmot,Qban,Vl)
%Pmot   : Potencia consumida por el motor
%Pter   : Potencia consumida por la carga termostatica
%Qmot   : Reactiva consumida por el motor
%Qban   : Reactiva aportada por el banco de condensadores
%Vl     : Tensión en bornes de la carga 
Iload = ( Psis -j*(Qmot - Qban) ) / Vl;

function Qt = Q_trafo(xtrafo,Il)
Qt = xtrafo*abs(Il)^2;

%función para evaluar la potencia reactiva aportada por el sistema

function Qsis = Q_sistema(Vth,Xth,Vrec,Ps)
Qsis = sqrt( (Vth*Vrec/Xth)^2 - Ps^2 ) -Vrec^2/Xth;

%Función para evaluar la tension en los enrrollados secundarios del transformador

function E2t = Esec_trafo(Xt,Icar,Vl)
E2t = j*Xt*Icar + Vl;


   
   
   
