FORMULACIÓN DEL PROBLEMA. LA FUNCIÓN LP...

FORMULACIÓN DEL PROBLEMA. LA FUNCIÓN LP SOLVE. Universidad de Sevilla.

Jorge Vázquez Ballesteros. Ingeniería Industrial.

43

CAPITULO 3: “FORMULACIÓN DEL PROBLEMA. LA

FUNCIÓN LP SOLVE”.

En este capítulo se hablará de los distintos nodos de la red y que parámetros hay que fijar en cada

uno de ellos. Se explicará detalladamente la función de MATLAB, lp solve, y las restricciones

técnicas incluidas que debe respetar. Por último se hablará de un parámetro que servirá como

condición de parada del algoritmo, al que se le ha llamado k_min, y se incluirá un diagrama de flujo

donde se explicará los pasos que realiza el algoritmo de optimización.

3.1 TIPOS DE NODOS DE LA RED.

En la red modelo podemos identificar tres tipos de nodos:

1. Nodos de cargas (o consumo).

Se encuentra las cargas y se fija la potencia activa y la potencia reactiva requeridas por cada carga.

2. Nodos generadores.

Son nodos donde se encuentran los generadores distribuidos en la red, controlables y fijos.

En estos nodos se fija la potencia activa y el valor eficaz de la tensión.

¿Por qué se fija la tensión y no la potencia reactiva?

- Fijando la tensión (generalmente entre el valor nominal y un ±5% del valor nominal), se fija

el nivel de tensión en otros nodos de referencia, a menudo dispersos por toda la red, por lo

tanto la solución de las ecuaciones proporciona ya desde el primer intento, una solución

aceptable. Además facilita la convergencia del procedimiento iterativo para la resolución de

las ecuaciones de red.

- Se conoce también que la potencia reactiva Q, de los generadores se puede variar entre un

límite máximo en sobrexcitación y un límite mínimo en subexcitación, variando la corriente

de excitación.



44

Por esto podemos dejar libre la potencia reactiva, y esta variará de forma que permita alcanzar la

tensión predeterminada.

3. Nodo de saldo (Slack Bus).

Asociaremos este nodo a un generador que suministre una potencia considerablemente alta.

La definición de este nodo es necesaria ya que hasta ahora hemos fijado la potencia activa en los

nodos de carga y en los nodos generadores. Esto nos lleva a pensar inmediatamente que el problema

estaría ya vinculado de manera que sabríamos las pérdidas en la red, y no podríamos intervenir para

modificar nada.

Las pérdidas en la red, lógicamente, no se conocen inicialmente. Éstas se calculan a la vez que el

flujo de cargas, una vez resueltas las ecuaciones de la red.

Es por esto que en el nodo de saldo se fija el valor eficaz y la fase de la tensión. Su potencia activa,

en cambio, se obtendrá como balance de las potencias activas de las cargas, de los generadores y de

las pérdidas en la red.

Fijaremos la fase de la tensión en el nodo de saldo, n=0, de esta forma las fases de todos los nodos

de la red se medirán respecto a ésta.

3.2 FUNCION LP_SOLVE.

(Solves mixed integer linear programming problems).

3.2.1 Finalidad.

En el capítulo 2, se muestra el proceso de linealización del problema. Ahora con las ecuaciones ya

linealizadas escogemos una función para minimizar nuestro problema. La función elegida para

resolver el problema de mínimo con restricciones se llama lp_solve. Es una función de Matlab, que

es el programa que se ha elegido para programar por las características que posee explicadas en el

capítulo precedente. Esta función posee dos características principales: es de tipo lineal y trabaja

con variables enteras y/o continuas.

Synopsis: [obj,x,duals,stat] = lp_solve(f,a,b,e,vlb,vub,xint,scalemode,keep)



45

Esta función de Matlab resuelve el siguiente problema de tipo MILP (mixed integer linear

program):

max v = f'*x

a*x <> b

vlb <= x <= vu

