Estabilidad transitoria en sistemas multimaquinas

Estabilidad Transitoria en Sistemas Multimaquinas Problema Estabilidad Transitoria: “Habilidad de las máquinas sincrónicas interconectadas de operar en sincronismo” ¿Que implica operar en sincronismo? G1 G2 Sistema eléctrico de potencia interconectado. Maq. referencia Eje rotatorio sincrónico de G1 (eje de referencia) Eje rotatorio sincrónico de G2 ws (Gp:) ws ws Cte. En régimen estacionario. =cte=0 (en régimen estacionario) En condiciones estacionarias las posiciones relativas de los rotores permanecen constantes y corresponde a la transferencia de potencia entre las máquinas y la red el sistema está en sincronismo.

¿Bajo que circunstancias el sistema puede perder el sincronismo? Cuando aparece una alteración abrupta en el régimen de transferencia de potencia, las posiciones relativas de los rotores se verán alteradas, pudiendo el sistema de ser capaz o no de autorestituirse a un nuevo estado de equilibrio. Ejemplo de alteración abrupta en el régimen de transferencia de potencia: Corto circuito y las maniobras sucesivas Dado el siguiente sistema eléctrico de potencia interconectado: G G G Slack Gen_2 Gen_3 Carga_5 Carga_6 Carga_4 c.c Se estudiará la variación de los ángulos de las máquinas tras la ocurrencia de un cortocircuito en la líneas Carga_5-Carga_6 próximo a la barra Carga_6, seguido de la apertura de dicha línea (eliminación de la falta). (Gp:) ws )) )) Eje rotórico en régimen oscilatorio Eje rotatorio sincrónico de G1 (eje de referencia)

Se grafica la variación de los ángulos de las diferentes máquinas respecto a la variación del ángulo de la máquina de referencia, para dos situaciones de duración del cortocircuito previo a la apertura de la línea Caso 1: tiempo c.c. 0.4s Caso 2: tiempo c.c. 0.5s régimen estacionario c.c 0.4s Apertura de la línea (eliminación de la falta) régimen estacionario c.c 0.5s Apertura de la línea (eliminación de la falta) Caso 1, caso estable, el sistema fue capaz de restituirse a una nueva situación de equilibrio. Caso 2, representa los casos de inestabilidad, donde una máquina pierde el sincronismo. Al estudio de la capacidad de una red interconectada de restituirse frente a una sucesión de modificaciones abruptas en la transferencia de potencia en función de: – los tiempos, tipos y magnitudes de dichas modificaciones – parámetros de la red lo llamamos: Estabilidad Transitoria en Sistemas Multimaquinas

Objetivo Estudio del comportamiento de un sistema de potencia multimáquina interconectado del punto de vista de la estabilidad transitoria (oscilaciones electromecánicas). Esto se logrará a través del desarrollo de una herramienta computacional que contemple todas las posibles maniobras que se suceden luego de una falta (intervalos de estudio), como también que tenga una razonable flexibilidad como para realizar estudios paramétricos. Modelo Dinámico Simplificado de la Máquina Sincrónica Comportamiento Mecánico Comportamiento Eléctrico (Gp:) ws )) )) Eje rotórico en régimen oscilatorio Eje rotatorio sincrónico de G1 (eje de referencia) La ley de variación de esta dada por la llamada ecuación de oscilación. En régimen estacionario En régimen oscilatorio asumimos Pm permanece constante durante todo el estudio como venía del régimen estacionario. i i’ G Ra Xd’ barra de la red barra interna (ficticia) P’ mi Pi , Qi P’ e i’ , Q’ ei’ Ii Asumimos:

