Despacho Optimo de la Generación Flujo de Carga: Estimamos valores razonables de Pgen de las barras PV adicionalmente Pgen de la barra Slack es calculada por: Despacho Optimo: Pgen de las barras PV e incluso de la slack se calculan tal que el costo total de la generación sea mínimo. Función objetivo (Gp:) Flujo de carga Despacho Min. costo Costos de la generación. Límites min. y max. Datos de la red no FIN si Despacho Optimo o (más general) Flujo de Carga Optimo
Optimización de una función sujeta a restricciones de igualdad El problema es minimizar la función costo: Sujeta a restricciones de igualdad Tales problemas pueden resolverse por el método de los multiplicadores de Lagrange. Se crea una función aumentada introduciendo un vector de k elementos ?: Los valores de que minimizan f sujeto a la igualdad g son los que resuelven las siguientes ecuaciones: Ejemplo 7.1: Hallar el mínimo de la función: (cuadrado de la distancia del origen hasta x,y). Sujeto a la restricción:
Formamos la función de Lagrange: Las ecuaciones a resolver son: En muchos problemas la solución directa no es posible por lo que las ecuaciones arriba son resueltas iterativamente. De las dos primeras ecuaciones, encontramos x e y: Sustituyendo en la tercera ecuación resulta en: La que puede ser resuelta por Newton-Raphson: CORREGIR
Empezando con un valor estimado de ?, un nuevo valor es encontrado. El proceso se repite en la dirección del gradiente decreciente hasta que ?f(?) es menor que un ? especificado. Este método es conocido como el método del gradiente. Para la función arriba el gradiente es: Utilizar la función ‘te6ej1’ para resolver la ecuación de f(?), luego calcular x e y. Hallar el mínimo o el máximo dependerá de la dirección del gradiente, ¿Para que rango de estimación inicial de ? hallaremos un mínimo y para cual un máximo? Optimización de una función sujeta a restricciones de igualdad y restricciones de desigualdad El problema es ahora minimizar la función costo: Sujeta a restricciones de igualdad Se trata de formular una extensión de los multiplicadores de Langrange a los efectos de incluirlas restricciones, este método generalizado se le conoce como condiciones necesarias de optimalidad de Kuhn-Tucker. En la expresión abajo se incluye entonces un vector ?j de m elementos indeterminados a los efectos de considerar las m restricciones de desigualdad: Y a restricciones de desigualdad
Siendo las condiciones necesarias las siguientes: Ejemplo: Hallar el mínimo de la función: (cuadrado de la distancia del origen hasta x,y) Sujeto a la restricción: Y a la desigualdad: Si el problema no está planteado de la misma forma los signos de los multiplicadores podrías ser diferentes: Planteando Las condiciones de Kuhn-Tucker son:
Si Las ecuaciones a resolver son: Tenemos que: Resolviendo por Newton-Raphson: Sabemos que de la resolución de las tres primera ecuaciones que x=4 e y=3 Lo que viola la condición de desigualdad de la cuarta ecuación, por lo tanto de la quinta ecuación se debe cumplir que:
COSTO OPERATIVO DE LAS CENTRALES TERMICAS En todos los casos prácticos el costo del generador i puede ser representado como: Una característica importante es la derivada del costo respecto a la potencia activa, lo que se conoce como costo incremental: $/h Pi MW ?i $/MWh Pi MW
Aplicando el método de los multiplicadores de Lagrange: Y planteando las respectivas ecuaciones : Despacho óptimo de las unidades de generación sin considerar pérdidas ni límites de generación. C1 P1 C2 P2 Cng Png PD Nuestra función objetivo es entonces: Sujeta a la restricción: La primera condición resulta en : Pero como : La condición para el despacho óptimo: todos los generadores tengan el mismo costo incremental
La segunda condición: Método analítico de resolución: Por un lado tenemos Para cada generador (i=1,…,ng) se las conoce como ecuaciones de coordinación. Tenemos que determinar el valor de , de la segunda condición: De donde: Ejemplo: El costo total de tres plantas térmicas en $/h está dada por: Donde P1, P2 y P3 están en MW. La demanda total PD es 800MW. Sin considerar pérdidas ni límites en la generación, encontrar el despacho óptimo y el costo total en $/h.
Sustituyendo ? en las ecuaciones de coordinación: El costo total es entonces: Interpretación gráfica: P, MW $/MWh 8.5 150 250 400
Ejemplo incluyendo límites en la generación: El costo total de tres plantas térmicas en $/h está dada por: Donde P1, P2 y P3 están en MW. La demanda total PD es 975 MW. Los límites de generación son: Sin considerar pérdidas, encontrar el despacho óptimo. P1 viola el ímite de Pmax, por lo que la “pego” al tope de 450MW y redespacho las otras dos
P, MW $/MWh 9.16 187 305 483 450 8.9 325 9.4 200 El costo total sin considerar las restricción de P3: Con la restricción:
Despacho Económico Optimo Incluyendo Restricciones en la Generación y Pérdidas Como ya hemos visto, en todoos los casos prácticos el costo del generador i puede ser representado como: Por lo tanto, la función aminimizar(función objetivo) es: Sujeta a la restricción de igualdad: Y a las desigualdades: Usando los multiplicadores de Lagrange y los terminos adicionales para incluir las desigualdades: Queda entendido que: o sea, si las restricciones de desiguladad no son violadas los correpondientes terminos no existen. Una práctica común para incluir el efecto de las pérdidas de la transmisión es expresar las pérdidas totales de la transmisión como una función cuadrática de las potencias de las unidades generadoras, cuya forma más general es: Se la conoce como la fórmula de Kron, y los coeficientes B son llamados coeficientes de pérdidas o coeficientes-B, más adelante se presenta la obtención de los mismos.
Los valores de que minimizan L son los que anulan las derivadas parciales: La primera condición, y resolviendo el problema sin considerar en primera instancia las restrcciones de desigualdad: resulta en: como: la condición resulta en: Incremental de perdidas de transmisión Incremental del costo de generación “Se activan” cuando alguna o algunas restricciones son violadas en uno o varios generadores: Es común reordenarla como: Factor de penalidad del generador i o El incremental de las pérdidas de transmisión vale: Además sabemos que: Sustituyendo respectivamente en la expresión arriba
Reordenando los término de la siguinete forma: Extendiendo la ecuación arriba a todas las plantas resulta en el siguiente sistema linear de ecuaciones representado en su forma matricial: O en su forma abreviada: En la práctica se resuelve: P=E D De la segunda condición: Siendo: Sustituyedo, nos queda: o: La resolvemos por Newton-Raphson, siendo entonces ?(0) la estimación inicial y ??(0) la pequeña desviación de la solución correcta tenemos:
Expandiendo en series de Taylor hasta el término de primer orden: o: finaemente: Y se repite el proceso hasta que: Es menor que un dado valor de precisión especificado. A partir de la segunda iteración los valores de P de las distintas unidades generadores se obtienen del sistema de ecuaciones lineares que resuelve la primera condición: P=E D Una vez que converge se verifica si alguna máquina viola alguno de sus límites de generación, si esto es así, la o las unidades correspondientes pasan a generar un valor igual al límite que correspondiente, y se vuelve a entrar en el algoritmo de Newton-Raphson, siendo entonces los valores de generación de estas máquinas parámetros dados y ya no incognitas.
Funciones matalb desarrolladas: despacho.m – función principal donde se implementa el algoritmo presentado. op2dat.m – función del estilo de red2.mat, desde donde se lee un archivo ascii con los datos de ls costos de las máquinas y sus límites operativos y se guardan en variables a ser usadas por las demás funciones. costogen.m – cálculo del costo total de la generación. costoB.m – se calcula los coeficientes B d la fórmula de Kron para el cálculo de las pérdidas en un sustema de transmisión.. daledes.m – rutina para corrida “facil” de la aplicación y realiza el procesos iterarativo flujo de carga despacho óptimo. function[]=daledes(archivo,archivo2) [N,pN,Barras]=red2mat(archivo); [mv,an,Pd,Qd,Pg,Qg,Qsh,maxerror,iter,Y]=flunrdr(N,pN); global Sb Zbus=full(inv(Y)); [B,B0,B00,PL]=coefB(pN,mv,an,Pd,Qd,Pg,Qg,Qsh,Zbus); [costo,mwlimites]=opt2dat(archivo2,N,pN,Barras); lambda=7; Pgg=Pg(pN(2,1):pN(3,1)); [costototal]=costogen(Pgg,costo) [Nopt,dpslack,lambda,Pgg,PL]=despacho(Pd,Pg,costo,B,B0,B00,pN,N,mwlimites,lambda); while dpslack>0.001, [mv,an,Pd,Qd,Pg,Qg,Qsh,maxerror,iter,Y]=flunrdr(Nopt,pN); [B,B0,B00,PL]=coefB(pN,mv,an,Pd,Qd,Pg,Qg,Qsh,Zbus); [Nopt,dpslack,lambda,Pgg,PL]=despacho(Pd,Pg,costo,B,B0,B00,pN,N,mwlimites,lambda); end tabbar [costototal]=costogen(Pgg,costo) save ejemplo5b.dat Nopt lambda PL Pgg -append Si interesa salvar las variables hay que cambiar a mano el nombre del archivo
Ejemplo Dado la red abajo, con los valores estimados de despacho de potencia reactiva, determinar el despacho óptimo. V1=1.06?0° |V2|=1.045 |V3|=1.03 30 MW 40 MW 1 2 3 4 5 20 MW 10 MVar 20 MW 15 MVar 50 MW 30 MVar 60 MW 40 MVar 0.08+j0.24 0.02+j0.06 0.04+j0.12 0.06+j0.18 0.06+j0.18 0.08+j0.24 0.01+j0.03 % DATOS DE BARRA % CARGA GENERACION min max Shunt % BARRA TENSION MW MVAR MW MVAR MVAR MVAR MVAr SL 1 1.06 0 0 0 0 0 0 0 PV 2 1.045 20 10 40 30 10 50 0 PV 3 1.03 20 15 30 10 10 40 0 PQ 4 1.00 50 30 0 0 0 0 0 PQ 5 1.00 60 40 0 0 0 0 0 % % DATOS DE LINEAS % BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIA Linea 1 2 0.02 0.06 0.060 Linea 1 3 0.08 0.24 0.050 Linea 2 3 0.06 0.18 0.040 Linea 2 4 0.06 0.18 0.040 Linea 2 5 0.04 0.12 0.030 Linea 3 4 0.01 0.03 0.020 Linea 4 5 0.08 0.24 0.050 % DATOS PARA DESPACHO OPTIMO DE LA GENERACION % archivo: ejemplo5b % % BARRA C1 C2 C3 Pmin Pmax Slack 200 7.0 0.008 10 85 Gen_1 180 6.3 0.009 10 80 Gen_2 140 6.8 0.007 10 70 Coeficientes de la función costo de la generación y límites operativos de los generadores
» daledes('ejemplo5b.m','ejemplo5b_opt.m') Flujo de carga no optimo Máximo error en la potencia = 0.058002 No. de Iteraciones = 4 Barra Tensión Angulo ——Carga—— —Generación— Shunt Mag. grados MW MVAr MW MVAr MVAr Carga_1 1.019 -3.248 50.0 30.0 0.0 0.0 0.0 Carga_2 0.990 -4.406 60.0 40.0 0.0 0.0 0.0 Gen_1 1.045 -1.782 20.0 10.0 40.0 39.4 0.0 Gen_2 1.030 -2.664 20.0 15.0 30.0 23.3 0.0 Slack 1.060 0.000 0.0 0.0 83.0 7.3 0.0 Total 150.0 95.0 153.0 70.0 0.0 costototal = 1.6332e+003 Flujo de carga optimo Máximo error en la potencia = 0.0593484 No. de Iteraciones = 4 Barra Tensión Angulo ——Carga—— —Generación— Shunt Mag. grados MW MVAr MW MVAr MVAr Carga_1 1.019 -1.199 50.0 30.0 0.0 0.0 0.0 Carga_2 0.990 -2.717 60.0 40.0 0.0 0.0 0.0 Gen_1 1.045 -0.270 20.0 10.0 69.8 28.3 0.0 Gen_2 1.030 -0.481 20.0 15.0 59.1 13.2 0.0 Slack 1.060 0.000 0.0 0.0 23.2 25.9 0.0 Total 150.0 95.0 152.1 67.3 0.0 costototal = 1.5973e+003
Método de Kron para obtenerlas perdidas del sistema en función de la potencia activa del parque generador: Matriz impedancia
PQ PV+Sl PV+Sl 1
C=C1 * C2 PV+Sl+1 PV+Sl+1 PQ PV+Sl
Los coeficientes B son en valores pu, cuando la potencia está expresada en MW, los coeficientes B valen: Bij= Bij pu/Sb, Boi= Boi , B00= B00 pu*Sb Donde Sb son los MVA Base
Función coefB : Esta función calcula los coeficientes de perdidads B, dada una red con su respectivo flujo de carga. Argumentos de entrada: • Matriz pN puntera de la matriz N. • Resultado del flujo de carga:. • Matriz Zbus, inversa de Ybus Argumentos de salída: • Coeficientes de perdidas B. • Perdidas totales en MW. function[B,B0,B00,PL]=coefB(pN,mv,an,Pd,Qd,Pg,Qg,Qsh,Zbus) global Sb; V=mv.*exp(j*deg2rad(an)); % Tensión compleja. Il=-1/Sb*(Pd-j*Qd)./conj(V); % Corrientes de carga en todas las barras. ID= sum(Il); % Sumatoria de las corrientes (ec. [13]). l=Il/ID; % Fracción de la corriente total (vector ec. [15]). sl=pN(3,1); % Ubicación de la barra slack. T=Zbus(sl,:)*l; % Cálculo de T (ec.[18]). nB=sl; % Número total de Barras. fPQ=pN(1,2); % Número de barras PQ. iPV=pN(2,1); % Ubicación inicio de las barras PV. ng=nB-fPQ; % Número total de barras de generación (Slack+PV). W(1:ng) = Zbus(sl,iPV:sl)/T; % Cálculo de w (arriba ec. [23], definiendo …). C1gg=eye(ng,ng); % Sub matrices que se concatenan para armar C1. C1g=[zeros(fPQ, ng);C1gg]; C1=[C1g,l]; % C1 C2gD=[C1gg;-W]; % Sub matrices que se concatenan para armar C2. C2D=zeros(ng,1); CnD=[C2D;-W(ng)]; C2=[C2gD,CnD]; % C2 C=C1*C2; % C al=(1-j*((Qg(iPV:sl)+Qsh(iPV:sl))./Pg(iPV:sl)))./conj(V(iPV:sl)); % Elementos al=al.'; % para armar la matriz alpha (ec. [28]). alp=[al, -V(sl)/Zbus(sl,sl)]; % Útimo elemento de la diagonal de la matriz alpha. alpha=diag(alp); % Obtención de la matriz alpha (segun ec. [30]). H = real(alpha*conj(C)'*real(Zbus)*conj(C)*conj(alpha)); % Cálculo de H (ec. [34]). B=H(1:ng,1:ng); % Partición de la matriz H conforme ecuación [36]. B0=2*H(ng+1,1:ng); B00=H(ng+1,ng+1); PL = Pg(iPV:sl)'*(B/Sb)*Pg(iPV:sl)+B0*Pg(iPV:sl)+B00*Sb; % Perdidas totales (ec.[36]) % convirtiendo los valores pu de los coeficiente B. (Gp:) 1
Comentarios: (Gp:) 1 (Gp:) – C1g C1gg C2gD C2D CnD
Ejemplo Dado la red abajo, calcular los coficientes B y las perdidas totales de la red. V1=1.06?0° |V2|=1.045 |V3|=1.03 30 MW 40 MW 1 2 3 4 5 20 MW 10 MVar 20 MW 15 MVar 50 MW 30 MVar 60 MW 40 MVar 0.08+j0.24 0.02+j0.06 0.04+j0.12 0.06+j0.18 0.06+j0.18 0.08+j0.24 0.01+j0.03 % DATOS DE BARRA % CARGA GENERACION min max Shunt % BARRA TENSION MW MVAR MW MVAR MVAR MVAR MVAr SL 1 1.06 0 0 0 0 0 0 0 PV 2 1.045 20 10 40 30 10 50 0 PV 3 1.03 20 15 30 10 10 40 0 PQ 4 1.00 50 30 0 0 0 0 0 PQ 5 1.00 60 40 0 0 0 0 0 % % DATOS DE LINEAS % BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIA Linea 1 2 0.02 0.06 0.060 Linea 1 3 0.08 0.24 0.050 Linea 2 3 0.06 0.18 0.040 Linea 2 4 0.06 0.18 0.040 Linea 2 5 0.04 0.12 0.030 Linea 3 4 0.01 0.03 0.020 Linea 4 5 0.08 0.24 0.050
clear [N,pN]=red2mat('ejemplo5b.m'); [mv,an,Pd,Qd,Pg,Qg,Qsh,maxerror,iter,Y]=flunrdr(N,pN); Zbus=full(inv(Y)); [B,B0,B00,PL]=coefB(pN,mv,an,Pd,Qd,Pg,Qg,Qsh,Zbus) B = 0.0228 0.0017 0.0093 0.0017 0.0179 0.0028 0.0093 0.0028 0.0218 B0 = 0.0031 0.0015 0.0003 B00 = 3.0523e-004 PL = 3.0525 Las perdidas totales de la red son de 3.0525 MW