x (int) = entero

Argumentos: Los cuatro primeros términos son obligatorios,

Variables

de

ingreso

Descripción

Dimensiones

f Vector de coeficientes de la función objetivo lineal.

n

a Matriz que representa las restricciones lineales.

[m,n]

b Vector de elementos conocidos del vínculo de desigualdad.

m

e

Vector que determina el tipo de vínculo (>,<,=)

e(i) = -1 si es <=

e(i) = 0 si es =

e(i) = 1 si es >=

m

vlb

Vector de cotas inferiores de las variables a optimizar.

(Si esta vacío o no está, el programa lo entenderá como un vector de

ceros).

m

vub Vector de cotas superiores de las variables a optimizar.

(Puede estar vacío o no estar).

m

xint Vector que indica cuales variables de x son enteras.

scalemode Bandera de escalera.

(Estará en “Off” cuando es cero o no está definida).

keep Bandera que mantiene el problema lp después de haber sido resuelto.

(Si se omite el lp se eliminará una vez resuelto).

Tabla3.2.1a Variables de ingreso de la función lp solve.



46

“Outputs”: (Si no se encuentra una solución se devolverá un vector vacío).

Variables

de salida Descripción

Dimensiones

x Vector de variables de optimización. m

obj Valor final de la función objetivo. escalar

duals Solución del problema dual.

Tabla 3.2.1b Variables de salida de la función lp solve.

Se ha realizado un programa para ingresar las variables de entrada a la función lp_solve, y la

posterior llamada a éste, ha sido el: lpsolve22, a continuación se explicará, las consideraciones

efectuadas para implementar este programa en Matlab (el programa completo está en el Apendice

B).

[obj,x]=lp_solve(fobj,Alin,blin,e,LB,UB,xint)

Y queremos que resuelva:

min fobj'*x

a*x <= b

LB <= x <= UB

xint = entero

3.2.2 Variables de entrada.

fobj: Es el vector que, multiplicado por las variables a optimizar del problema, constituye la

función a maximizar (que se multiplicará por -1 para obtener un mínimo en lugar de un máximo).



47

Será un vector columna con tantas filas como variables tenga el problema.

El vector fobj lo definiremos como:

[

]

donde

g es el número de generadores controlables.

T es el número de transformadores regulables.

nd es el númoer de nodos de la red.

nl es el número de ramas de la red.

Alin: Es la matriz de coeficientes que, multiplicada por las variables a optimizar del problema,

constituye el término izquierdo de las ecuaciones de vínculos.

Esta matriz de matrices sería la siguiente:

[

]

3nl (3*número de lineas)

3nd (3*números de nodos)

1

g

2g+T



48

Matriz A C H

Dimensiones (2g+4+1) x (2g+T) (3nd+1+g) x (2g+T) 3nl x (2g+T)

Componentes Vínculos técnicos.

(Ver apartado 3.2.4)

Coeficientes

de sensibilidad de las

tensiones y de las

pérdidas.

Coeficientes de

sensibilidad de las

corrientes.

Tabla 3.2.2a Descomposición de la matriz de entrada Alin.

[

]

[

]

[

]

[

]

[

]



49

[

]

[

]

Donde cada uno de sus términos ha sido calculado y explicado en el capítulo 2.

Matriz

Dimensiones (2g+4+1)

x

(3nd+1+g)

(2g+4+1)

x

(3nl)

(3nd+1+g)

x

(3nd+1+g)

(3nd+1+g)

x

(3nl)

(3nl)

x

(3nd+1+g)

(3nl)

x

(3nl)

Componentes Ceros. Ceros. Diagonal de

-1.

Ceros. Ceros. Diagonal de

-1.

Tabla 3.2.2b Descomposición de la matriz de entrada Alin.

blin: Es el vector que contiene los términos conocidos de las ecuaciones de vínculos.

[

]



50

[

]

[

]

[

]

[

]

[

]