1 – Cada máquina síncrona es representada por una fuente de tensión de módulo constante atrás de una impedancia (resistencia de armadura más reactancia transitoria directa), conforme figura arriba.) 2 – La tensión detrás de la reactancia transitoria se considera constante durante todo el intervalo de estudio y el coincide con el ángulo mecánico del rotor. 3 – Se asume que la potencia mecánica de entrada de la máquina es constante durante todo el periodo de la simulación, no son considerados acciones de reguladores. 4 – Se desprecia potencia disipada en los arrollamientos amortiguadores. 5 – Usando las tensiones pre-falta, todas las cargas son convertidas en admitancias a tierra y se asumen constantes durante toda la simulación. Hipótesis simplificatorias 1) Mediante flujo de carga se determina para cada máquina: – potencia mecánica – tensiones en bornes – en régimen para cada máquina. 2) Para cada intervalo de estudio (falta, apertura línea, reenganche): – Determino la matriz admitancia para estudios de estabilidad – Resuelvo el sistema de ecuaciones diferenciales (son tantas ecuaciones como máquinas): Metodología para resolver el problema de estabilidad transitoria H, f, |E|, Pm : constantes

Potencia mecánica y tensión en bornes de las máquinas i i’ G Ra Xd’ barra de la red barra interna ficticia P’mi Pi , Qi P’ei’ , Q’ei’ Las magnitudes eléctricas (potencia y tensión) en la barra i, son conocidas como resultado de la corrida del flujo de carga. La corriente Ii entonces está dada por: Ii Siendo la tensión en bornes del generador: Entonces la potencia activa eléctrica de la máquina: En régimen tenemos un equilibro entre la potencia mecánica de entrada de la máquina y la potencia activa eléctrica en bornes de la misma: Este valor de potencia mecánica una vez determinado lo asumimos constante durante toda la simulación (hipótesis simplificatoria nro. 3), asimismo, el módulo de la tensión en bornes también permanece constante (hipótesis simplificatoria nro. 2). cte. durante todo el estudio Condición inicial del siguiente intervalo de estudio (falta). cte. durante todo el estudio

Potencia Eléctrica G G G Slack Gen_2 Gen_3 Carga_5 Carga_6 Carga_4 Barra de referencia G G G Slack Gen_2 Gen_3 Carga_5 Carga_6 Carga_4 y60 y40 y42 y14 y46 y16 y56 y35 y50 y’G1 G’1 G’3 y’G3 Y’G3 G’2 y14 y16 y56 y46 y15 y15 Dada la siguiente configuración de red: El modelado para estudios de estabilidad es el siguiente:

Diferencias respecto al modelado para estudios de flujo de carga – Aparecen las reactancias transitorias de las máquinas y una barra interna (ficticia) detrás detrás de las reactancias. – Todas las cargas son convertidas en admitancias a tierras, calculadas en base a las tensiones pre-falta, se asume que las admitancias permanecen constantes durante todo el intervalo de estudio (hipótesis simplificatoria nro. 5). En el modelado para estudios de transitorios, a todas las barras de la red (excepto las ficticias) concurren solamente elementos pasivos, entonces: G Slack y14 y16 y15 Slack y14 y16 y’G1 G’1 y15 G Modelo flujo de carga Modelo estabilidad transitoria Analogamente para las barras de carga La ecuación nodal para una red de n barras y ng máquinas queda: Corrientes en las barras internas maquinas Corrientes en las barras =0

Siendo entonces cero las corrientes en las barras podemos re-escribir la ecuación nodal: El vector Vn puede ser eliminado por las siguientes substituciones: Sustituyendo Vn en la segunda ecuación: Donde tiene dimensión igual al número de generadores. La potencia eléctrica de salida de cada máquina puede ahora expresarse en términos de su tensión interna: Reescribiendo la tensión y las componentes de la matriz admitancia en la forma polar:

Representación de la red en condición de falta y maniobras sucesivas (intervalos de estudio) Se pueden estudiar tres situaciones: 1 – falta y apertura definitiva 2 – falta, apertura y reenganche (falta extinguida). 3 – falta, apertura, reenganche (falta permanece) y apertura definitiva Red con la falta Red con la falta Red con la falta Red sin la línea que limpia la falta Red restablecida, falta extinguida Red sin la línea que limpia la falta Red sin la línea que limpia la falta Red con la falta Red sin la línea que limpia la falta En cada intervalo se calculará la potencia eléctrica de cada máquina utilizando la Yred que corresponda: – red original. – red con una barra en falta trifásica a tierra. – red sin la línea que limpia la falta. Los estudios se hacen para faltas trifásicas próximas a la barra. Por ejemplo durante una falta próxima a la barra Carga_6, la configuración de la red queda: (Gp:) Barra de referencia (Gp:) G (Gp:) G (Gp:) G (Gp:) Slack (Gp:) Gen_2 (Gp:) Gen_3 (Gp:) Carga_5 (Gp:) Carga_4 (Gp:) y40 (Gp:) y42 (Gp:) y14 (Gp:) y16 (Gp:) y56 (Gp:) y35 (Gp:) y50 (Gp:) y’G1 (Gp:) G’1 (Gp:) G’3 (Gp:) y’G3 (Gp:) Y’G3 (Gp:) G’2 (Gp:) y46 (Gp:) y15 Eliminando fila y columna de la Yred original obtengo Ycc, luego calculo la matriz reducida

Suponiendo que sea la apertura de la línea Carga_5-Carga_6 la que elimina la falta, durante el intervalo que la línea esta abierta la configuración de la red queda: Barra de referencia G G G Slack Gen_2 Gen_3 Carga_5 Carga_6 Carga_4 y60 y40 y42 y14 y16 y35 y50 y’G1 G’1 G’3 y’G3 Y’G3 G’2 y46 y15 Elimino de la matriz original Yred los aportes de admitancia y suceptancia correspondientes a la línea obteniendo Ycl, luego calculo la matriz reducida Resolución de la ecuación de oscilación La ecuación de oscilación de una maquina i está dada por: Donde H es la constante de inercia de la máquina expresada en una base MVA común SB. Si HGi es la constante de inercia de la máquina i expresada en base a la potencia SGi de la misma entonces Hi esta dado por: La misma se resuelve por métodos numéricos de resolución de ecuaciones diferenciales, para lo que hay que representarla en la forma de variables de estado, esto implica representar una ecuación diferencial de orden n en n ecuaciones diferenciales de orden 1 mediante convenientes cambios de variable. En este caso siendo la ecuación de oscilación de segundo orden: Tendremos entonces un sistema de 2xNro._de_máquinas ecuaciones diferenciales de orden uno.

Software desarrollado Archivo datos de la red G G G Slack Gen_2 Gen_3 Carga_5 Carga_6 Carga_4 % DATOS DE BARRA % CARGA GENERACION min max Shunt Shunt % BARRA TENSION MW MVAr MW MVAR MVAr MVAr MVAr Suceptancia SL Slack 1.06 0 0 0 0 0 0 0 0 PV Gen_2 1.04 0 0 150 0 0 140 0 0 PV Gen_3 1.03 0 0 100 0 0 90 0 0 PQ Carga_4 1 100 70 0 0 0 0 0 0 PQ Carga_5 1 90 30 0 0 0 0 0 0 PQ Carga_6 1 160 110 0 0 0 0 0 0 % % % DATOS DE LINEAS % BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIA Linea Slack Carga_4 0.035 0.225 0.013 Linea Slack Carga_5 0.025 0.105 0.009 Linea Slack Carga_6 0.040 0.215 0.011 Linea Carga_4 Carga_6 0.028 0.125 0.007 Linea Carga_5 Carga_6 0.026 0.175 0.06 % % % DATOS DE TRANSFORMADORES % BARRA_1 BARRA_2 RESISTENCIA REACTANCIA TAP Trafo Gen_2 Carga_4 0.00 0.035 1 Trafo Gen_3 Carga_5 0.00 0.042 1 % % DATOS DE LOS GENERADORES (SOLO PARA CALCULO DE CC y Est.) % BARRA_1 Ra X' H Gen Slack 0.00 0.20 20 Gen Gen_2 0.00 0.15 4 Gen Gen_3 0.00 0.25 5