Donde el vector b se refiere a los términos conocidos de los vínculos técnicos (referido a la matriz

A), el vector dtr es el vector de las diferencias entre los sucesivos valores de las variables respecto

de los valores iniciales, de las pérdidas activas de la red y de las diferencias entre los valores de las

potencias activas de los generadores respecto a los valores de set point (este vector hace referencia

a la matriz C) e i_tr es el vector de los valores absolutos de las corrientes en las líneas (se refiere a

la matriz H).



51

e: Vector cuyas componentes son -1 porque queremos que la desigualdad tome el signo <=. Este

vector tendrá la misma dimensión del vector x.

[

]

LB: Es el vector que contiene los valores inferiores que podrán tomar las variables del problema.

[

]

[

(

)

(

) ]

[ ]

[ ]

El vector lb se refiere a los valores mínimos que podrán tomar las P,Q y n. El lbs indica los valores

mínimos que podrán tomar los coeficientes de sensibilidad de la matriz C. Y el vector lbs1 se refiere

a aquellos valores mínimos en módulo, que podrán tomas las corrientes en cada rama.

UB: Es el vector que contiene los valores superiores que podrán tomar las variables del problema.

[

]

2g

T

3nd+1+g

3nl



52

[

(

)

(

) ]

[

]

[

]

El vector ub se refiere a los valores máximos que podrán tomar las P,Q y n. El ubs indica los

valores máximos que podrán tomar los coeficientes de sensibilidad de la matriz C. Y el vector ubs1

se refiere a aquellos valores máximos en módulo, que podrán tomas las corrientes en cada rama.

xint: Es un vector que contiene los índices de línea de los elementos que pertenecen al vector de

variables x, que queremos que tengan valores enteros.

3nd+1+g

3nl

T

2g



53

3.2.3 Variables de salida.

obj: Es el valor que tomará la función objetivo una vez que las variables hayan sido optimizadas.

En nuestro caso será:

∑ | |

∑ | |

son los pesos que se le dará a cada término de la función.

x: es el vector solución del problema.

con i=1, 2,…, g (número de generadores controlables).

con i=1, 2,…,g

con i=1, 2,…, T (número de transformadores regulables).

es la variación de posición del transformador.

Utilizamos el vector y no el , porque el programa nos da la posibilidad que algunas variables

asuman valores enteros, y esto es muy útil en lo referente a encontrar el set point óptimo para los

transformadores regulables.

Veamos cuales serían las relaciones que nos permiten utilizar el vector en lugar de . La

tensión en vacío del secundario se puede expresar como:

(3.2.3.1)

donde la relación entre y es la inversa de la relación de espiras.



54

Las variables referidas al circuito primario son constantes, pero las concernientes al secundario se

pueden variar, de manera que para la regulación podemos cambiar el número de espiras , de la

siguiente forma:

(3.2.3.2)

podrá tomar valores positivos o negativos pero siempre enteros.

El incremento de tensión en el secundario asociado a un incremento de la posición del

transformador, lo implementamos directamente en lpsolve22 mediante la fórmula:

(3.2.3.3)

donde es la tensión en vacío del secundario del transformador, (voltage per tap) es la

variación de tensión del secundario respecto a un escalón unitario y positivo aplicado a la posición

del transformador, expresado en tanto por ciento respecto al valor de tensión nominal.

Imponiendo vpt =1.5% y V2n= 20 kV y sabiendo que la posición del transformador puede variar

entre -8 y +8, la tensión del secundario del transformador podrá tomar los siguientes valores:

V2n [kV] V2n [kV]

-8 17.6 1 20.3

-7 17.9 2 20.6

-6 18.2 3 20.9

-5 18.5 4 21.2

-4 18.8 5 21.5

-3 19.1 6 21.8

-2 19.4 7 22.1



55

-1 19.7 8 22.4

0 20

Tabla 3.2.3 Valores posibles para la tensión en vacío del secundario de los transformadores.