Ejecución del programa Se ejecuta dalestabil.m, es un script desarrollado para que quede “cómodo” el ingreso de los datos específicos para el estudio de estabilidad y procesamiento de los resultados: % funcion para corrida del programa estabilidad estabil.m. % R. Hirsch, Junio 2002 global Sb f f=60; % Frecuencia. Sb=100; % Potencia base. archivo='pag516.m'; % Archivo de datos de la red barracc='Carga_6'; % Barra en falta barrab1='Carga_6'; % Extremo 1 de la linea que limpia la falta barrab2='Carga_5'; % Idem extremo 2 % tadata=(tc1,tr,tc2,tf) vector tiempos de los sucesos, arranca en 0s. % tc1 : Clearing time 1, la falta se mantiene de 0 a tc1. % tr : Reenganche (opcional) desde tc1 la red esta sin la linea en falta % a partir de tr reengancho la linea (*) % tc2 : Si se especifica tc2 se asume que la falta se sustenta luego del % reenganche, la misma dura desde tr a tc2, donde hago una apertura % definitiva. (*) % tf : Duracion total de la simulacion % (*) Opcionales, no se puede especificar tc2 sin haber especificado tr. % Ejemplo con tdata=(tc1,tf), falta y apertura. tdata=[0.4 2]; % Ejemplo con tdata=(tc1,tr,tc2,tf), falta,apertura,recierre,falta,apertura. % tdata=[0.3 0.4 0.4+0.6 10]; [delta,t,Barrasgen]=estabil(archivo,tdata,barracc,barrab1,barrab2); plot(t,delta) ylabel ('Angulo maquinas respecto a la slack en grados') xlabel ('t en segundos') legend(Barrasgen) grid El estudio de arriba esta hecho para el archivo pag516.m, caso falta trifásica en la línea Carga_6-Carga_5, sobre la barra Carga_6, se analiza la situación 1 con apertura de línea a los 0.4s. Definiendo el tiempo total de simulación en 2s.

Lo que da como resultado: Donde se ve que la red mantiene la estabilidad. Podemos aumentar el tiempo de duración de la falta cambiando dentro de destabil.m los parámetros de tdata, por ejemplo para que la falta dure 0.5s: tdata=[0.5 2]; Resulta evidente la perdida de estabilidad para esta situación.

tdata=[0.4 2] La misma situación del primer caso , donde el resultado era estable, pero aumentando en un 50% carga y generación se pierde la estabilidad: Obs.: En versiones anteriores a Matlab 6.1 puede haber problemas en la presentación de recuadro con los nombres de los generadores en el gráfico.

Asimismo para simular un evento completo: 3 – falta, apertura, reenganche (falta permanece) y apertura definitiva Red con la falta Red sin la línea que limpia la falta Red con la falta 0 0.3 0.4 1 15 Entonces tdata=[0.3 0.4 1 15]; Por otro lado si dejamos más tiempo antes de la apertura definitiva: tdata=[0.3 0.4 1.4 15]; Red sin la línea que limpia la falta El sistema deja de ser estable.

Descripción de las funciones estabil.m Función principal, para estudio de estabilidad flunrdr.m Función clásica para flujo de carga utilizando Newton-Raohson desacoplado rápido. red2mat.m * Esta función convierte un archivo ASCII con los datos de la red, en matrices utilizables por el Matlab, verifica conectividad, barras aisladas. Además de crear las barras internas de los generadores yest.m Calcula las tres matrices admitancias reducidas. es.m Función de entrada a la función del Matlab ode23.m para resolución numérica de ecuaciones diferenciales ordinarias. * ex fcm2dat.m funciones desarrolladas previamente funciones “nuevas”

Listado de las funciones function[delta,t,Barrasgen]=estabil(archivo,tdata,barracc,barrab1,barrab2) % Funcion para calculo de estabilidad transitoria (oscilaciones electromecanicas) % de sistemas electricos de Potencia. % % [delta,t,Barrasgen]=estabil(archivo,tdata,barracc,barrab1,barrab2) % % Argumentos de entrada % archivo : Datos de la red. % tadata=(tc1,tr,tc2,tf) vector tiempos de los sucesos, arranca en 0s. % tc1 : Clearing time 1, la falta se mantiene de 0 a tc1. % tr : Reenganche (opcional) desde tc1 la red esta sin la linea en falta % a partir de tr reengancho la linea (*) % tc2 : Si se especifica tc2 se asume que la falta se sustenta luego del % reenganche, la misma dura desde tr a tc2, donde hago una apertura % definitiva. (*) % tf : Duracion total de la simulacion % (*) Opcionales, no se puede especificar tc2 sin haber especificado tr. % barracc : Barra en falta trifasica tierra % barrab1 : Extremo 1 de la linea que limpia la falta. % barrab2 : Idem extremo 2. % % Argumentos de salida % delta : Vector diferencia de angulo en bornes de la maquina con respecto al % angulo de la maquina slack. % t : Vector tiempo. % Barrasgen : Nombre de las maquinas correspondientes a los respectivos deltas. % % R. Hirsch Junio 2002 global Sb N pN mv Barras f Y th ng H mEg Pm [N,pN,Barras]=red2mat(archivo); % Traigo los datos del archivo de la red. [mv,an,Pd,Qd,Pg,Qg,Qsh,maxerror,iter,Y]=flunrdr(N,pN); % Ejecuto flujo de carga. [Yrred,Yrcc,Yrcl,Yest]=Yest(Y,barracc,barrab1,barrab2); % Calculo matrices admitancias. % Yrred : Matriz admitancia reducida de la red original. % Yrcc : Matriz admitancia reducida con la barracc en falta trifasica (puesta a tierra). % Yrcl : Matriz admitancia reducida sin la linea que limpia la falta. estabil.m