3.2.4 Restricciones técnicas (vínculos de desigualdad).

Los vínculos de desigualdad impuestos al problema, se refieren a los límites entre los cuales pueden

variar los valores de set point asignados a las diferentes variabales controlables.

Estos vínculos se expresarán como:

1. Restricción sobre la potencia activa suministrada por los generadores.

(3.2.4.1)

Esta restricción permitirá a los generadores trabajar siempre, para cada variación impuesta por el

optimizador, entre una potencia activa mínima y una máxima.

k es la iteración actual, entonces:

, restando a cada miembro la potencia

obtenida en la iteracion precedente:

(3.2.4.2)

Donde todos los términos son conocidos menos el de la potencia activa de la iteración k.

Los términos y

(potencia activa mínima y máxima de generador genérico j-ésimo) son

leídos por el JavaScript, el cual está enlazado con los datos que hemos impuesto en la interfaz

gráfica del EMTP-RV.

2. Restricción sobre la potencia reactiva suministrada por los generadores.

(3.2.4.3)



56

Los valores mínimos y máximos de la potencia reactiva de los generdores no son datos, a priori, que

el JavaScript pueda leerlos, ya que estos valores no se escriben en el EMTP-RV, la forma de

calcularlos es:

(3.2.4.4)

(3.2.4.5)

es el factor de potencia. Y ademas:

(3.2.4.6)

La curva de Capability estaría simplificada a un rectángulo de valores posibles para P y Q:

Figura 3.2.4 Representación de los límites de potencia activa y reactiva en la curva de “Capability”.

Al variar la potencia activa de uno de los generadores, ∆Pj , la variación de potencia reactiva

consecuente será:

(3.2.4.7)

y la variación de potencia reactiva total no podrá ser superior a:

(3.2.4.8)

Por lo tanto la restricción :

, (3.2.4.9)

Que vendrá impuesta como:

(3.2.4.10)

Y separando los términos conocidos:

(3.2.4.11)

se tiene :

P

Q

Pmáx

Qmáx

φmín



57

(3.2.4.12)

3. Restricción de la tensión en vacío del secundario de los transformadores.

(3.2.4.13)

Para referirnos al paso iterativo k:

(3.2.4.14)

Y desarrollando tenemos:

(3.2.4.15)

En este caso y

se calculan imponiendo dos valores límites a la relación de espiras (-8 y

8).

La posición inicial de n, la relación de espiras, es n=0. Esto quiere decir que la tensión en vacío en

el secundario de los transfomadores será igual a la nominal, que en nuestro caso será V2n = 20 kV.

Los nuevos valores de tensiones al variar n, serán calculados de la siguiente forma:

(

) (3.2.4.16)

Donde vpt es la variación de tensión asociada a un incremento unitario de la relación de espiras,

expresada en tanto por ciento respecto a la nominal.

Entonces los límites de tensión, los calcularemos:

(

) (3.2.4.17)

(

) (3.2.4.18)

Los valores mínimo y máximo de n, serán escritos en la interfaz gráfica del EMTP-RV.



58

4. Restricción de la potencia máxima intercambiable con la red de transmisión.

(3.2.4.19)

Esta restricción se refiere al balance de potencia que se hace entre la red de transimsión y la de

distribución. La potencia del nodo de saldo, llamado slack, podrá ser absorbida o suministrada.

La potencia del nodo de saldo se calcula siempre como un balance de tres términos, la potencia de

los generadores controlables y de aquellos no controlables, la potencia requerida por las cargas de la

red y las pérdidas de potencia activa, en definitiva:

∑

∑

∑

(3.2.4.20)

donde

nc es el número de cargas.

g es el número de generadores controlables.

G es el número total de generadores (los controlables más los no controlables). Por tanto el

último término se refiere a la potencia activa suministrada por los generadores no

controlables.

El nodo de saldo se ha considerado como otro generador. Si la Pslack >0 entonces hablaremos de

potencia absorbida y si Pslack <0 entonces será suministrada.