S=(Pg-Pd)+j*(Qg+Qsh-Qd); % Potencia aparente en las barras. an=an*pi/180; [vx vy]=pol2cart(an,mv); % Paso a polar las tensiones. V0=vx+j*vy; I=(conj(S)/Sb)./conj(V0); % Corriente en las barras. % Calculo la tension y la potencia mecanica en bornes del generador % esto es, detras de la impedancia interna de la maquina. ng=0; % Inicializo contador de generadores. ngss=0; % Inicializo contador de generadores sin slack for i=pN(6,1):pN(6,2), ng=ng+1; Eg(ng)=V0(N(i,1))+I(N(i,1))*(N(i,3)+j*N(i,4)); % Tension en bornes de la maquina Pm(ng)=real(Eg(ng)*I(N(i,1))'); % Potencia de la maquina (asumida cte. durante toda la simulacion) H(ng)=N(i,5); if N(i,1)~=N(pN(3,1)), ngss=ngss+1; Barrasgen(ngss)=Barras(N(i,1)); % Nombre de las barras (sin slack) que corresponderan % correlativamente a los delta lambda que se calcularan. else ngslack=ng; end end t0=0; % Tiempo inicial tf=tdata(end); % Tiempo final de simulacion. tc1=tdata(1); % Clearing time 1 (primera apertura) t2=[]; % Inicializo como vectores vacios todoa los t3=[]; % potenciales resultados de tiempo y delta lambda t4=[]; % esto es, los sucesivos resultados de la ecuacion t5=[]; % swing. xf2=[]; xf3=[]; xf4=[]; xf5=[]; if length(tdata)>2, tr=tdata(2); % Si hay mas de dos parametros es porque se definio reenganche else tr=0; end if length(tdata)==4, tc2=tdata(3); % Si hay cuatro parametros es porque se definio segunda apertura definitiva % si no se define este parametro se asume que al hacer el reenganche se extinguio % la falta, caso contrario la falta subsiste luego del reenganche por lo tanto tengo % apertura trifasica definitiva. else tc2=0; end

% Primer tramo de simulacion, comun a todas las situaciones: tiempo de duracion % de la falta. clear t x x0=angle(Eg); mEg=abs(Eg); % Modulo tension en bornes (asumido cte. durante toda la simulacion) x0=[x0 0*x0]; % Vector condiciones iniciales, [angulos(sale del flujo de carga) velocidad angular (0)]. tspan=[0,tc1]; Y=abs(Yrcc); % Matriz reducida en condicion de falta. th=angle(Yrcc); [t1,xf1]=ode23('es',tspan,x0); % El sugundo tramo, red sin la linea en falta, dura hasta el final de la simulacion % o hasta el tiempo de reenganche si fue este especificado. if tr==0; ts=tf; fin_sim=1; else ts=tr; end tspan=[tc1,ts]; Y=abs(Yrcl); % Matriz reducida red sin linea en falta. th=angle(Yrcl); x0c=xf1(end,:); [t2,xf2]=ode23('es',tspan,x0c); if tr~=0 & tc2==0; % Ocurrio un reenganche y la falta se extinguio, o sea, % no especifique tc2. tspan=[tr,tf]; Y=abs(Yrred); % Matriz red reestablecida. th=angle(Yrred); x0c=xf2(end,:); [t3,xf3]=ode23('es',tspan,x0c); end % Ocurrio reenganche y se especifico tc2, esto es, la falta se mantiene % entre tr y tc2 y luego entre tc2 y tf apertura definitiva. if tc2~=0; % Implica necesariamente que tr es diferente de cero. tspan=[tr,tc2]; Y=abs(Yrcc); % Matriz en condicion de cortociruito. th=angle(Yrcc); x0c=xf2(end,:); [t4,xf4]=ode23('es',tspan,x0c); % Ultimo tramo apertura trifasica definitiva tspan=[tc2,tf]; Y=abs(Yrcl); % Matriz sin la linea en falta th=angle(Yrcl); x0c=xf4(end,:); [t5,xf5]=ode23('es',tspan,x0c); end t=[t1;t2;t3;t4;t5]; % Concateno todos los resultados de los x=[xf1;xf2;xf3;xf4;xf5]; % sucesivos tramos

% Calculo la diferencia entre el angulo de los generadores y el angulo del % generador slack. ii=0; for i=1:ng, if i~=ngslack, ii=ii+1; delta(:,ii)=180/pi*(x(:,i)-x(:,ngslack)); else, end end

function[Yrred,Yrcc,Yrcl,Yest]=Yest(Y,barracc,barrab1,barrab2) % Funcion para el calculo de las matrices admitancias reducidas para diferentes % estados de la red. % % [Yrred,Yrcc,Yrcl,Yest]=Yest(Y,barracc,barrab1,barrab2,N,pN,mv,Barras) % % Argumentos de entreda: % Ver fcm2dat, flunrdr y estabil % % Argumentos de salida: % Yrred : Matriz admitancia reducida de la red original. % Yrcc : Matriz admitancia reducida con la barracc en falta trifasica (barra puesta a tierra). % Yrcl : Matriz admitancia reducida sin la linea que limpia la falta. % Yest : Matriz red original sin reducir. % R. Hirsch 10 de Junio de 2002. global Sb N pN mv Barras Yest=Y; % Se carga la matriz admitancia clasica para flujos de carga. % Se le suma los valores de la impedancia de las máquinas (Ra +jXd'). for i=pN(6,1):pN(6,2), b1=N(i,1); % Nombre barra 1 b2=N(i,2); % Nombre barra 2 y=1/(N(i,3)+j*N(i,4)); Yest(b1,b2)=-y; Yest(b2,b1)=-y; Yest(b1,b1)=Yest(b1,b1)+y; Yest(b2,b2)=Yest(b2,b2)+y; end % y también las cargas (asumiendo impedancia constante). for i=pN(1,1):pN(3,1), b1=N(i,1); y=(N(i,4)-j*N(i,5)+j*N(i,10))/(Sb*mv(i)^2); Yest(b1,b1)=Yest(b1,b1)+y; end nBtotal=max(N(:,2)); nB=pN(3,1); yest.m

% Calculo matriz reducida red original Yrred=Yest(nB+1:end,nB+1:end)-Yest(nB+1:end,1:nB)*inv(Yest(1:nB,1:nB))*Yest(1:nB,nB+1:end); % Falta trifasica a tierra implica que barracc es parte del sistema de referencia de tierra, esto % es, no aporta ecuacion por lo que a priori al calculo de la matriz reducida en condicion de % falta se elimina la fila y columna correspondiente a barracc. nbcc=find(strcmpi(Barras,barracc)); Yestcc=Yest; Yestcc(nbcc,:)=[]; Yestcc(:,nbcc)=[]; Yrcc=Yestcc(nB:end,nB:end)-Yestcc(nB:end,1:nB-1)*inv(Yestcc(1:nB-1,1:nB-1))*Yestcc(1:nB-1,nB:end); % Para calcular la matriz reducida sin la limpia que limpia la falta tengo que % primero eliminar de la matriz original los aportes de admitancia y suceptancia % correspondientes a esta linea. Yestcl=Yest; nbl1=find(strcmpi(Barras,barrab1)); % Se busca el número de barra que corresponde b1 nbl2=find(strcmpi(Barras,barrab2)); % Se busca el número que le corresponde a la barra b2 fb=(N(:,1:2)==nbl1); % fb es un matriz de dos columnas x filas de N, donde hay unos % cuando coincide el nombre de la barra 1 fb=fb+2*(N(:,1:2)==nbl2); % A la matriz anterior le sumo dos cuando hay coincidencia % con el nombre de la barra 2. fN=find(sum(fb')==3); % Las filas de fN que sumen tres son las buscadas. a=1; % En principio estoy considerando solo lineas bl=N(fN,5)/2; % calculo la suceptancia y=1/(N(fN,3)+j*N(fN,4)); % Admitancia serie. Yestcl(nbl1,nbl1)=Yestcl(nbl1,nbl1)-j*bl-y; % Subtraigo elementos de la diagonal Yestcl(nbl2,nbl2)=Yestcl(nbl2,nbl2)-j*bl-y/(a*a); Yestcl(nbl1,nbl2)=Yestcl(nbl1,nbl2)+y/a; % Subtraigo (sumo) Elementos fuera de la diagonal. Yestcl(nbl2,nbl1)=Yestcl(nbl1,nbl2); Yrcl=Yestcl(nB+1:end,nB+1:end)-Yestcl(nB+1:end,1:nB)*inv(Yestcl(1:nB,1:nB))*Yestcl(1:nB,nB+1:end);

function xpri = es(t,x) % Representacion de la ecuacion swing de un sistema % multimaquina en la forma espacio-estado. % Potencia mecanica y modulo de la tension son asumidos constantes. global Sb N pN mv Barras f Y th ng H mEg Pm % Calculo la potencia electrica en bornes de la maquina para un % dado estado la red representada por la matriz admitancia. Pe=zeros(1, ng); for ii = 1:ng for jj = 1:ng Pe(ii) = Pe(ii) + mEg(ii)*mEg(jj)*Y(ii, jj)*cos(th(ii, jj)-x(ii)+x(jj)); end end % Ecuacion swing for k=1:ng xpri(k)=x(k+ng); xpri(k+ng)=(pi*f)/H(k)*(Pm(k)-Pe(k)); end xpri=xpri'; es.m