Tanto las potencias de los generadores controlables como la potencia de pérdidas, son variables de

nuestro problema, por lo tanto variarán en cada paso iterativo k:

∑

∑

Separando términos conocidos:

∑

∑

∑

∑

(3.2.4.21)



59

∑

∑

∑

∑

(3.2.4.22)

En lo concerniente a la restricción sobre la potencia reactiva,

(3.2.4.23)

procedemos de la misma forma, obteniendo así:

∑

∑

∑

∑

(3.2.4.24)

∑

∑

∑

∑

(3.2.4.25)

La potencia reactiva del nodo de saldo se toma directamente del EMTP-RV, a diferencia de lo que

sucede con en el caso de los generadores dispersos, que había que implementar la ecuación (3.2.4.4),

en el Matlab para poder conocer el valor.

Tenemos finalmente todos las componentes de la matriz A, que como referimos anteriormente,

correspondía a los vínculos técnicos del problema.

3.3 DEFINICIÓN DEL PARÁMETRO Kmin Y SU FUNCIÓN.

En el programa que se desarrolló en JavaScript nos encontramos con un parámetro llamado k_min,

y que es preciso comentar ya que será una de las condiciones de parada del algoritmo.

Como sabemos, considerar lineal el problema de optimización objetivo de este trabajo es en

realidad una aproximación, ya que las relaciones entre las variables controlables y las variables del

problema no son lineales.

Si tenemos una función no lineal que tenga un solo mínimo en un determinado intervalo, determinar

éste mediante un método lineal, significa que, dado un punto inicial, tenemos que seguir la tangente

a la curva.

En el caso de un problema de mínimo con restricciones, como es nuestro caso, la optimización

lineal continuará por la tangente hasta encontrarse con una de las restricciones, y en ese punto se

parará y lo tomará como solución.



60

Si nos fijamos en la figura 2.2, el punto A representa las condiciones iniciales del problema, el

optimizador nos dirá que la solución se encuentra en el punto C, mientras que en realidad se

encuentra en el B.

Figura 3.3a Principio de resolución de un optimizador lineal.

De esta forma si es la solución ofrecida por el optimizador, podemos calcular un coeficiente que

la reduzca. Este coeficiente lo llamaremos kmin y será tal que:

El valor de este parámetro se calcula con un ciclo de optimización que se basa en el método de

Kiefer.

Método de Kiefer: dada una función objetivo estrictamente monótona en el intervalo de interés, es

posible eliminar partes sucesivas de un intervalo dado, conociendo el valor de la función en cuatro

puntos diversos del intervalo.

La elección de estos puntos puede efectuarse de diferentes formas, si usamos por ejemplo la serie de

Fibonacci es posible disminuir el intervalo a la longitud deseada en el menor número de iteraciones

posibles.

La función de optimización del kmin se implementa en JavaScript, refiriéndonos al esquema de la

Figura 3.3a, el cálculo se realiza después de la salida del optimizador lineal y antes de la

actualización del set point. El procedimiento de cálculo es el siguiente:

A

C

B

vincolo

F

minF



61

1. Acotamos 0 ≤ kmin ≤ 1, ya que estamos buscando un coeficiente que reduzca la

solución. Esto se hace definiendo: kup = 1, klow = 0.

2. Definimos algunas variables:

golden ratio: τ = 0.38197

tolerancia:

(3.3.1)

número máximo de iteraciones:

(3.3.2)

Que deberá ser aproximada a un número entero, en nuestro caso a 5.

El golden ratio es el factor que nos permite determinar los otros dos puntos del

intervalo, además de kup e klow , que nos sirven para calcular kmin.

(3.3.3)

(3.3.4)

El intervalo completo será:

3. Comenzamos con la primera iteración y calculamos la función objetivo utilizando

primero k1 y despues k2, en lugar kmin.

Si f( k1) >f(k2) entonces {

Esto querrá decir que el mínimo se encuentra a la derecha de k1.

Si f( k2) >f(k1) entonces {

kup k2 klow k1



62

El mínimo se encontrará a la izquierda de k2.

4. Con estos resultados se realiza una nueva iteración, y así sucesivamente hasta llegar

al número máximo de iteraciones, en nuestro caso 5. El resultado de las

iteraciones serán los valores actualizados de klow e kmin, es decir habremos acortado el

intervalo inicial.

En este momento se calcular el coeficiente kmin como el punto medio del intervalo

obtenido:

(3.3.5)

Una de las condiciones para salir del ciclo de optimizador es que kmin<0.05. Tomar como solución

una parte así de pequeña es como decir que el punto de mínimo ha sido encontrado en la iteración

precedente.

En el siguiente diagrama de bloques se explica el proceso de optimización:



63

LOAD FLOW TRIFASICO

(condiciones iniciales)

EMTP-RV

SI

- Valores

Iniciales:

Pset e Vnom

- Opciones de

Simulación.

- Parámetros del

problema de

óptimo.

MODELLO DELLA

RETE

Set point inalterados ¿Condición

Respetada?

NO

SOLUCION DEL PROBLEMA

DE MINIMO

Algoritmo de

recuperación de la

linealización

¿Condición

Respetada?

ACTUALIZACION DEL

SET POINT

SI

NO

Ciclo de optimización

Cálculo coeficientes de

sensibilidad

Matlab

Javascript

Load flow

EMTP-RV

Figura 3.3b Representación esquemática del proceso de optimización.

Javascript



64

La lógica del esquema es la siguiente:

Definición de la configuración:

- P de los generadores controlables.

- e .

- Opciones de simulación.

- Parámetros del problema.

Cálculo del Flujo de Cargas: A partir de los datos iniciales el programa EMTP_RV realiza

un primer cálculo de la repartición del flujo de potencia en estado de régimen permanente.

Si los resultados obtenidos respetan los vínculos se finaliza el proceso, devolviendo esos resultados

como óptimos. En caso contrario se activa el proceso de optimización.

Ciclo de optimización:

- Cálculo de los coeficientes de sensibilidad: linealizan el problema y son necesarios

para el funcionamiento del optimizador.

- Ejecución del algoritmo de optimización lineal con la función lp solve. Esta es la parte

fundamental del algoritmo (siendo además el punto más complejo de programar de todo

el proyecto). A partir de las condiciones iniciales la función lp solve devuelve en cada

iteración, una solución de óptimo considerando las hipótesis de simulación.

Además de minimizar la función según las hipótesis de simulación, el programa,

lpsolve22, desarrollado en este proyecto, devuelve un vector con los valores óptimos de

la potencia activa y reactiva de los generadores y los valores de regulación de los

transformadores.

- Corrección de los resultados obtenidos por el algoritmo de optimización lineal, a

través de un proceso de recuperación de la linealización mediante sucesivos cálculos de

LF (ver apartado 3.3).



65

- Se seguirán realizando nuevas iteraciones mientras no se respeten las condiciones

establecidas del problema.

Se ha realizado el programa para que durante las n-iteraciones vaya creando un hoja de cálculo en

Excel para poder visionar los resultados intermedios que se han ido obteniendo. Finalmente cuando

se cumple el criterio de convergencia se crea un nuevo archivo en Excel con los resultados finales.

Se han impuesto tres condiciones de parada del algoritmo iterativo:

La diferencia entre los valores de dos funciones objetivo sucesivas es inferior a un

determinado valor de tolerancia.

El parámetro

El número de iteraciones es mayor que 10.

Estos criterios de parada se han implementado en el programa lanzado directamente del EMPT-RV,

llamado nuovo31_dwj, que se adjunta en el Apéndice B.

.

FORMULACIÓN DEL PROBLEMA. LA FUNCIÓN LP...

Documents

Transcript of FORMULACIÓN DEL PROBLEMA. LA FUNCIÓN LP...