• Author / Uploaded
• Rafa Hernandez

Citation preview

Optimización y Simulación de Procesos Enrique Eduardo Tarifa Facultad de Ingeniería - Universidad Nacional de Jujuy

Introducción de Simulación Introducción El curso de simulación está basado en el libro “Chemical Engineering Dynamics” de Ingham et al. (1995) y en los manuales del HYSYS. Como se planteo anteriormente una parte importante de la simulación es el modelado del sistema en estudio. Kapur (1988) planteó treinta y seis principios del modelado, los más importantes son: • El modelo matemático sólo puede ser una aproximación del sistema real, el cual puede ser extremadamente complejo y aún no comprendido del todo. La complejidad del modelo estará determinada por los objetivos que persigue su construcción. El mejor modelo es el más simple que satisface las expectativas del ingeniero. • El modelado es un proceso continuo. Se debe comenzar modelando los fenómenos principales y luego ir agregando los restantes si es que son necesarios. Es decir, se va de menor a mayor, de un prototipo al modelo final. Tratar de enfrentar el problema de una sola vez puede traer problemas. • El modelado es un arte, no hay reglas fijas, pero el premio es un profundo conocimiento del sistema en estudio. Muchas veces, el problema que se quiere resolver con el simulador se resuelve en la etapa del modelado; es decir, antes de completar el simulador. • El modelo debe ser realista y robusto. Las predicciones del modelo deben estar de acuerdo con las observaciones y no debe ser demasiado sensible a cambios en los parámetros o variables de entradas. Los pasos para encarar las construcción de un modelo son: 1. Definición de los objetivos del modelo: Es de suma importancia establecer claramente qué es lo que se espera del modelo. Esto implica la definición de las fronteras del sistema a modelar, de la exactitud necesaria, de la necesidad de conocer la evolución en el tiempo, etc. 2. Formulación de un modelo físico: Se identifican los fenómenos más importantes que se llevan a cabo en el sistema. 3. Balances (estacionarios o dinámicos): Se escriben los balances de materia, energía, etc.. 4. Propiedades: Se obtienen o estiman los valores para las propiedades físico químicas que intervienen en los balances. 5. Suposiciones considerando los objetivos: Se realizan las simplificaciones en el modelo considerando los objetivos. 6. Consistencia matemática: El grado de libertad del modelo. 7. Se resuelve el modelo utilizando algún método matemático analítico o numérico (simulación). 8. Perfeccionamiento: Si los resultados obtenidos con el modelo no son satisfactorios, Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

1

entonces se repite el ciclo adoptando suposiciones menos restrictivas. En esta etapa se puede sugerir la adquisición de nuevos datos experimentales. El modelo puede ser obtenido empíricamente o utilizando los primeros principios. En el primer caso el tiempo de desarrollo es menor, no hay que pensar mucho, sólo se trata de ajustar parámetros para que las predicciones del modelo coincidan con los datos experimentales. Las redes neuronales son un ejemplo de este tipo de modelos. Sin embargo, la dificultad aparece a la hora de seleccionar los datos experimentales que se utilizarán para realizar el ajuste ya que se necesita una gran cantidad si es que se quiere que la zona de validez del modelo no sea muy pequeña. Por otra parte, el otro enfoque ofrece la ventaja de requerir menos datos experimentales, produce modelos más robustos; pero requieren de un gran esfuerzo para su desarrollo. Las herramientas que se disponen actualmente para asistir al ingeniero en el desarrollo de un simulador son muy variadas: 1. Lenguajes de programación: a. Fortran. b. Pascal. c. Delphi. 2. Lenguajes de simulación: a. GPROMS. b. ISIM. 3. Utilitarios matemáticos: a. MatLab. b. Simulink. c. MathCad. 4. Simuladores: a. Aspen. b. HYSYS. c. ChemCad.

Formulación de balances Un balances expresa el principio de conservación de una determinada propiedad en un dado volumen de control. Su formulación general es: {velocidad de acumulación} = {velocidad de entrada} + {velocidad de generación} - {velocidad de salida} Unidades: [{velocidad de acumulación}] = [Propiedad].[tiempo]-1 Los balances pueden ser: • Estacionario @ ({velocidad de acumulación} = 0). • Dinámico @ ({velocidad de acumulación} g 0). Se puede decir que Equilibrio < Estacionario pero no vale la inversa. Por ejemplo, en un reactor tubular funcionando en su estado estacionario no existe equilibrio dada la presencia de Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

2

gradientes de temperaturas, presiones y composiciones.

Ecuación de continuidad Para el caso de conservación de materia se debe considerar que no existe generación en el campo de aplicación de la ingeniería química, por lo tanto el balance debe tener la siguiente forma: {velocidad de acumulación} = {velocidad de entrada} - {velocidad de salida} Unidades: [{Acumula}] = [Masa].[tiempo]-1 Sólo se puede escribir una única ecuación de este tipo. Veamos la aplicación para un sistema de parámetros concentrados, un tanque de homogeneización:

d(ρV ) = F0 ρ 0 − Fρ dt

Figura 1: Tanque de homogeneización. Para la correcta formulación de los balances se puede seguir el siguiente procedimiento: • Elija como volumen de control una región donde las variables permanecen constantes o varían muy poco. En este ejemplo, se tomó como volumen de control el volumen líquido. • Identificar las corrientes que atraviesan las fronteras del volumen de control. En este caso son las corrientes de alimentación y de descarga del reactor; es decir que sólo existen flujos convectivos. • Escribir el balance con palabras, tal como se hizo al inicio de esta sección. • Expresar cada término en forma matemática utilizando variables medibles. A continuación se dan algunas formas matemáticas de los términos que generalmente aparecen en un balance: Velocidad de acumulación: d ( ρV ) (1) dt Flujo convectivo:

Fρ Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

(2) 3

Flujo difusivo: j i A = − 'i

dC i A dz

(3)

donde A es el área normal de flujo y Di es la difusividad del compuesto i. La fuerza impulsora es un gradiente de propiedad. Flujo en la interfaz:

J i = K i A ∆C i

(4)

donde A es el área total de transferencia, ûCi es la diferencia entre la concentración de equilibrio y la real del compuesto i (no es la diferencia entre las concentraciones de las fases) y Ki es el coeficiente global de transferencia.

Velocidad de generación:

R i = ri V

(5)

La velocidad de generación es positiva cuando i es un producto, y es negativa cuando se trata de un reactivo. Veamos un balance aplicado a un sistema con parámetros distribuidos, una tubería de transporte:

∂ ρ ∂ t

+

∂ ( ρv ) =0 ∂ z

Figura 2: Tubería de transporte.

Balance por componentes En forma general, este balance se puede expresar como: {velocidad de acumulación} = {velocidad de entrada} + {velocidad de generación} - {velocidad de salida} Unidades: [{velocidad de acumulación}] = [Masa].[tiempo]-1 Cantidad = NC - 1 (si es que se escribe el balance global) Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

4

A continuación se da un ejemplo con un reactor CSTR (Continuous Stirred Tank Reactor):

d(VC A ) = F0 C A 0 − FC A − VkC A dt

Figura 3: Reactor CSTR. Veamos un ejemplo aplicado a un sistema con parámetros distribuidos, un reactor tubular:

∂ CA ∂ t

+

∂ CA  ∂ ( vC A ) ∂   4 A  + kC A = ∂ z ∂ z ∂ z  Figura 4: Reactor tubular.

Otro ejemplo puede ese el que se presenta en la extracción con solvente desde un sólido. En este caso se desea disolver un sólido contenido en una matriz sólida como se muestra en la Figura 5.

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

5

Figura 5: Extractor sólido-líquido. Como se trata de un reactor batch, no existen salidas ni entradas durante su operación. En este caso, los flujos entre las fases son los relevantes. En palabras, el balance debe ser: {velocidad de acumulación en el solvente} = {velocidad de transferencia del sólido} Matemáticamente, el balance en la fase líquida es: dC L = k L A (C L* − C L ) VL dt

(6)

donde CL* es la concentración de equilibrio del lado líquido y es función de la concentración real del lado sólido, por lo tanto necesitamos plantear el balance también del lado sólido: dC S (7) = −k L A (C L* − C L ) VS dt donde CL* es función de CS.

Balance de energía El balance de energía es importante cuando la temperatura afecta las variables que nos interesan estudiar. Por ejemplo, en un reactor la conversión será afectada por la temperatura a través de la constante cinética; a su vez, la reacción afectará a la temperatura mediante la generación de calor. La ley de conservación de energía es planteado por el primer principio de la termodinámica. La energía interna depende no sólo de la temperatura sino también de la materia del sistema y su composición. Por este motivo, el balance de materia siempre está asociado al de energía. Para un sistema abierto que intercambia energía a través de su frontera, el balance de energía puede ser escrito como: {velocidad de acumulación} = {velocidad de entrada} - {velocidad de salida} Unidades: Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

6

[{velocidad de acumulación}] = [Energía].[tiempo]-1 Cantidad = 1 Algunos términos empleados a la hora de escribir el balance de energía son: • Formas Energía: potencial, cinética, interna. • Transporte de Energía: calor, trabajo. Tanto el calor como el trabajo no son energías sino formas de transporte. • “Generación”: La energía no se genera (por lo menos en las condiciones de las plantas químicas), el calor de reacción tiene otro origen que será analizado más adelante. Con estas aclaraciones, podemos escribir matemáticamente el balance de la siguiente forma: N N dE = F0 ∑ C i 0 E i 0 − F1 ∑ C i 1 E i 1 + Q − W dt i =1 i =1

(8)

donde E es la energía total del sistema Ei es la energía total por mol de compuesto i, F es el caudal volumétrico, C es la concentración molar, N es el número total de compuestos incluyendo los inertes, Q es el calor suministrado al sistema, mientras W es el trabajo entregado por el sistema. El término de trabajo puede ser separado considerando el trabajo de flujo y el trabajo neto WS: N

N

i =1

i =1

− W = F0 ∑ C i 0 P V i 0 − F1 ∑ C i 1 P V i 1 − W S

(9)

La energía E es la suma de la energía interna U, la energía cinética, la energía potencial y toda otra energía. De todos ellas, la energía interna es la predominante en los procesos químicos comunes. Entonces, haciendo E U y recordando que la entalpía se define como H = U + PV, y suponiendo el trabajo entregado igual a cero, podemos rescribir el balance de energía de la siguiente forma: d  V dt 

N N  = − C h F C h F  ∑1 i i 1  0 ∑1 i 0 i 0 1 ∑1 C i 1 hi 1 + Q − W s i= i= i= N

(10)

La ecuación anterior vendría a ser el “balance térmico” o de “energía química” (energía interna, calor), donde se consideran sólo los efectos térmicos. Otro balance que suele plantearse es el de energía mecánica pura (potencial, cinética y trabajo). En rigor, debería plantearse un único balance, pero por conveniencia se plantean en forma separada vinculados por términos de intercambios, por ejemplo el término de disipación viscos que transforma la energía mecánica en energía química, y el término de trabajo útil que transforma la energía térmica en trabajo mecánico.

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

7

El gran problema con el balance de energía es la determinación de las entalpías. Los valores absolutos de las mismas son imposibles de medir (lo mismo ocurre con la energía potencial y cinética); lo que sí se pueden medir son los cambios de entalpías. Es así como aparecen los estados de referencia y los ûh. Generalmente se utilizan las capacidades caloríficas específicas para calcular los cambios de entalpía como sigue:

∆h = h1 − h 0 =

∫T C p ( T ) d T T1

(11)

0

Mientras es común encontrar en la bibliografía los coeficientes para la siguiente expresión de cálculo de los Cp: (12) C p (T ) = a + b T + c T 2 El calor de reacción Note que en las ecuaciones anteriormente escritas no aparece el término de calor de reacción. Este término no existe físicamente porque la energía no se genera. El calor que genera la reacción es realidad surge de una transformación de energía que ya venía en los reactivos. A continuación veremos entonces cómo es que aparece el término de calor de reacción. En primer lugar, rescribiremos el término de acumulación como sigue:

d  V dt 

 ∑1 C i hi 1  = i= N

d (V C i 1 ) N dh ∑1 hi 1 d t + ∑1 V C i 1 d ti 1 i= i= N

(13)

Por otra parte, tenemos el balance de materia para cada componente:

d (V C i 1 ) = F 0 C i 0 − F1 C i 1 + ri V dt

(14)

Multiplicando cada una de estas ecuaciones por su correspondiente hi1 y sumando para todos los componentes, se tiene: N N d (V C i 1 ) ∑1 hi 1 dt = F0 ∑1 hi 1 C i 0 − F1 ∑1 hi 1 C i 1 + V i= i= i= N

N

∑r h i =1

i

i1

(15)

Recordando que ri = .i rg, donde rg es la velocidad de reacción general y .i es el coeficiente estequeométrico del compuesto i (positivo para los productos y negativo para los reactivos), tenemos que el último término es el calor de reacción a la temperatura T1, esto es:

V

N

N

i =1

i =1

∑ ri hi 1 = V r g ∑ α i hi 1 = V r g ∆ H

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

(16)

8

Restando la ecu.(15) de la ec.(10), considerando todas las ecuaciones intermedias, y suponiendo R reacciones, tenemos el balance: N

N d hi 1 = F0 ∑ C i 0 ( h i 0 − h i 1 ) + V V ∑ C i1 dt i =1 i =1

R

∑ rg j ( − ∆ H j 1 ) + Q − W s

(17)

j =1

Trabajando con capacidades caloríficas tenemos: N  N  dT V  ∑ C i 1 C p i 1  1 = F0 ∑ C i 0  i =1  dt i =1

∫T

T0

C p i d Ti + V

1

R

∑ rg j ( − ∆ H j 1 ) + Q − W s

(18)

j =1

A continuación se verán otras maneras empleadas para escribir los términos de los balances de energía. Término de acumulación Para cambios moderados de temperaturas, los Cp pueden ser considerados constantes, por lo tanto: V ρ1 C p 1

dT1 dt

(19)

donde !1 es la densidad en el sistema. Términos de flujo Con las mismas consideraciones, tenemos: F 0 ρ0 C p ( T0 − T1 )

(20)

U A ( Ta − T1 )

(21)

Términos de transferencia

donde Ta es la temperatura ambiente, U es el coeficiente global de transferencia y A es el área de transferencia. Términos de reacción La definición de la velocidad de reacción para el compuesto i es: ri =

1 dn i V dt

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

(22)

9

donde ni es el número de moles del compuesto i. También tenemos la velocidad de reacción general que es independiente del compuesto:

ri αi

rg =

(23)

donde .i es el coeficiente estequeométrico del compuesto i. La velocidad de una reacción irreversible puede ser expresada como: NR

rg = k ∏C iβ

(24)

i

i =1

donde k es la constante cinética y es calculada por la fórmula de Arrhenius: k = Z e−

E RT

(25)

donde Z y E (energía de activación) son determinados experimentalmente. Otros términos Dependiendo del caso en particular, pueden tener efectos importantes el calor generado por la agitación, pérdidas de calor por convección, evaporación, radiación o mezclado. Balance simplificado Suponiendo constante el Cp, se tiene:

V ρ1 C p

d T1 = F0 ρ0 C p ( T0 − T1 ) + V dt

R

rg j ( − ∆ H j 1 ) + U ∑ j 1 =

A ( Ta − T1 ) − W s

(26)

A continuación se muestra un balance en un tanque de calefacción:

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

10

d(ρVU ) = F0 ρ 0 h0 − F ρ h + Q dt

ρC p

d (V T ) = ρ C p ( F 0 T0 − F T ) + Q dt

F0,!0,T0 Q

V,! T

F,!,T Figura 6: Tanque calefaccionado. El balance simplificado por la combinación del balance de materia es:

V ρ Cp

dT = F 0 ρ 0 C p ( T0 − T ) + Q dt

(27)

Otro ejemplo, es el balance aplicado a un sistema de dos fases: d( ρ vVv H + ρVL h) = F0 ρ 0h0 − Fρh − Fv ρ v H + Q dt d[ ρ vVv (C p T + λ v ) + ρVL C p T ] dt F0,!0,T0

= F0 ρ 0C p T0 − FρC p T − Fv ρ v (C p T + λ v ) + Q Vv,!v,Tv,P VL,!,T

Fv,!v,Tv

F,!,T

Q

Figura 7: Evaporador. Veamos a continuación un ejemplo sobre un sistema con parámetros concentrados, un reactor CSTR calefaccionado:

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

11

d(ρVU ) = F0 ρ 0h0 − F ρ h + Q − λ V k C d t

ρC

p

d (V T ) = ρC d t

p

A

( F0 T0 − F T ) + Q − λ V k C

A

F0,!0,CA0,T0

V,! CA,T

Q

F,!,CA,T

Figura 8: Reactor CSTR. Ahora veamos el caso de un sistema con parámetros distribuidos, un tubo calefaccionado:

∂ ( ρC p T ) ∂ t

+

∂ (vρC p T ) 4 hT ∂  ∂ T  [T  + ( T − TM ) = ∂ z ∂ z  ∂ z  D

Figura 9: Tubo calefaccionado.

Conservación de cantidad de movimiento La forma general es: {velocidad de acumulación de cantidad de movimiento} = {velocidad de entrada de c.m } - {velocidad de salida de c.m.} + {fuerzas} Cantidad: 1. La ecuación es de carácter vectorial. Su estudio es todo un campo, el de la Mecánica de los Fluidos, que está fuera del alcance de este curso. Sólo se verán algunos ejemplos para dar una idea general. La expresión matemática simplificada de este balance es: &

& dM v & & = F 0 ρ0 v 0 + F1 ρ1 v 1 + ∑ F dt

(28)

Veamos un ejemplo para un tubería, el balance correspondiente es:

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

12

dv 1 K = ( P0 − P1 ) − F v 2 ρA d t Lρ

Figura 10: Tubería.

donde se modeló el término de fricción como proporcional al cuadrado de la velocidad del fluido y a la longitud del tubo. El balance para un sistema de parámetros distribuido es:

∂ vz

=

∂ t

µ 1 ∂  ∂ v z  ∆P r − ρ r ∂ r  ∂ r  ρL

Figura 11: Tubería.

Ecuaciones complementarias Junto con los balances se necesitan ecuaciones complementarias como ya vimos. A continuación se listan las más comunes. Ecuaciones de Transporte Forma General {Flujo}=-{Constante}{Gradiente} Unidades: [{Propiedad}].[tiempo]-1.[longitud]-2 Ejemplos: • Ley de Furier (q). • Ley de Fick (NA). • Ley de Newton (2rz). • Moleculares (Gradientes). • Globales (Diferencias). Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

13

Ecuaciones de Estado Ejemplos: • !L = f(P,T,x) • !v = f(P,T,y) • h = f(P,T,x) • H = f(P,T,y) Ecuaciones de Equilibrio Ejemplos: • Químico. • Entre fases: • Ley de Rault. • Ley de Henry. • Volatilidad Relativa (.ij). • Valores K. • Coeficiente de actividad. Ecuaciones de Cinética Ejemplos: • Ecuación de Arrhenius. • Ley de acción de masas.

Optimización y Simulación de Procesos. Simulación Introducción Enrique Eduardo Tarifa - Universdad Nacional de Jujuy

14

Optimización y Simulación de Procesos Enrique Eduardo Tarifa Facultad de Ingeniería - Universidad Nacional de Jujuy

Métodos Analíticos La Ingeniería y la Optimización "Para definirlo ruda pero no inapropiadamente, ingeniería es el arte de hacer bien con un dolar aquello que cualquier chapucero haría con dos dólares" ("To define it rudely but not inaptly, engineering is the art of doing that well with one dollar which any bungler can do with two dollars.") -Arthur M. Wellington, The Economic Theory of the Location of Railways, Introduction (6th ed., 1900). La optimización está en todo lo que el ingeniero hace. Cuando alguien propone la construcción de un túnel, el ingeniero aporta diciendo que se inicie la construcción desde los dos extremos para ahorrar tiempo y reducir el recorrido de los materiales de deshecho. Una de las mayores diferencias entre la Ingeniería con respecto a la Ciencia y la Tecnología es precisamente su sentido de eficacia y eficiencia. Mientras la Ciencia busca conocer, y la Tecnología busca resolver problemas a cualquier costo, la Ingeniería usa la Ciencia para perfeccionar las soluciones obtenidas por la Tecnología, a la vez que reduce sus costos y potenciales peligros. Cualquiera puede hacer Tecnología, no es necesario ser un científico ni un ingeniero, Edison es un buen ejemplo de ello. Sin embargo, la ingeniería sólo puede ser hecha por una persona con sólidos conocimientos teóricos y un sentido de lo factible y práctico muy desarrollado. Si la Ciencia es el arte de hacer posible lo imposible, la Ingeniería es el arte de hacer práctico lo posible (Luyben, 1973). La formación teórica es lo que distingue a un ingeniero de un técnico. Podemos así definir una cadena de evolución del conocimiento, éste se inicia en la Ciencia o en la Tecnología, y evoluciona hasta llegar a la Ingeniería donde se vuelve de utilidad general para toda la sociedad. No siempre es el mismo grupo de investigación, ni siquiera el mismo país, el que realiza la cadena completa. Tomemos como ejemplo la electrónica, los europeos dieron los primeros pasos de esta ciencia, pero fueron los americanos los que desarrollaron la tecnología con aplicaciones como las válvulas y los transistores. Sin embargo, los Japoneses hicieron la ingeniería, “optimizaron” el tamaño y el costo, y hoy son líderes en electrónica. En Química, fueron nuevamente los europeos los que desarrollaron la ciencia y la tecnología. El precio que pagaron por este desarrollo a cualquier costa fue la contaminación ambiental, la pérdida de recursos naturales, conflictos sociales, etc. Todavía hoy seguimos pagando esos costos, pero estamos trabajando en la ingeniería que resolverá esos problemas.

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

1

Un Ingeniero debe estar capacitado para reconocer los problemas y encarar su solución con eficacia y eficiencia. Es en la eficiencia donde está involucrada la optimización; por ejemplo, suponga que tenemos una corriente para enfriar y otra para calentar, cualquiera puede proponer que calentemos la corriente fría con vapor, y que enfriemos la caliente con agua de enfriamiento. En cambio, el ingeniero analizará la posibilidad de combinar ambas corrientes en un intercambiador para minimizar el uso de los servicios (tecnología Pinch). Veamos un ejemplo, suponga se quiere diseñar el sector de producción de NH3 de una planta. Se dispone del gas de síntesis desde otro sector. Lo que un técnico nos diría es que sólo necesitamos un reactor y tal vez un separador. En cambio, un ingeniero vería que mucho de los reactivos quedan sin reaccionar, por lo tanto propondría un reciclo. Debido a que también hay inertes, el ingeniero sabe que deberá purgar algo del reciclo. Además, nota que tiene corrientes que deben calentárse y otras que deben enfriarse, por lo tanto propone una integración

Figura 1: Planta de NH3. energética. Se da cuenta también que existe un reciclo y una purga óptima. Nuestro ingeniero no se queda tranquilo, y continúa con el análisis de los materiales, del diseño de los equipos, de la operabilidad, etc., siempre buscando la mejor solución.

Un ingeniero debe resolver en forma óptima los problemas planteados en los siguientes niveles: Al nivel de gerenciamiento: Evaluación de un dado proyecto, comparación entre 1. proyectos, asignación de capitales a diferentes proyectos. Nuevos proyectos: Elección del producto, elección de la capacidad de la planta, 2. elección de la ruta química, elección de las materias primas, elección de la localización. Process and flowsheeteing: Selección del proceso, selección del flowsheet, 3. ordenamiento de los flowsheet. Especificaciones de los equipos: Selección del tipo de equipos para un dado flowsheet, 4. selección de los materiales de construcción, dimensionamiento de los equipos, layout. Condiciones operativas: Unidades de procesos dependientes e independientes del 5. tiempo, sistemas de control, política de mantenimiento y reemplazos de equipos, scheduling de la producción, localización de materias primas, stock de materias primas, de productos intermedios y productos; despacho, transporte, y distribución.

Definición del problema Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

2

Hasta aquí, tenemos una idea de lo que es Optimización, pero hasta ahora no la hemos definido. Podemos decir que Optimización significa seleccionar el “mejor” curso de acción de las posibles “alternativas”. Esto nos plantea el problema de definir qué es lo “mejor” y cuáles son las “alternativas”. Para seleccionar lo mejor, primero hay que probar y comparar. Para probar, necesitamos una metodología que nos permita encontrar lo mejor con la mínima cantidad de ensayos considerando las alternativas disponibles. Para comparar, necesitamos tener alguna forma de ponderar la bondad de una dada solución. Esquemáticamente, podemos plantear el problema como sigue: Maximizar {Objetivo} Alternativas

s.a (sujeto a): Alternativas 0 {posible, práctico}

La solución Los dos elementos que componen la solución del problema de optimización son: • La Alternativa Optima: es la Alternativa que cumple con todas las restricciones y que a la vez genera el valor máximo del Objetivo. • El Valor Optimo: es el valor que adopta el Objetivo cuando se aplica la Alternativa Optima.

Función objetivo Es la función matemática que se empleará para ponderar la bondad de una solución; podemos emplear distintos criterios: • Económicos: Beneficios, costos, retorno de capital, TIR, VAN. • Técnicos: cantidad de producto, tiempo de producción, servicios utilizados, calidad del producto. • Seguridad: probabilidad de accidentes, consecuencias de los accidentes. • Ambiental: contaminación, agotamiento de recursos naturales. La solución del problema de optimización depende del criterio que se adopte en la función objetivo. Muchas veces es necesario trabajar con varios criterios a la vez y no es raro que estos criterios compiten entre sí. En tales casos se tienen tres opciones: 1. Plantear un problema con múltiples objetivos: se plantea una función objetivo para cada criterio. Este campo está fuera del alcance del presente curso. 2. Plantear una función objetivo combinada: se construye una función objetivo que combina los criterios a seleccionados. Por ejemplo, el criterio ambiental puede ser combinado con uno económico si se tiene en cuenta las multas que pagaría la empresa por arrojar efluentes que violen los límites especificados por la legislación del lugar de Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

3

3.

operación. Agregar los objetivos secundarios como restricciones: se selecciona un criterio principal para la función objetivo, mientras los restantes criterios se utilizan para plantear restricciones. Por ejemplo, un criterio de seguridad podría limitar la temperatura que puede alcanzar una reacción.

Restricciones Se debe observar que las alternativas a explorar están sujetas a un conjunto de restricciones. Estas restricciones surgen tanto de la naturaleza propia del sistema como de los criterios ingenieriles aplicables. Las primeras, son las que constituyen el modelo del sistema. En efecto, el sistema plantea relaciones entre las variables que no pueden ser ignoradas al plantear las alternativas, por ejemplo: la velocidad de reacción aumenta con la temperatura. Las otras restricciones surgen del criterio del ingeniero, ya sea para acotar la solución o para facilitar la solución del problema. Entre las primeras se puede citar: el tanque deberá ser cilíndrico porque es más fácil de construir, el producto debe tener determinada concentración, la presión no podrá superar un dado valor porque sino se necesitará de otra tecnología. Entre las segundas se puede citar, la composición será una variable positiva, el reciclo será menor que el caudal de salida del reactor, etc.

Alternativas Antes de resolver el problema, primero el ingeniero tiene que detectarlo. Es fundamental que el ingeniero sepa cuándo es posible optimizar y cuándo no. Por ejemplo, si nos plantean que dibujemos una línea recta que pase por dos puntos dados, no hay alternativas. Sin embargo, si sólo nos dan un punto, existirán infinitas soluciones y se presentará la oportunidad para optimizar. Primero, se deberá elegir la función objetivo; supongamos que el área bajo la recta representa la energía consumida en el proceso, analizando la misma descubrimos que la misma varía con la pendiente de la recta y que existe un mínimo. Entonces, elegiremos la pendiente que minimice el área bajo la recta. Supongamos ahora que nos dan diez puntos no alineados. No existe recta que pase por todos ellos. Sin embargo, podemos optimizar; ya que no es posible tocar todos los puntos, por lo menos podemos pasar lo más cerca posible de todo ellos (regresión). En resumen, podemos optimizar cuando el sistema tiene alguna variable libre o cuando el problema está sobredeterminado; en definitiva, cuando el problema no tiene una única solución o no tiene solución exacta.

Un caso de estudio Dado una tubería de vapor, deseamos determinar óptimamente el espesor del aislante. Podemos hacer esto gráficamente considerando el costo total CT (costo fijo y variable). Esto nos dará un punto óptimo sin restricciones expresas para la variable espesor. El costo fijo CF estará dado por el costo de instalación del aislante, mientras que el costo variable CV estará dado por el costo del gas empleado para generar el vapor. Estos costos deben ser comparables, por lo tanto se Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

4

adopta una base anual. La Figura 2 muestra la determinación del punto óptimo en forma gráfica correspondiente a un espesor óptimo xo con un costo mínimo CTo.

Costos de aislación 15

CVx () [$/año]

CFx ( ) 10 CTx () CTo

5

0

0

2

4 x [cm]

Figura 3: Espesor óptimo.

Antes de tomar esta solución, un ingeniero plantea que según la normas de seguridad, la temperatura externa no debe superar los 60°C. Esto obliga a aumentar el espesor del aislante a un valor mínimo xT; es decir, se agrega la restricción x  xT, la solución es ahora la mejor que se puede obtener dada la restricción de seguridad. Además, puede existir un problema de espacio donde la distancia a la pared hace que el espesor no pueda superar el valor máximo xP; es decir, se agrega la restricción x  xP Esta es una restricción física que también puede afectar el valor óptimo. Con estas restricciones es interesante analizar los casos que podrían presentarse en cuanto a qué restricción se vuelve activa (hace valer la igualdad): xT < xo < xP: En este caso, ambas restricciones son inactivas, y es como si no existieran. 1. xP < xo: La restricción activa que representa la pared hace que tengamos que renunciar 2. al óptimo xo y conformarnos con xP. xo < xT: La restricción activa que representa la temperatura mínima en la pared hace que 3. tengamos que renunciar al óptimo eo y adoptar el valor eT. eP < eT: El problema no se puede resolver. 4.

Magnitud del Problema Los problemas que enfrenta el ingeniero son de una magnitud tal que no siempre es posible resolverlos en su totalidad, por ello se plantean técnicas para descomponer el problema, en otros Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

5

más pequeños y manejables; sin embargo, la solución óptima de los subproblemas no siempre está de acuerdo con la solución del problema total. Analicemos un ejemplo, deseamos separar ABCD, ordenados de acuerdo a sus volatilidades relativas, A es el más volátil. Si sólo nos restringimos a la operación de destilación, existen cinco alternativas posibles; en la Figura 3 se muestran sólo dos de ellas. Para decidir cuál de todas las alternativas es la mejor es necesario avanzar con el diseño de cada una hasta determinar el más mínimo detalle. Esto es, será necesario analizar la instrumentación y control óptimo, la forma de operación óptima, analizar la integración energética, etc. Esto implica explorar un árbol hasta llegar a sus hojas finales para poder compararlas entre ellas y recién elegir la mejor. Esta tarea es monumental, y aún no fue totalmente resuelta. Lo más que se puede hacer actualmente es particionar el problema, por ejemplo, se optimiza el diseño de los equipos por un lado, el control por otro, la puesta en marcha por otro, etc. Sin embargo, como ya se mencionó, de esta manera no se asegura encontrar la solución óptima para el sistema en su conjunto.

Figura 3: Síntesis de procesos.

Formulación Es sorprendente que todos los problemas de optimización tengan la misma forma, independientemente del área en que se esté trabajando. La formulación matemática general es (Ray y Szekely, 1973):

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

6

Max I(x1,x2,...,x n) i 1,2,...,n dx lx i l 1,2,...,m3 f l(x1,x2,...,x n) dt x l(t0) xl0 t0  t  t f

(1) (4)

h j(x1,x2,...,xn)

0

j

1,2,...,m1

(2)

g k(x1,x2,...,xn)  0

k

1,2,...,m2

(3)

La ec. 1 es la Función Objetivo (FO), es la que se quiere maximizar. El problema es determinar cuáles son los valores de los xi que maximizan la FO. Aunque existen trabajos que consideran más de una FO, en este curso nos restringiremos a estudiar problemas con un único objetivo. Las ecs. 2-4 son las restricciones del problema. Las ecs. 2 son las restricciones algebraicas. Cada una de estas restricciones reduce en una unidad el grado de libertad del problema; esto es, cuántas de las xi pueden variar independientemente. Las ecs.3 son las restricciones de desigualdad. Estas ecuaciones no afectan el grado de libertad pero pueden afectar a la bondad de la solución. Es siempre interesante analizar cuáles de todas las restricciones se vuelven activas en el óptimo (hacen valer la igualdad). De este análisis se puede tomar decisiones a fin de dirigir futuras inversiones para relajar las restricciones activas más importantes. Las ecs. 4 son las restricciones diferenciales o integrales de igualdad. La región determinada por todos los puntos que satisfacen todas las restricciones se denomina región factible. La búsqueda del óptimo debe restringirse a esta región. Tradicionalmente estos problemas fueron clasificados en: • Programación matemática. • Optimización de trayectoria.

Programación matemática Es el problema que no tiene restricción alguna, o sólo tiene del tipo 2 y 3, es decir: i

Max I(x1,x2,...,x n)

1,2,...,n

xi

h j(x1,x2,...,xn)

(5)

0

j

1,2,...,m1

(6)

g k(x1,x2,...,xn)  0

k

1,2,...,m2

(7)

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

7

Generalmente, para facilitar la resolución del problema, se suele agrupar las variables xi en variables de estado y variables de control o decisión. Esto surge de las restricciones de igualdad, como tenemos m1 ecuaciones de este tipo, frecuentemente es posible despejar una variable de cada una de dichas ecuaciones en función de las n-m1 restantes, es decir: xj

F j(xm 1,xm 2,...,xn) 1

j

1,2,...,m1

1

(8)

Hecho esto, podemos definir los siguientes vectores:

y

u

y1

x1

y2

x2

...

...

ym

xm

1

(9)

1

u1

xm

1

1

u2

xm

1

2

... un m

(10)

... xn

1

donde y es el vector de las variables de estado, mientras u es el vector de las variables de decisión. Luego, la ec. 8 se transforma en: yj

j

F j(u1,u2,...,un m ) 1

1,2,...,m1

(11)

Finalmente, el problema puede ser reformulado como: Max I(y1,y2,...,y m ,u1,u2,...,un m ) 1

ui

yj

i

F j(u1,u2,...,un m )

j

1

g k(y1,y2,...,y m ,u1,u2,...,un m )  0 1

1,2,...,n m1

1

1

1,2,...,m1

k

1,2,...,m2

(12)

(13)

(14)

Esta formulación es más fácil de resolver que la primera. Sin embargo, para problemas grandes Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

8

no siempre es posible despejar las variables, por lo que debemos resolver la primera formulación.

Optimización de Trayectoria Agregan las restricciones del tipo integral o diferencial. Analicemos el caso en que sólo está presente este tipo de restricciones (sin las de igualdades ni las de desigualdades): Max I(x1,x2,...,x n) t t x i(t)

dxl dt

i

f

l

f l(x1,x2,...,x n) xl0

x l(0)

1,2,...,n

1,2,...,m3

t0  t  t f

(15)

(16)

Como no hay problemas en despejar las variables de estado, siempre se reformula el problema como sigue: Max I(y1,y2,...,y m ,u1,u2,...,un m ) t t u i(t)

dy l dt

3

3

f

fl(y1,y2,...,y m ,u1,u2,...,un m ) 3

yl(t0)

3

yl0

i

1,2,...,n m3

l

1,2,...,m3

t0  t  t f

(17)

(18)

Clasificación de los modelos De acuerdo a la naturaleza de la FO y las restricciones, los modelos se clasifican en (Scenna et al., 1999): 1. Programación Lineal (LP, Linear Programming): tanto la FO como todas las restricciones son lineales. 2. Programación No Lineal (NLP, Nonlinear Programning): al menos una de las funciones anteriores es no lineal. 3. Programación Lineal Entero (ILP, Integer Linear Programming): el modelo es lineal pero las variables son enteras. 4. Programación No Lineal Entero (INLP, Integer Nonlinear Programming): el modelo es no lineal pero las variables son enteras. 5. Programación Lineal Entero Mixto (MILP, Mixed Integer Linear Programming): el modelo es lineal pero involucra variables continuas y enteras mezcladas. 6. Programación No Lineal Entero Mixto (MINLP, Mixed Integer Nonlinear Programming): el modelo es no lineal pero involucra variables continuas y enteras mezcladas. Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

9

Planteo y solución del problema Lo primero que debe hacer el ingeniero es detectar si la situación presente permite plantear un problema de optimización, algunas de estas situaciones son: 1. Ventas limitadas por producción: se debe maximizar la producción. 2. Ventas limitadas por el mercado: se debe minimizar los costos de producción. 3. Grandes volúmenes de producción: pequeños ahorros en los costos de producción tienen un gran impacto, es el caso de las destilerías y plantas químicas. 4. Alto consumo de materia prima o energía: optimizar los equipo con mayor consumo. 5. La calidad de los productos es superior a la requerida por el mercado: minimizar los costos reduciendo la calidad hasta el límite de aceptación. 6. Pérdida de componentes valiosos en los efluentes: minimizar las pérdidas teniendo en cuenta la legislación ambiental. 7. Alto costo de mano de obra: pasar de batch a continuo, modificar el scheduling, etc. Dos reglas prácticas para identificar la oportunidad de plantear un problema de optimización son: • Descubrir efectos que se opongan entre sí (trade-off); por ejemplo, el costo de instalación de aislante en una tubería de vapor varía en forma inversa al costo de gas utilizado para producir vapor. Entonces, debe existir una situación intermedia donde los dos efectos no son tan desfavorables. • Descubrir una variable cuyo valores extremos conducen a situaciones inconvenientes; por ejemplo, instalar una capa de aislante con un espesor de 1 m reducirá los costos de operación pero elevará sensiblemente los costos de instalación; por otra parte, un espesor de 1 mm reducirá los costos de instalación pero elevará los costos de operación. Entonces, debe existir una situación intermedia que de lugar a un valor óptimo. El procedimiento general para plantear y resolver el problema es: 1. Examine el proceso y haga una lista de las variables y las características de interés. 2. Determine el criterio de optimización; esto es, construya una función objetivo que involucre las variables del punto anterior. 3. Construya el modelo del proceso; esto es, el conjunto de funciones matemáticas que relacionan las variables del punto anterior. Este conjunto serán las restricciones. Identifique las variables independientes y dependientes para determinar el grado de libertad. 4. Si la formulación del problema es muy grande: 5. Particione el problema y/o 6. Simplifique el problema. 7. Aplique un adecuado método matemático para resolver el problema. 8. Verifique la solución, analice la sensibilidad de la solución con respecto a los parámetros y suposiciones. Nunca se debe olvidar que se está resolviendo el modelo planteado y no el sistema real. Esto significa que si el modelo está mal, por más que ejecutemos correctamente los restantes pasos, la solución no será adecuada. Esto se debe tener en cuenta a la hora de hacer simplificaciones. Supongamos que determinamos el perfil óptimo de temperatura para un reactor tubular (Figura Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

10

4 ), la implementación del sistema de calefacción podría ser tan compleja que haría perder todas las ventajas logradas; en cambio, se podría implementar dos reactores tubulares isotérmicos que aproximen el perfil óptimo sin tanto problemas y con un efecto despreciable sobre el rendimiento final. Esta situación se presentó debido a que, por una simplificación, la función objetivo no contempló los costos de implementación ni tampoco el modelo contempló más de un reactor, sólo se consideró el avance de la reacción.

Figura 4: Perfil de temperatura en un reactor tubular.

Ejemplos Fábrica de contenedores En este ejemplo analizaremos los seis pasos de resolución explicados en la sección anterior. Tenemos una planta que produce contenedores de plástico, la demanda anual es de Q = 100000 contenedores distribuida homogéneamente a lo largo del año, y el precio se mantiene fijo. El problema es determinar el schedule de producción (Edgar y Himmelblau, 1988). Paso 1: Una opción es producir todo Q al principio del año; pero el costo de inventario sería enorme. Este costo se reduciría si en cambio produjéramos Q en varias operaciones distribuidas en el año. Sin embargo, para el extremo de encender y apagar la planta cada hora, los costos de producción son inaceptables. Por lo tanto, debe haber un óptimo. Las variables que parecen ser importantes son: D: el número de unidades producida en cada corrida, es un número entero. n: el número de corridas por año, no es necesariamente entera. Q: el total de unidades producidas por año. Resta hacer un análisis de los costos para plantear el modelo e identificar las variables de decisión. Paso 2: Existen dos costos: el de inventario y el de producción. El costo de inventario incluye los costos Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

11

de productos y de almacenaje, y puede ser modelado como K1 D. Por otra parte, el costo de producción puede ser modelado como K2 + K3 D (costo fijo mas costo variable). Entonces, el costo total anual C es: C

K1 D

n (K2

K3 D)

(19)

Paso 3: Tenemos la siguiente restricción de igualdad: n

Q D

(20)

Esta restricción reduce un grado de libertad. Podemos elegir D como la variable independiente y n como la dependiente, siempre conviene elegir la variable entera como variable independiente. Debemos recordar que D está restringida a valores enteros y positivos. Eliminando n de la función objetivo, obtenemos: C

K1 D

K2 Q D

K3 Q

(21)

Este paso tiene la ventaja de eliminar la única restricción explícita que existía. El mismo método se puede aplicar para eliminar múltiples restricciones de igualdad.

Paso 4: No necesario. Paso 5: Este problema se puede resolver analíticamente. Diferenciando con respecto a D e igualando a cero tenemos: dC dD

K1

K2 Q D2

0

(22)

despejando D, tenemos: D opt

K2 Q K1

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(23)

12

Note que si los valores de los parámetros cambian, la solución óptima se obtiene simplemente recalculando la ec. 23, ésta es una clara ventaja del método analítico sobre los numéricos. Podemos verificar las condiciones de suficiencia para un mínimo: d 2C

2 K2 Q

dD 2

D3

> 0

(24)

Otra ventaja de la solución analítica es que permite formular reglas empíricas. En este caso, note que K3 no afecta la determinación del óptimo (ec. 23). También, note que si estamos en una schedule óptimo y deseamos cambiar a otro valor de Q, no es necesario conocer los valores de los parámetros para establecer el nuevo óptimo, basta con saber que varía con la raíz cuadrada de Q, esto podría plantearse como un eurístico. Paso 6: Finalmente, debemos analizar la sensibilidad de la solución. Dado que los parámetros, y el modelo, no corresponden totalmente a la realidad, es necesario estimar cual será el impacto de estos errores sobre la solución. Para este ejemplo podemos realizar el análisis analíticamente. Primero, sustituimos la solución en la función objetivo: C opt

2 K1 K2 Q

K3 Q

(25)

y luego calculamos la sensibilidades absolutas a través de las derivadas parciales: 0C opt 0K 1

K2 Q

0C opt 0K 2

K1 Q

(27)

K2

0C opt 0K 3

0C opt 0Q

(26)

K1

Q

K1 K2 Q

(28)

K3

(29)

Lo mismo podemos hacer para D: Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

13

opt 0D opt 0D 1 0 K2 Q K1 0K1 0K32 K1

(30) (32)

Q KK22 Q 11 22 K KK11 Q2

(31) (33)

0D opt 0D opt 0K 2 0Q

De mayor utilidad son las sensibilidades relativas, a modo de ejemplo: C

S K1

0C opt/C opt 0K1/K1

0ln(D opt) 0ln(K1)

K2 Q

K1

K1

C opt

(34)

Reconciliación de lecturas I Sea el sistema de la Figura 5, en un momento dado se tiene las siguientes lecturas: S1 [m3/h]

S2 [m3/h]

S3 [m3/h]

10

15

27

Figura 5: Mezclador. Aparentemente, el balance de materia no cierra. Sin embargo, esto se debe a errores en las lecturas; pero ¿a qué sensor debemos creer? A ninguno, se debe plantear el siguiente problema de optimización:

M (F 3

Min

F1,F2,F3 i 1

F1

F2

i

Si)2

F3

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(35)

(36)

14

Fj  0

j

1,2,3

(37)

donde los Fi son los caudales verdaderos, mientras que los Si son los indicados por los sensores. Para este caso, la solución es: F1 [m3/h]

F2 [m3/h]

F3 [m3/h]

10.67

15.67

26.34

Reconciliación de lecturas II Veamos un caso más complejo. En la Figura 6 se muestra un sistema de tuberías con caudalímetros Si. Cada caudalímetro tiene un error menor a 1.5 m3/h. Si las mediciones en un momento dado son: S1 [m3/h]

S2 [m3/h]

S3 [m3/h]

S4 [m3/h]

S5 [m3/h]

S6 [m3/h]

10

11

15

14

24

25

Figura 6: Mezclador. Al tener el dato del error, podemos detectar los sensores fallados y eliminarlos del cálculo. Todo esto se hace resolviendo el siguiente problema:

M A (F 6

Min

Fj, A i i 1

i

Si)2

i

MA

j

1,2,...,6

(38)

6

Max Fj, A i

i

j

1,2,...,6

(39)

i 1

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

15

F1

F3

F5

(40)

F1

F2

(41)

F3

F4

(42)

F5

F6

(43)

Aj (Fj S j)  1.5

j

1,2,...,6

(44)

A j (F j Sj)  1.5

j

1,2,...,6

(45)

A j  Binario Fj  0

j

j

1,2,...,6

1,2,...,6

(46)

(47)

Note que las restricciones han sido manipuladas para linealizar lo más posible el problema; por ejemplo, se evitó el uso de la función valor absoluto en las ecuaciones que representan el error de los sensores. Esto es conveniente para facilitar la resolución del problema. La variable binaria Aj debe valer 1 si el sensor j funciona correctamente, de lo contrario debe valer 0. De esta forma logramos que los sensores fallados no sean tenidos en cuenta en la estimación de los caudales reales. Sin embargo, cuando un sensor es normal tanto el valor 1 como el 0 satisfacen todas las restricciones; es aquí donde se ve la importancia de la segunda función objetivo que obliga a adoptar el valor 1. Por lo tanto, tenemos un problema con dos objetivos, minimizar el error de la estimación pero utilizando la mayor cantidad de sensores. Para resolver este problema vamos a combinar los objetivos en una única Función Objetivo: A (F M Min MA 6

i

i

Si)2

i 1

Fj, A i

j

1,2,...,6

6

i

(48)

4

i 1

El cuatro del denominador es para decir que aceptamos tan sólo un sensor fallado. La solución Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

16

para este caso es: F1 [m3/h]

F2 [m3/h]

F3 [m3/h]

F4 [m3/h]

F5 [m3/h]

F6 [m3/h]

10.33

10.33

14.33

14.33

24.66

24.66

A1

A2

A3

A4

A5

A6

1

1

1

1

1

1

S1 [m3/h]

S2 [m3/h]

S3 [m3/h]

S4 [m3/h]

S5 [m3/h]

S6 [m3/h]

10

11

15

14

28

25

F1 [m3/h]

F2 [m3/h]

F3 [m3/h]

F4 [m3/h]

F5 [m3/h]

F6 [m3/h]

10.5

10.5

14.5

14.5

25

25

A1

A2

A3

A4

A5

A6

1

1

1

1

0

1

Para el caso siguiente:

La solución es:

Encontramos que un sensor de la salida estaba fallado, el sensor S5, y es fácil ver que las F5 está fuera de la banda normal de este sensor. Sin embargo, si mantenemos la ec. 38 como función objetivo y convertimos la ec. 39 en una restricción de igualdad (igualada a 6), es posible encontrar una solución que mantiene activo el sensor S5. Como vimos anteriormente, este una forma alternativa de manejar un problema con múltiples objetivos. Transformamos el objetivo secundario en una restricción y le exigimos que cumpla con el máximo valor posible que es seis; si el problema no puede resolverse, probamos igualándola a cinco y así sucesivamente. Esto es:

MA 6

i

6

(49)

i 1

La solución encontrada para este caso es: F1 [m3/h]

F2 [m3/h]

F3 [m3/h]

F4 [m3/h]

F5 [m3/h]

F6 [m3/h]

11.25

11.25

15.25

15.25

26.5

26.5

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

17

A1

A2

A3

A4

A5

A6

1

1

1

1

1

1

La explicación a esta aparente contradicción es que si bien era posible encontrar la última solución, el balance expresado por la ec. 48 entre la suma de los errores y la cantidad de lecturas dio preferencia a la minimización de la suma de los errores. Planta de oxígeno Este es un ejemplo de Reklaitis et al. (1983). En la Figura 7 se muestra una planta de producción de acero que utiliza un horno tipo BOF (Basic Oxigen Furnace) operado en modo batch. La demanda de oxígeno es graficada en la Figura 8 este ciclo está determinado por el horno. El oxígeno es suministrado por una planta diseñada para tal fin. La misma está altamente automatizada y produce un caudal constante de oxígeno. Para vincular la planta continua con el horno batch, debemos diseñar un sistema de almacenamiento adecuado. Vamos a considerar que este sistema está formado por una bomba y un tanque. El problema es determinar la producción de la planta de oxígeno y dimensionar el sistema de almacenamiento. Existen varias alternativas, la más simple es fijar la producción de la planta para absorber el pico de demanda y ventilar al ambiente durante el resto del ciclo; sin embargo, esta alternativa tiene un alto costo operativo. En el otro extremo, podemos dimensionar un gran tanque de almacenamiento que permita trabajar a baja producción de oxígeno; pero en este caso el costo de inversión se vuelve inaceptable. Por lo tanto, de este trade-off deducimos que puede existir un óptimo. Otro trade-off se presenta entre la potencia del compresor y las dimensiones del tanque.

Figura 7: Planta de acero.

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

18

Figura 8: Ciclo de producción del horno. Las variables de decisión serán: F: la capacidad de la planta de oxígeno [kg/h]. H: la potencia del compresor [kW]. V: el volumen del tanque [m3]. Para plantear el modelo consideramos primero la ley corregida de los gases: V

Imax R T z M p

(50)

donde R es la constante de los gases, T es la temperatura, M es el peso molecular del oxígeno, p es la presión del tanque, z es el factor de compresibilidad, Imax es la máxima cantidad de oxígeno que puede ser almacenada. De la Figura 8 podemos deducir: Imax

(D1

F) (t2

t1)

(51)

t 1) R T z p

(52)

combinando con la ecuación anterior, tenemos:

V

(D1

F) (t2 M

Asumiendo compresión isotérmica, la potencia del compresor se calcula como:

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

19

H

F) (t2

(D1

t1

t1) R T p ln k1 k2 p0

(53)

donde el primer término es el caudal que bombea el compresor durante el llenado del tanque, k1 es un factor de conversión de unidades, k2 es la eficiencia del compresor, p0 es la presión de la planta productora de oxígeno. Además, la capacidad de la planta de oxígeno debe satisfacer la demanda de todo el sistema. F 

D0 t1

D1 (t2

t 1)

(54)

t2

Más aún: p  p0

(55)

Como función objetivo escogemos el costo anual. Este será la suma de: 1. Costo de operación de la planta de oxígeno (primer término). 2. Costo de capital del tanque (segundo término). 3. Costo de capital del compresor (tercer término). 4. Costo de operación del compresor (cuarto término). C[$/año]

[a1

b

a2 F]

b

d [b3 H 4]

d [b1 V 2]

[N b5 t1 H]

(56)

donde N es el número de ciclos por año y d es el factor de costo anual. En resumen, el problema es: Min C[$/año] F,V,H,p

b

a2 F]

[a1

d [b1 V 2]

b

d [b3 H 4]

[N b5 t1 H]

(57)

sujeto a: V

H

(D1

F) (t2

(D1

M

F) (t2 t1

t 1) R T z p

(58)

t1) R T p ln k1 k2 p0

(59)

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

20

F 

D0 tp1  D p01 (t2 t2

t 1)

(61) (60)

Reactor batch Deseamos determinar la política óptima de suministro de calor a un reactor batch, para minimizar las reacciones colaterales. La formulación esquemática es: tf

Min I Q(t)

Pt

F(c,T,Q,t) dt

(62)

0

sujeto a: dc f(c,T) dt c(t0) c0

(63)

dT g(c,T,Q) dt T(t0) T0

(64)

Q  Qmax

(65)

Q  Qmin

(66)

Tanques buffers Se debe dimensionar los tanques de almacenamiento de una planta química. Los tanques se instalan para compensar las diferencias en las capacidades de producción de los sectores que separan. El tanque i separa el sector i-1con capacidad de producción qi-1 [m3/h] del sector i con una capacidad qi como se muestra en la Figura 9.

Figura 9: Tanques de buffer.

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

21

El modelo matemático es el que se plantea a continuación, donde se consideran N tanques con volumen Vi. Los tiempos de operación de las etapas son ti que no pueden ser menores a un Tmin. Se desprecian los costos de arranque y apagado. Las restricciones dependen de la función que cumple cada tanque, almacenamiento o suministro.

MV N

Min t i, V i, Almacena i

i

i 1

s.a. t i > Tmin Almacenai

~i

(qi 1 > qi)

i > 0

Almacena i  Binario Almacena i (qi 1

[Llenado] Vi [Vaciado] V i

Almacena i q i (t i

qi) ti 1 t i 1)

~i

(1 Almacena i) qi 1 (ti 1 (1 Almacena i) (q i

(67)

t i) i >0

q i 1) t i

i >0

Para entender este modelo conviene analizar los siguientes diagramas de Gantt que presentan un ciclo de trabajo para los dos casos posibles de la variable Almacenai:

Optimización y Simulación de procesos. Optimización Introducción Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

22

Optimización y Simulación de Procesos Enrique Eduardo Tarifa Facultad de Ingeniería - Universidad Nacional de Jujuy

Métodos Analíticos Introducción Hasta aquí hemos analizado la forma de plantear un problema de optimización, ahora es el turno de resolverlo. Si bien existen numerosos software para ese fin (Solver, LINGO, GAMS), poco se podría hacer con ellos si no se tiene una noción por lo menos de la forma en que trabajan. Por ejemplo, suponga que quiere resolver el siguiente simple problema: Max f(x)

x3

9 x2

26 x

24

x

s.a.

x  1.5 x  4.75

(1)

Utilizando el Solver de Excel, la solución reportada (inicializando con x = 1.5) es xopt = 2.42 con yopt = 0.38. Cualquier usuario sin experiencia podría tomar este resultado y utilizarlo ignorando que en realidad no es el óptimo. En efecto, graficando f(x) vemos cuál es el problema:

Figura 1: Extremos relativos y absolutos. El Excel nos da una pista cuando le pedimos el informe de Límites, la solución óptima es xopt = 4.75 con yopt = 3.61, ¡bastante diferente de la reportada al principio!. Este problema aumenta geométricamente con el número de variables. Por otra parte, Solver tiene varios parámetros para regular, sin un conocimiento cabal de sus significados dependeríamos peligrosamente de la capacidad de los programadores para anticipar Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

1

todo tipo de problemas. La Figura 2 muestra el cuadro de Opciones de Solver.

Figura 2: Opciones de Solver (Excel 97). Otro tanto ocurre con LINGO. Este programa nos facilita enormemente la tarea de ingresar el modelo. Por ejemplo, el problema de los sensores planteados en el capítulo anterior se formula de la siguiente manera: MODEL: TITLE Validacion_II; !En este ejemplo se considera la posibilidad de sensores fallados; SETS: !Cada corriente tiene una valor verdadero otro medido y una bandera de estado del sensor. Si la bandera es 1, el sensor funciona bien; Corrientes/Entrada1, Entrada2, Entrada3, Entrada4 Salida1, Salida2/: Caudal, Medido, Activo; ENDSETS DATA: !Lecturas de los sensores; Medido = 10, 11, 15, 14, 24, 25; ENDDATA INIT: Caudal = 10, 11, 15, 14, 24, 25; Activo = 1, 1, 1, 1, 1, 1; ENDINIT !Queremos minimizar el error y maximizar los sensores; [Objetivo] MIN = ((@SUM( Corrientes: Activo*(Caudal-Medido)^2))/ (@SUM( Corrientes: Activo)-4)); !Sujeto a; [Balance_Global] Caudal(1) + Caudal(3) = Caudal(5); [Balance_Local1] Caudal(1) = Caudal(2); [Balance_Local2] Caudal(3) = Caudal(4); Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

2

[Balance_Local3] Caudal(5) = Caudal(6); @FOR( Corrientes: [Positivo] Caudal > 0; [Errores_sup] Activo*(Caudal-Medido) = -1.5; [Binarias] @BIN(Activo)); END

Observe que hay un bloque de inicialización, ¿qué tan sensible es LINGO a los valores especificados en esta sección. Por otra parte, cuál es el significado del reporte de LINGO: Rows= 23 Vars= 12 No. integer vars= 6 Nonlinear rows= 13 Nonlinear vars= 12 Nonlinear constraints= 12 Nonzeros= 63 Constraint nonz= 39 Density=0.211 No. < : 6 No. =: 4 No. > : 12, Obj=MIN Single cols= 0 Local optimal solution found at step: Objective value: Branch count:

22 0.8333333 0

Model Title: VALIDACION_II Variable CAUDAL( ENTRADA1) CAUDAL( ENTRADA2) CAUDAL( ENTRADA3) CAUDAL( ENTRADA4) CAUDAL( SALIDA1) CAUDAL( SALIDA2) MEDIDO( ENTRADA1) MEDIDO( ENTRADA2) MEDIDO( ENTRADA3) MEDIDO( ENTRADA4) MEDIDO( SALIDA1) MEDIDO( SALIDA2) ACTIVO( ENTRADA1) ACTIVO( ENTRADA2) ACTIVO( ENTRADA3) ACTIVO( ENTRADA4) ACTIVO( SALIDA1) ACTIVO( SALIDA2) Row OBJETIVO BALANCE_GLOBAL BALANCE_LOCAL1 BALANCE_LOCAL2 BALANCE_LOCAL3 POSITIVO( ENTRADA1) ERRORES_SUP( ENTRADA1) ERRORES_INF( ENTRADA1) POSITIVO( ENTRADA2) ERRORES_SUP( ENTRADA2) ERRORES_INF( ENTRADA2) POSITIVO( ENTRADA3) ERRORES_SUP( ENTRADA3) ERRORES_INF( ENTRADA3) POSITIVO( ENTRADA4) ERRORES_SUP( ENTRADA4) ERRORES_INF( ENTRADA4) POSITIVO( SALIDA1)

Value 10.33333 10.33333 14.33333 14.33333 24.66667 24.66667 10.00000 11.00000 15.00000 14.00000 24.00000 25.00000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000

Reduced Cost 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -0.3611107 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

Slack or Surplus 0.8333333 -0.3552714E-14 0.0000000 0.0000000 0.0000000 10.33333 1.166667 1.833333 10.33333 2.166667 0.8333328 14.33333 2.166667 0.8333326 14.33333 1.166667 1.833333 24.66667

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

Dual Price -0.3333333 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

3

ERRORES_SUP( ERRORES_INF( POSITIVO( ERRORES_SUP( ERRORES_INF(

SALIDA1) SALIDA1) SALIDA2) SALIDA2) SALIDA2)

0.8333346 2.166665 24.66667 1.833335 1.166665

0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

Para lograr un uso correcto de estas potentes herramientas, es necesario dominar algunos conceptos teóricos que serán tratados en este capítulo.

Funciones de una sola variable Introducción La optimización de funciones con una única variable independiente es frecuente en el campo de la ingeniería. Además, varios métodos diseñados para múltiples variables utilizan optimizaciones internas univariables. Por este motivo, son muchos los métodos desarrollados para las funciones de una única variable. Antes de resolver cualquier problema, debemos estudiar las funciones a fin de saber qué esperar de sus comportamientos. Funciones continua Una función f(x) es continua en a si: 1. Existe f(a). 2. Existe el Lim f(x) para x tendiendo a a. 3. El valor de (1) es igual al valor de (2). Función monotónica Una función f(x) es monotónica si para todo par de puntos x1, x2 con x1  x2, se cumple f(x1)  f(x2) (monótona creciente), o se cumple f(x1)  f(x2) (monótona decreciente). No es necesario que se trate de una función continua. Función unimodal Un función es unimodal en a si es monótona creciente para x  a y monótona decreciente para x  a. También vale para la situación inversa. Definición Una función es unimodal en el intervalo [a, b] si y sólo si es monotónica en ambos lados del único punto óptimo xopt en el intervalo. En otras palabras, si xopt es un punto mínimo de f(x) en dicho intervalo, luego la función es unimodal si y sólo si para todo punto x1 y x2, x opt  x1  x2 < f(x opt)  f(x1)  f(x2)

(2)

x opt  x1  x2 < f(x opt)  f(x1)  f(x2)

(3)

y

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

4

Criterio de optimalidad Cuando se está resolviendo un problema de optimización, aparecen dos situaciones: Situación estática: ¿Cómo podemos determinar si un dado punto x* es el la solución • óptima? Situación dinámica: Si el punto x* no es el óptimo, entonces ¿cómo llegamos a la • solución óptima? En este capítulo trataremos de contestar la primera pregunta. Optimo global y local Sea f(x) definida en la región R, el punto xopt  R es un mínimo global si y sólo si f(xopt)  f(x) para todo x  R. Sea f(x) definida en la región R, el punto xopt  R es un mínimo local si y sólo si existe un valor / >0 tal que f(xopt)  f(x) para todo x que cumpla con |x - xopt| < /. Observaciones: • Invirtiendo las desigualdades tenemos las definiciones para los máximos. • Si la función es unimodal, el mínimo local es automáticamente el mínimo global. • Cuando la función no es unimodal, pueden existir múltiples puntos mínimos locales, y el mínimo global sólo puede ser determinado buscando todos los locales y seleccionando el mejor. Funciones cóncavas y convexas Como demostraron Wilde y Beightler (1967) la propiedad de unimodalidad es difícil de establecer analíticamente. Entonces, debemos contar con algún criterio más fácil de establecer que nos permita determinar si el problema tendrá un único óptimo o no. Este criterio es el que analizaremos ahora. Una función es cóncava sobre una región R si la siguiente relación vale para todo par de puntos xa y xb (pueden ser vectores) pertenecientes a R: f[ xa

(1

) x b]   f(xa)   [0,1]

(1

) f(x b)

(4)

Si se prohíbe la igualdad, entonces tenemos una función estrictamente cóncava. Cambiando el sentido de la desigualdad podemos definir las funciones convexas. En la Figura 3 se muestran algunos ejemplos. Observe, que la función todas son también unimodales, sin embargo la última no es ni cóncava ni convexa. Otra observación interesante es que una lía recta es cóncava y también convexa, pero no satisface las condiciones de estrictamente cóncava ni convexa. Además, la suma de dos funciones mantiene el carácter. Las ecuaciones utilizadas para definir las funciones cóncavas y convexas no son muy cómodas Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

5

para calcular. En cambio, para este fin, se utiliza la segunda derivada. Una función es cóncava sobre una región R si f”(x)  0 para todo x  R. Esta definición se puede adaptar para las funciones convexas y las estrictamente cóncavas y convexas con sólo modificar la desigualdad.

Figura 3: Funciones unimodales. Región convexa Una región R es convexa si para todo par de puntos xa y xb (pueden ser vectores) pertenecientes a R se cumple: x

 x a (1 ) x b   [0,1] x  R

(5)

En la Figura 4 se muestran algunos ejemplos para dos dimensiones.

Figura 4: Regiones convexas y no convexas. Es muy importante contar con un método para determinar si una región es o no es convexa porque esto afecta la solución del problema; en la Figura 5 se muestran dos casos donde una función cóncava está siendo maximizada. En el primer caso, el dominio es convexo, por lo tanto basta con encontrar el único máximo relativo (fácil de encontrar) que automáticamente será también el máximo absoluto. En cambio, en el segundo caso misma función es maximizada en un dominio no convexo; ahora, es necesario estudiar lo que ocurre en las fronteras. Se puede determinar si una región es convexa o no estudiando las restricciones que la determinan. Si todas las restricciones g(x) son funciones cóncavas y además g(x)  0; entonces, tenemos una región cerrada convexa. Lo mismo vale si todas las restricciones con convexas pero ahora g(x)  0. Recuerde que las líneas rectas son tanto cóncavas como convexas. En la Figura 6 se muestran los ejemplos para una sola variable, note que lo que se está graficando ahora son Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

6

las restricciones, no la función objetivo. En la Figura 7 se combinan los gráficos, f(x) es la función objetivo y g(x) es la restricción.

Figura 5: Efecto de la región sobre el óptimo.

Figura 6: Efecto de la restricción sobre la región.

Figura 7: Efecto de la región sobre el óptimo. Condiciones para un punto óptimo Los valores óptimos de una función limitada a una región sólo pueden darse en limitadas Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

7

localizaciones. Por eso, es importante identificar estas posiciones para limitar la búsqueda a ellas. Las localizaciones posibles de un óptimo son: 1. Los puntos estacionarios, donde f’(x) = 0. 2. En los puntos de discontinuidad de la f’(x). 3. En los puntos de discontinuidad de la f(x). 4. En las fronteras. La Figura 8 muestra los puntos que deben ser explorados en la búsqueda de un máximo en el intervalo [a,b].

Figura 8: Puntos extremos. Condiciones para un óptimo interno o sin restricciones Un óptimo interno a una región, o un óptimo de un problema sin restricciones puede estar dado por un punto estacionario. Aquí veremos la condiciones necesarias y suficientes para que un dado punto sea considerado un punto óptimo. Si f(x) es continua y sus derivadas también, entonces podemos de desarrollar la serie de Taylor alrededor del punto candidato a: f(a h)

f(a)

h f (a)

h2 f (a) 2!

...

h j (j) f () j!

(6)

donde h = x - a, y   [a, x]. La condición necesaria para que en a exista un punto óptimo se deriva de la siguiente manera. Supongamos que en x = a existe un mínimo, debido a que para valores pequeños de h es el primer término de la mano derecha el que domina, tenemos que el signo de la diferencia f(a+h) f(a) depende del signo de h, hecho que se opone con la definición de mínimo. Para asegurar que la definición de mínimo no es violada no queda otra alternativa que f’(a) = 0. Esta es la condición necesaria para que un punto interno o sin restricciones sea mínimo. De la misma manera se puede demostrar para un punto máximo. Sabiendo que el primer término debe ser nulo, ahora nos concentramos en el segundo término. Para h pequeños, este término predomina sobre el resto y fija el signo de la diferencia f(a+h) Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

8

f(a). Se puede ver que si f”(a) > 0, nos aseguramos que se cumpla la definición de mínimo independientemente del signo de h. Esta es una condición suficiente para la existencia de un mínimo. Para un máximo pediremos f”(a) < 0. En forma más general, sea n es el orden de la primera derivada que no es nula en el punto x = a, entonces: • Si n es par y f(n)(a) > 0, entonces tenemos un mínimo. • Si n es par y f(n)(a) < 0, entonces tenemos un máximo. • Si n es impar, entonces tenemos un punto de inflexión o ensilladura.

Funciones de varias variables Introducción Las mismas situaciones consideradas para las funciones de un sola variables se presentan también aquí, pero con el agravante que los inconvenientes crecen geométricamente con el número de variables a considerar. La extensión de las mayorías de las definiciones antes realizadas son triviales, basta con pensar que la variable x ahora es un vector. Aquí sólo se repasarán las definiciones que incorporen algún grado de dificultad. En primer lugar, la f’(x) será reemplazada por el gradiente /f(x), mientras las segunda derivada será reemplazada por /2f(x) o Hessiano H(x) definido como:

H(x)

02f(x) 02f(x) 0x12 0x10x2 02f(x) 02f(x) 0x20x1 0x22

(7)

Funciones cóncavas y convexas Primero debemos plantear las siguientes definiciones: • H(x) es definida positiva si y sólo si vT H v > 0 para todo v g 0. • H(x) es definida negativa si y sólo si vT H v < 0 para todo v g 0. • H(x) es indefinida si y sólo si vT H v es < 0 para algunos v y > 0 para otros. • Cuando se contempla la igualdad, H es semidefinida. Se puede demostrar, utilizando la serie de Taylor que: • Si H(x) es definida negativa en la región R, la función f(x) es estrictamente cóncava en esa región. • Si H(x) es definida positiva en la región R, la función f(x) es estrictamente convexa en esa región. • Si H(x) es semidefinida, entonces f(x) será sólo cóncava o convexa. Para clasificar H(x) como definida positiva, se debe cumplir: • Todos los elementos de la diagonal deben ser positivos y los determinantes de H(x) y de Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

9

•

sus principales menores deben ser positivos. O, todos los autovalores de H(x) deben ser positivos.

Para clasificar H(x) como definida negativa, se debe cumplir: • Todos los elementos de la diagonal deben ser negativos y los determinantes de H(x) y de sus principales menores deben tener signos alternos comenzando con det(M1(H(x))) < 0, donde M1 es la matriz principal 1 (el elemento de la esquina superior izquierda). • O, todos los autovalores de H(x) deben ser negativos. Para los H semidefinidos, basta incorporara la igualdad. Cuando los autovalores tienen signos mezclados, la función no es ni cóncava ni convexa. Los autovalores () y autovectores (v) se definen de la siguiente manera:  v

H v

(8)

Los autovalores se obtienen de resolver la siguiente ecuación: |H

 I|

0

(9)

Mientras que los autovectores se obtienen de resolver el siguiente sistema de ecuaciones, una vez por cada autovalor (la solución no es única): i I) v i

(H

0

(10)

Ejemplo Clasifique la siguiente función: f(x)

2

2 x1

3 x1 x2

2

2 x2

(11)

Calculando el Hessiano: 0f(x) 0x 1

4 x1

3 x2

(12)

0f(x) 0x 2

3 x1

4 x2

(13)

02f(x) 0x12

4

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(14)

10

02f(x)

4

0x22

02f(x) 0x10x

(15)

02f(x) 0x20x

2

3

(16)

1

Entonces, el Hessiano es: 4

H(x)

3

(17)

3 4

Los menores principales son: M1 (orden 1) = 4 det(M1) = 4 M2 (orden 2) = H det(M2) = 7 por lo tanto, H(x) es definida positiva y f(x) es estrictamente convexa. Ejemplo Clasifique la siguiente función utilizando autovalores: f(x)

2

2 x1

2

2 x1 x2

1.5 x2

7 x1

8 x2

24

(18)

El Hessiano es: 4 2

H(x)

(19)

2 3

Para determinar los autovalores planteamos: |H(x)

 I|



4 2

2 3



(4

) (3

)

4

0

(20)

los autovalores son 5.56 y 1.44, ambos positivos; por lo tanto, la función es estrictamente convexa. Regiones convexas Se mantiene la misma definición pero ahora se utiliza el Hessiano. Veamos un ejemplo.

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

11

Ejemplo Determine si las siguientes restricciones forman una región convexa: x2  1

(21)

x2  2

(22)

2

x1

x1

Primero escribimos las restricciones en la forma g(x) 0: g1(x)

x1

2

x2

1  0

(23)

g2(x)

x1

x2

2  0

(24)

Debemos demostrar que las dos funciones son cóncavas, entonces calculamos los Hessianos: H[g1(x)]

H[g2(x)]

2 0 0 0

0 0 0 0

(25)

(26)

como todos los autovalores son nulos o negativos, las matrices son semidefinidas negativas; por lo tanto, las funciones son cóncavas, y la región es convexa. La Figura 9 muestra esta región.

Figura 9: Región convexa. En la Figura 10 se puede ver los problemas que trae una región no convexa, se pueden Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

12

determinar dos puntos óptimos.

Figura 10: Región no convexa. Condiciones necesarias y suficientes para un óptimo sin restricciones La demostración es parecida a la anterior, sólo que ahora se utiliza el gradiente y el Hessiano. La serie de Taylor alrededor del punto xopt es: f(x)

1 (ûX)T /2f(x opt) 2

/Tf(x opt) ûx

f(x opt)

ûx

O3(ûx)

...

(27)

Siguiendo el mismo razonamiento que antes se tienen las condiciones necesarias (1 y 2) y la suficiente (de 1 a 3): f(x) es dos veces diferenciables en xopt. 1. /f(xopt) = 0, esto es xopt es un punto estacionario. 2. 3. H(xopt) es definida positiva para que un mínimo exista en xopt, y definida negativa para que exista un máximo en xopt.

Ejemplo Calcule los puntos extremos de la siguiente función: f(x)

4

4.5 x1

4 x2

2

x1

2

2 x2

2 x1 x2

4

2

2 x1 x2

x1

(28)

Haciendo el gradiente cero, tenemos el siguiente sistema de ecuaciones:

0f(x) 0x 1

4.5 2 x1

2 x2

3

4 x1

4 x1 x2

0

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(29)

13

0f(x) 0x 2

4

4 x2

2

2 x1

2 x1

0

(30)

La solución numérica de este sistema y su clasificación se da en la siguiente tabla: Puntos estacionarios (x1, x2)

f(x)

1

2

Clasificación

(1.941, 3.854)

0.9855

37.03

0.97

Mínimo local

(-1.053,1.028)

-0.5134

10.5

3.5

Mínimo local

(0.6117,1.4929 )

2.83

7.0

-2.56

Ensilladura

Optimo con restricciones Ahora analicemos cómo afectan las restricciones al problema de encontrar un punto óptimo. Consideremos el problema general: x s.a. h j(x)

Min f(x) [ x1, x2, ..., x n]T j

0

gj(x)  0

j

m

1, 2, ..., m 1, m

(31)

2, ..., p

Las desigualdades (restricciones débiles) se pueden transformar en igualdades por un método que veremos luego, por lo tanto el problema se reduce a manejar las igualdades (restricciones fuertes). La forma más aconsejable de manejar las restricciones de igualdades es utilizar las mismas para eliminar una cantidad equivalente de variables (método de sustitución) y transformar el problema original en un problema sin restricciones y con menor cantidad de variables. Este método ya fue explicado en un capítulo anterior. A continuación veremos un ejemplo: Ejemplo Min f(x) s.a. 2 x1

2

4 x1

2

5 x2

(32) 3 x2

6

Utilizamos la igualdad para despejar x1:

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

14

x1

3 x2

6

(33)

2

sustituyendo en la función objetivo, tenemos la nueva función objetivo de una única variable y sin restricciones: Min f(x2)

2

14 x2

36 x2

36

(34)

Como es una función de una única variable podemos resolverlo con el método estudiado en la sección anterior: 0f(x2) 0x 2

28 x2

36

0

(35)

entonces, x2opt = 1.286, y por lo tanto x1opt = 1.071. Multiplicadores de Lagrange Para aplicar el método de los multiplicadores de Lagrange necesitamos llevar el problema original a una forma donde sólo existan restricciones fuertes (igualdades). A fin de analizar el fundamento del método, veamos el siguiente problema: Min f(x1, x2) s.a. h(x1, x2) 0

(36)

para que un punto sea considerado óptimo deberá satisfacer las siguientes ecuaciones: df

0

dh

0

0f dx1 0x 1 0h dx1 0x 1

0f dx2 0x 2 0h dx2 0x 2

(37)

Debido a que las variables x1 y x2 no son independientes (están vinculadas por la restricción), no es necesario pedir que ambas derivadas parciales de f sean nulas, sino que simplemente pediremos que ambos término se cancelen mutuamente asegurando un df nulo (condición para ser un óptimo). Por otra parte, la segunda igualdad se debe cumplir porque h está igualada a una constante (cero). Las ecs. 37 constituyen un sistema de ecuaciones homogéneo, el cual tendrá solución no trivial (dx1 y dx2 no nulos) cuando el determinante sea igual a cero, esto es: Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

15

0f 0h 0x 1 0x 2

0f 0h 0x 2 0x 1

0

(38)

esto se puede escribir también de la siguiente forma: 0f 0x 2

0h 0x 2

0f 0x 1

0h 0x 1

0f 0x 2

0f 0x 1

0h 0x 2

0h 0x 1

(39)

o también:

(40)

o lo que es equivalente: 0f 0x 1 0f 0x 2

0h 0x 1 0h  0x 2 

0 (41) 0

donde  relaciona las derivadas de f y h y es llamado multiplicador de Lagrange. Por comodidad, podemos plantear la siguiente función ampliada (Lagrangeana): L(x, )

f(x)

 h(x)

(42)

note que h(x) = 0, Si aplicamos a esta función las condiciones para un óptimo sin restricciones, obtendremos el siguiente sistema de ecuaciones:

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

16

0L 0x 1 0L 0x 2

0f  0x 1 0f  0x 2 0L h(x) 0

0h 0x 1 0h 0x 2

0 0

(43)

0

con incógnitas , x1 y x2. Observe que estas ecuaciones son las que obtuvimos en la ec. 41 y en la ec. 36 como las condiciones que debía satisfacer un punto óptimo. En la Figura 11 se muestra la interpretación gráfica de la condición expresada en la ec. 39. La conclusión es que en el punto óptimo, la restricción se hace tangente a las curvas de nivel de la función objetivo. Recuerde que el gradiente de una función es un vector perpendicular a la curva de nivel e indica hacia el máximo de la función.

Figura 11: Multiplicadores de Lagrange. Ejemplo 2

2

Min f(x) 4 x1 5 x2 s.a. h(x) 2 x1 3 x2 6 0

(44)

La función ampliada será entonces: L(x,)

2

4 x1

2

5 x2

 (2 x1

3 x2

6)

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(45)

17

aplicando las condiciones necesarias: 0L 8 x1 2  0 0x 1 0L 10 x2 3  0 0x 2 0L 2 x1 3 x2 6 0 0

(46)

la solución es opt = - 4.286, x1opt = 1.071 y x2opt = 1.286. En este caso la solución fue encontrada con tan sólo resolver un sistema lineal de ecuaciones. Sin embargo, en los casos en que el sistema a resolver es no lineal deberemos sopesar la conveniencia de utilizar el método de los multiplicadores debido a que la resolución de un sistema de ecuaciones lineales es tan o más complejo que el problema original de optimización. Se debe siempre verificar el xopt debido a que la función Lagraneana exhibe un punto de ensilladura con respecto a x y . Generalización del método de los multiplicadores de Lagrange A continuación generalizaremos el método para más de dos variables e incluiremos restricciones débiles. Dado el modelo general planteado anteriormente: Min f(x) [ x1, x2, ..., x n]T

x s.a. h j(x)

j

0

gj(x)  0

j

(47)

1, 2, ..., m

m

1, m

2, ..., p

Lo rescribimos para adaptarlo al método de los multiplicadores de Lagrange: Min f(x) [ x1, x2, ..., x n]T

x s.a. h j(x) g j(x)

j

0

12j

j

0

m

(48)

1, 2, ..., m 1, m

2, ..., p

Note que al escribir 12 evitamos incorporar las restricciones 1  0. La Lagrangeana es:

M m

L(x, )

f(x)

j 1

M  (g (x) p

j h j(x)

j

j

12j )

(49)

j m 1

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

18

Las condiciones necesarias y suficientes que debe cumplir xopt son: • f(x) debe ser convexa. • La región debe ser convexas. • Satisfacer las siguientes ecuaciones: 0L 0x i

0

i

1, ..., n

0L 0j

0

j

1, ..., p

0L 01j

(50)

2 j 1j j

0

0

j

m

1, ..., p

j

m

1, ..., p

Para un máximo de f(x) se debe cambiar el signo de la última desigualdad por j 0. Observe que el valor de j se elige de tal modo que penaliza el no cumplimiento de las restricción de desigualdad. Ejemplo Resolvamos el siguiente problema: Min f(x) s.a. g(x)

25

x1 x2

2 x1

2 x2

 (25

x1

(51)

0

La Lagrangeana es: L(x,)

x1 x2

2

2

x2

1 2)

(52)

Planteando las condiciones necesarias:

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

19

0L 0

0L 0x 1 0L 0x 2

x2

2  x1

0

x1

2  x2

0

2

x1

25 0L 01

2

x2

2  1 

12

(53)

0

0

0

El sistema es homogéneo, por lo tanto las ecuaciones son dependiente, se debe probar para  = 0 y para  no nulo. La solución para estos casos es: 

x1

x2

1

f(x)

Comentarios

0

0

0

5

0

Ensilladura

-0.5

3.54

-3.54

0

-12.5

Mínimo

-0.5

-3.54

3.54

0

-12.5

Mínimo

0.5

3.54

3.54

0

12.5

Máximo

0.5

-3.54

-3.54

0

12.5

Máximo

Interpretación económica de los multiplicadores Para una dada restricción, el correspondiente multiplicador indica cuánto cambiará la función objetivo cuando se cambie diferencialmente la constante e de la restricción h*(x) = e. Esto es:

0L 0e

0 [f(x) 0e

 (h (x)

e)]



(54)

y debido a que en óptimo h*(xopt) = e y xopt es función de e, tenemos que:

0L 0e

0f 0e

(55)

Entonces, el multiplicador indica la sensibilidad de tanto L(x) como de f(x) con respecto a e. El parámetro e puede ser interpretado como los recursos existentes, mientras el correspondiente multiplicador será el precio sombra. El estudio de estos multiplicadores nos permite decidir dónde conviene colocar inversiones adicionales, se colocará en las restricciones que afecten más favorablemente la función objetivo. Condiciones necesarias y suficiente para un mínimo local, Kuhn-Tucker La Lagrangeana considerando las desigualdades es: Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

20

M

M

m

L(x, , µ)

f(x)

j

p

j h j(x) j

1

µ j (gj(x)

12j )

(56)

m 1

Los  pueden ser positivos o negativos, pero los µ deben ser positivos. Si la desigualdad de las restricciones fuera del tipo menor o igual a cero, debería ir sumando en lugar de ir restando el la sumatoria correspondiente a las restricciones débiles. Entonces, siempre pediremos µ 0. Una dirección válida de búsqueda s, debe asegurarnos que disminuirá el valor de f(x). Esto se verifica si se cumple /Tf(x) s < 0. Recuerde que el gradiente de una función indica hacia el máximo, por lo tanto una dirección cuya proyección sobre el gradiente es negativa debe necesariamente conducir hacia un mínimo. Analizando la Figura 12, vemos que el único modo de que no exista una dirección s que posibilite continuar minimizando (y por lo tanto ya estaríamos en un mínimo) es que el gradiente sea una combinación lineal de los gradientes de las restricciones activas, es decir:

/f(x opt)

µ 1 [ /g1(x opt)] opt

µ 2 [ /g2(x opt)] opt

(57)

Otra forma de decir esto es, la intersección del cono de direcciones descendentes (toda s con proyección negativa sobre /f(x)) con el cono de direcciones factibles (todas las s con proyección positivas sobre los /g(x) de las restricciones activas) es vacía. Generalinzando, podemos afirmar que debe cumplirse lo siguiente en un punto óptimo:

M p

/f(x opt) j

opt

µj

[ /g j(x opt)]

m 1 opt µj

0

(58)

Sin embargo, como la relación debe cumplirse sólo para las restricciones activas, imponemos: opt

0

opt

0

µj

µj

si gj(x opt)

0

si gj(x opt) >0

(59)

o lo que es lo mismo: opt

µj

g j(x opt)

0

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(60)

21

Figura 12: Condición de dependencia lineal. Por otra parte, de la definición de la Lagrangeana, se cumple:

/L(x opt, µ opt)

/f(x opt)

M m

j

1

opt

j

M p

/hj(x opt) j

opt

µj

/g j(x opt)

(61)

m 1

Note que no figuran las variables 1 porque se hacen nulas para las restricciones activas, mientras que para las restricciones no activas µ es nulo. Dada la relación deducida anteriormente entre los gradientes de las restricciones y de la función, concluimos que /L(xopt) debe ser nulo. En resumen, las condiciones necesarias para que xopt sea un mínimo son: f(x), hj(x) y gj(x) son todas diferenciables dos veces en xopt. 1. 2. Las restricciones son satisfechas (tanto las fuertes como las débiles). 3. Vale la “calificación de restricciones de segundo orden”. Para esto, es suficiente con que los gradientes de las restricciones activas (/gj(xopt)) y de las igualdades (/hj(xopt)) sean linealmente independientes. 4. Vale µ jopt gj(xopt) = 0. 5. Vale /Lx(xopt, opt, µ opt) = 0. 6. Los µ jopt son no negativos. 7. La Hessiana de L sea semidefinida positiva para todo v para el cual vT*/gj(xopt) = 0 (de las activas) y vT*/hj(xopt) = 0: vT*/2[L(xopt, opt, µ opt)]*v  0 Las condiciones suficientes para que xopt sea un mínimo son: 1. De la 1 a la 6 de las necesarias. 2. La Hessiana de L es positiva definida para todo v para el cual: T v */gj(xopt) = 0, para las activas con µ jopt > 0. vT*/hj(xopt) = 0 Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

22

vT*/gj(xopt)  0, para las activas con µ jopt = 0. vT*/2[L(xopt, opt, µ opt)]*v > 0 Observaciones: si f(x) es convexa, las hj(x) son lineales, las gj(x) son cóncavas y xi 0, entonces un mínimo local también es un mínimo global. Análogamente, un máximo local es un máximo global si la función objetivo es cóncava y las restricciones forman un conjunto convexo. Note que todas las condiciones que vimos pueden ser utilizadas para: 1. Verificar si un candidato dado es un punto óptimo. 2. Para encontrar candidatos. Si se hace esto, hemos transformado el problema de Optimización en un problema de resolución de un sistema de ecuaciones. Generalmente, este nuevo problema puede ser tanto o más complejo que el problema inicial, y por lo tanto conviene la utilización de otro método. Ejemplo Determinar si el punto (1.00, 4.90) es un mínimo local para el problema. Min f(x)

4 x1

2

x2

12

s.a. g1(x)

h1(x)

25

10 x1

2 x1

g2(x)

2

x1

10 x2

2

x2

0 2 x2 2

34  0

(62)

(x1 3) (x2 1)  0 g3(x) x1  0 g4(x) x2  0 2

En la Figura 13 se muestra gráficamente este problema. 1. 2. 3.

Por inspección, las funciones son dos veces diferenciables. Reemplazando el valor de xopt en las restricciones vemos que todas las restricciones se satisfacen y que g1(xopt) se vuelve activa (hace valer la igualdad). Ahora verificamos la independencia lineal entre los gradientes de las restricciones activas, /g1(xopt) y /h1(xopt). Para ello, analizamos el determinante del sistema de ecuaciones que surge de plantear la combinación lineal. Como el determinante del sistema es distinto de cero, la única solución es el vector nulo (c = (0,0)), por lo tanto los gradientes son linealmente independientes.

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

23

Figura 13: Condiciones de Kuhn-Tucker.

opt

c2 ( 2 x1 )

opt

c2 ( 2 x2 )

c1 (10

2 x1 )

c1 (10

2 x2 )

opt

0

opt

0

(63)

es decir: 2 c2 0 9.8 c2 0

8 c1 0.2 c1

4. 5.

(64)

µ 1 puede tener cualquier valor porque g1(xopt) = 0, pero µ 2, µ 3 y µ 4 deberán ser nulos. Debemos verificar que el gradiente de la Lagrangeana es nulo. opt

1

4 2

opt x2

opt

opt

( 2 x1 )

opt 1

(2

u1

opt x2 )

opt

2 x1 )

(10

opt u1

(10

2

opt

0

opt x2 )

0

0 0

0

(65)

es decir: 4 9.8

opt

2 1 9.8

opt 1

8 u1

0.2

opt u1

0

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(66)

24

La solución de este sistema es 1opt = -1.016 y µ 1opt = 0.754. 6. 7.

Estos son valores permitidos por las condiciones que estamos verificando ya que µ 1opt debe ser positivo, mientras 1opt puede tener cualquier valor. Calculamos la Hessiana de la Lagrangeana: 02L

2 1

2

0x 1 02L 2

0x 2

2 µ1 (67)

2

2 1

2 µ1

Las derivadas cruzadas son nulas, por lo tanto el Hessiano en xopt es: H(x opt)

3.54

0

0

1.54

(68)

es definida positiva. Sin embargo, el único vector v que satisface vT*/g1(xopt) = 0 y vT*/h1(xopt) = 0, es el vector nulo; por lo tanto, sólo podemos demostrar la condición necesaria vT*/2[L(xopt, opt, µ opt)]*v  0, pero no la suficiente vT*/2[L(xopt, opt, µ opt)]*v > 0. Lo que ocurrió es que para este problema, en xopt no tenemos grados de libertad porque tenemos dos variables y dos restricciones fuertes.

Optimización y Simulación de procesos. Métodos Analíticos Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

25

Optimización y Simulación de Procesos Enrique Eduardo Tarifa Facultad de Ingeniería - Universidad Nacional de Jujuy

Modelos y Función Objetivo Introducción Este capítulo está basado en los capítulos II y III del libro “Optimization of Chemical Processes” (Edgar y Himmelblau, 1988), debe completarse con la lectura del capítulo II y III del libro “Optimization: Theory and Practice” (Beveridge y Schechter, 1970). En este capítulo nos referiremos a la formulación de las restricciones y de la función objetivo necesarias para plantear un problema de optimización. Ambos temas son demasiados extensos, cada uno de ello justifica un curso aparte. Sin embargo, dado que necesitamos tener algunas nociones de estos temas para poder plantear problemas de optimización, aquí se abordará someramente cada uno de ellos.

Restricciones y modelos de un sistema Las restricciones que figuran en la formulación general de un problema de optimización, por lo general, surgen del modelo del sistema que está siendo optimizado. Cuando vamos a desarrollar un modelo debemos considerar los siguientes aspectos: 1. ¿Cuál es el nivel de abstracción más adecuado, microscópico o macroscópico; y cuál es el grado de detalle requerido? 2. ¿Se puede describir el sistema utilizando principios de físico-química? 3. ¿Cuál es la exactitud requerida y cómo afecta ésta a la aplicación final? 4. ¿Cuáles son las variables medidas y cuáles los datos disponibles para la validación del modelo? 5. ¿Se puede subdividir el proceso en subprocesos más simples? Si el modelo es muy simple, los resultados no serán confiables; pero si el modelo es demasiado complicado, la optimización será también muy difícil sino imposible. Por lo tanto, existe un modelo óptimo que podemos definirlo como aquel que es el más simple pero que satisface nuestras necesidades.

Clasificación de Modelos Podemos clasificar los modelos en las siguientes categorías: 1. Basados en los primeros principios vs. Empíricos. 2. Lineales vs. No lineales. 3. Estacionarios vs. Dinámicos. 4. Parámetros concentrados vs. Parámetros distribuidos. 5. Con variables continuas vs. Discretas. 6. Determinísticos vs. Estocásticos. Un ejemplo del primer tipo es el siguiente modelo para un reactor tubular isotérmico con Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

1

dispersión axial: 0c 0t

Dz

0 2c 0z

v

2

0c 0z

r

(1)

No existe único modelo, una alternativa posible es ver el reactor como un conjunto de p reactores agitados continuos en serie como muestra la Figura 1.

Figura 1: Reactor tubular. Como modelo empírico podemos citar los siguientes modelos propuestos para un precipitador electrostático: 



b1 A

100 1

1

b2

e

(2)

2 A

1 3 A

(3)

donde  es la eficiencia y A es la superficie de recolección. El principal inconveniente de los modelos empíricos es que están restringido al dominio que contiene los datos utilizados para ajustar el modelo. Además, el sistema a modelar debe existir necesariamente para poder obtener los datos requeridos. Todo esto está atenuado en los modelos con base teórica, sin embargo, el costo de desarrollo puede ser grande. Los modelos lineales tienen la ventaja del principio de la superposición; esto es, dado un modelo lineal J, se cumple: J(a x1

b x2)

a J(x1)

b J(x2)

(4)

El primer modelo utilizado para el precipitador es lineal, el segundo no lo es. Existe una teoría matemática ampliamente desarrollada para tratar modelos lineales. No ocurre lo mismo con los modelos no lineales. Por lo tanto, siempre que sea posible, se prefiere un modelo del tipo lineal. En los modelos estacionario no se considera la evolución de las variables con el tiempo; es decir, que las derivadas parciales con respecto al tiempo son nulas. Este tipo de modelos es ampliamente utilizado para optimizar. Cuando se utilizan modelos dinámicos, entramos en el Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

2

campo del “control óptimo”. Cuando se desprecia la variación de las propiedades con respecto a las coordenadas espaciales, se dice que el modelo es del tipo parámetro concentrado; por ejemplo, un tanque perfectamente agitado. Si las propiedades del sistema están en función de algunas de las coordenadas espaciales, tenemos un modelo con parámetro distribuido; por ejemplo, un reactor tubular. Cuando la respuesta de un modelo está perfectamente determinada por las entradas y los parámetros estamos frente un modelo determinístico (por ejemplo, cualquiera de los modelos anteriores). Cuando a pesar de conocer las entradas y los parámetros existe incertidumbre en la respuesta se está frente un modelo estocástico (por ejemplo, el modelo estadístico utilizado para predecir la población de un país).

Desarrollo de un modelo En la Figura 2 se esquematizan los pasos involucrados en el desarrollo de un modelo.

Experiencia, realidad

Objetivo, planificación

Objetivos de la empresa

Definición del problema

Determinación de variables y primeros principios

Simulación

Desarrollo del modelo

Observación, datos

Diseño Estimación de parámetros

Evaluación y verificación

Evaluación

Aplicaciónd del modelo

Figura 2: Desarrollo de un modelo.

Estimación de parámetros, método de los mínimos cuadrados Uno de los métodos más utilizados para ajustar los parámetros de un modelo es el método de los mínimos cuadrados. Se trata de un problema optimización donde se plantea un modelo lineal en los parámetros  de la forma: Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

3

M n

y

i x i

(5)

i 0

x0

1

por otra parte tenemos los vectores de datos X e Y con p puntos experimentales. El problema de ajustar los coeficiente se formula como sigue:

M

M

p

Min 

n

Yj

j 1

2

i Xji

(6)

i 0

X j0

0

Este problema puede ser resulto analíticamente. Para ello, derivando la ecuación anterior con respecto a los  e igualando a cero (aplicando el operador 0/0i = 0), tenemos:

0 0 0

M M MX X  M MX X  M p

2

1

Xj0

j 1 p

p

Xj0Xj1

j 1

p

j1 j0

1

j 1 p

j2 j0

1

2

M M p

j 1

Xj2Xj1

MX X

jn j0

1

M

Xj0Xj2

j 1 p

2

j 1

p

j 1

2

Xj1

j 1 p

j 1

0

2

Xj1Xj2

M p

2

Xj2

... n ... n ... n

j 1

M M M p

j 1 p

p

XjnXj1

j 1

2

M

Xj1Xjn

XjnXj2

... n

j 1

YjXj1

j 1 p

Xj2Xjn

j 1

p

YjXj0

j 1 p

j 1 p

...

M M M p

Xj0Xjn

YjXj2

(7)

j 1

M p

2 Xjn

j 1

M p

YjXjn

j 1

1 X11 X12 ... X1n X

1 X21 X22 ... X2n ... ...

... ... ...

(8)

1 X p1 Xp2 ... Xpn

Usando notación matricial se puede escribir en una forma más compacta: Entonces: XT X 

XT Y

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

(9)

4

0 1 

2

(10)

...

n

despejando (note que la matriz X no es cuadra y por lo tanto no tiene inversa):

 (X T X) 1 X T Y

(11)

Y1 Y

Y2

(12)

... Yp

Ejemplo. Costos de fabricación de intercambiadores Se necesita estimar los costos de fabricación de un nuevo tipo de intercambiadores de calor en base al conocimiento de la estructura de costos histórica. Se asume que la curva de costo tiene la siguiente forma lineal: C

0

1 N

2 A

(13)

donde N es el número de tubos, y A es la superficie de la coraza. Los datos correspondientes a intercambiadores de acero templado con cabezal flotante para 0 a 500 psig (manométrico) (Shahbenderian, 1961):

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

5

Y

X1

X2

Mano de obra $

Area

Número de tubos

310

120

550

300

130

600

275

108

520

250

110

420

220

84

400

200

90

300

190

80

230

150

55

120

140

64

190

100

50

100

XT X = 10

891

3430

891

86241

349120

3430

349120

1472700

XT Y = 2135 207290 844800

Resolviendo, tenemos: = 38.177 1.164 0.209

Note que la función debe ser lineal en los coeficiente, el resto no afecta el desarrollo. Por ejemplo, la siguiente función: y

0

1 x

2 x 2

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

(14)

6

se puede resolver tomando x0 = 1, x1 = x, x2 = x2.

Diseño Factorial El cálculo de los coeficientes puede ser facilitado y su exactitud aumentada eligiendo adecuadamente los valores de X; es decir, teniendo un plan adecuado de experimentación. El diseño factorial permite alcanzar estos objetivos; por ejemplo, para la función: y

0

1 x 1

2 x 2

(15)

el método recomienda desarrollar los siguientes experimentos (ubicados en las esquinas del cuadrado unitario y en el origen): Experimento Respuesta Diseño Experimental n

Y

X1

X2

1

Y1

-1

-1

2

Y2

1

-1

3

Y3

-1

1

4

Y4

1

1

5

Y5

0

0

Se puede verificar que las sumatorias que aparecen en la ec. 7 adoptan ahora los siguientes valores:

MX 5

j1

0

(16)

j2

0

(17)

j 1

MX 5

j 1

MX X 5

0

j1 j2

(18)

j 1

MX 5

2 j1

4

(19)

j 1

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

7

MX 5

2 j2

4

(20)

j 1

con lo cual la ec. 7 simplifica a: 5 0 0 XT X

(21)

0 4 0 0 0 4

MY

j

XT Y

Y1

Y2

Y3

Y4

Y1

Y2

Y3

Y4

(22)

Una de las ventajas de este método es que para invertir la ec. 21 basta con tomar la recíproca de los elementos de la diagonal; pero más importante aún, es que la covarianza entre los coeficientes es mínima. Ejemplo. Modelo de un reactor Tenemos un reactor operando a 220°F y 3 atm. Queremos desarrollar un modelo lineal en T y P para predecir el efecto sobre la producción de perturbaciones de ± 20°F y ± 2 atm. Primero normalizamos las variables: x1

x2

T

220 20

P

(23)

3

(24)

2

Entonces, tomamos los siguientes datos: Y (Producción)

X1

X2

20.50

-1

-1

60.14

1

-1

58.89

-1

1

67.71

1

1

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

8

77.87

0

0

78.93

0

0

70.10

0

0

Los datos extras tomados en (0,0) son utilizados para medir el error involucrado en el experimento. Calculando, tenemos: 7 0 0 X

T

X

(25)

0 4 0 0 0 4

(X T X) 1

1/7

0

0

0

1/4

0

0

0

1/4

(26)

434.1 XT Y

48.46

(27)

45.96 por lo tanto = 62.02 12.11 11.49

Función objetivo En esta sección consideraremos las funciones objetivos que pueden ser planteadas con criterio económico. De acuerdo a los factores considerados, tenemos: 1. Sólo costos de operación e ingresos: corresponden a los problema de supervisor control. Los costos de capital permanecen fijos. 2. Sólo costos de capital: diseños mecánico que no afectan los costos de operación. 3. Ambos costos. Es importante analizar también algunos casos que se presentan con frecuencia dada la forma de la función objetivo. Sea la función f(x) con un punto máximo en xopt, este punto también será un Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

9

punto óptimo de la función g(x) definida como:

. f(x)

g(x)



(28)

Si . > 0, entonces: Max{g(x)}

. Max{f(x)}



(29)

Min{g(x)}

. Max{f(x)}



(30)

Si . < 0, entonces:

Un caso particular de esta situación es: Min{ f(x)}

Max{f(x)}

(31)

Otra forma de definir g(x) es: g(x)

h(f(x))



(32)

Si la derivada de h(x) es siempre positiva, entonces: Max{g(x)}

h(Max{f(x)})



(33)



(34)

Si la derivada de h(x) es siempre negativa, entonces: Min{g(x)}

h(Max{f(x)})

Ejemplo para 1. Sólo costos de operación. Una planta química produce tres productos (D, E, F) utilizando (A, B, C). Los procesos son paralelos. Las reacciones involucradas son: A A 3A

B  D B  E 2B C  F

(35)

Los datos del proceso son: Materia Prima

Disponibilidad [lb/día]

Costo [c/lb]

A

40000

1.5

B

30000

2.0

C

25000

2.5

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

10

lb reactivo / lb producto

Costo del proceso [c / lb producto]

Precio de venta [c/ lb producto]

D

2/3 A, 1/3 B

1.5

4.0

E

2/3 A, 1/3 B

0.5

3.3

F

½ A, 1/6 B, 1/3 C

1.0

3.8

 x2  3000(36) Producto

Nota: c es centavo de dolar. La masa se conserva. Se elije como función objetivo a maximizar el beneficio por día. Definiendo el vector x como la producción o consumo por día (lb/día) de (A, B, C, D, E, F), tenemos que los ingresos por días son: 0.04 x4

0.033 x5

0.038 x6

(37)

0.02 x2

0.025 x3

(38)

0.005 x5

0.01 x6

(39)

Los costos de materias primas son: 0.015 x1 y los costos de procesamiento: 0.015 x4

Entonces, los beneficios diarios son (ingresos - egresos): f(x)

0.025 x4

0.028 x5

0.028 x6

0.015 x1

0.02 x2

0.025 x3

(40)

Las restricciones surgen del balance de materia: x1 x2

0.667 x4 0.333 x4

0.667 x5 0.333 x5

0.5 x6

(41)

0.167 x6

(42)

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

11

x3

0.333 x6

(43)

El balance global no se plantea porque no es independiente. Las restricciones de disponibilidad son: 0  x1  40000

(44)

0  x3  25000

(45)

Ejemplo para 2. Sólo costos de capital. Se quiere determinar la mejor relación L/D para un tanque cilíndrico y compararla con el heurístico que plantea que la mejor relación es 3. Suponiendo que: 1. Que las dos tapas son planas. 2. Que el espesor (t) y la densidad de las paredes son constantes y que el espesor no es función de la presión. 3. Que el costo de fabricación y de materia es el mismo para las tapas como para el cuerpo y es igual a S ($/lb). 4. Que no hay desperdicios durante la fabricación. Podemos plantear las siguientes funciones objetivos, área: f1

ΠD2 2

ΠD L

ΠD2 2

ΠD L

(46)

peso: f2

!

(47)

t

costo: f3

S !

ΠD2 2

ΠD L

t

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

(48)

12

Las tres son equivalentes, por simplicidad tomaremos la primera. La restricción surge cuando especificamos el volumen:

ΠD2 L 4

V

(49)

Eliminando L, tenemos la nueva función objetivo:

ΠD2 2

f4

4 V D

(50)

Derivando y despejando, tenemos: D opt

4 V Œ

1/3

(51)

Entonces f4 es proporcional a V2/3, es decir que cumple con la conocida regla del 0.6; sin embargo: L opt

4 V Œ

1/3

(52)

dando el sorprendente resultado (L/D)opt = 1. Esto está lejos del 3 planteado por la experiencia, y se debe a que las suposiciones realizadas no se cumplen en la realidad. Ejemplo para 3. Ambos costos Se quiere determinar el espesor óptimo de una capa aislante para un recipiente que contiene un fluido caliente. El calor perdido es: A ûT x/k 1/h c

Q

(53)

donde x es el espesor del aislante. Se supone que el coeficiente externo hc es el limitante. El costo de instalación por unidad de área del aislante viene dado por: F0

F1 x

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

(54)

13

El aislante tiene un vida útil de cinco años. Suponga que tomamos un préstamo y pagamos anualmente (en cinco años) una fracción r del capital. El beneficio obtenido del calor ahorrado es Ht (10-6 $/Btu), y las horas de operación en un año son Y. El calor ahorrado es: Q

Q0

hc A

û

A ûT x/k 1/h c

T

(55)

La función objetivo que debe ser minimizada ($/año) es: f

(Q0

Q) Y Ht 10 6

(F0

F1 x) A r

(56)

Derivando y despejando, tenemos: x

opt

k

Y Ht

ûT

1/2

6

10 k F1 r

1 hc

(57)

Valor del dinero en el tiempo Hasta aquí no hemos considerado cómo varía el valor del dinero con el tiempo, frecuentemente esta variación debe ser contemplada en la función objetivo. El valor que damos al dinero varía con el tiempo, no es lo mismo recibir $100 ahora que recibirlos el año que viene. Este grado de importancia se puede expresar matemáticamente recurriendo a distintas fórmulas. Ellas serán de utilidad cuando debamos comparar costos de producción con costos de inversión ya que se debe comparar la suma que se invierte hoy con la que se pretende ganar en el futuro. Valor Futuro Es el valor futuro Fn que un determinado capital actual P tendrá en el año n con una tasa anual unitaria i capitalizando (esto es, transformar los intereses en capital como ocurre cuando se renueva un plazo fijo) q veces por año será: Fn

P 1

i q

n q

(58)

Si q = 12, tenemos un interés mensual i/q. Para q  , la ecuación anterior tiende a: Fn

P en i

(59)

Valor presente Si sólo debemos convertir un único monto futuro, la fórmula es:

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

14

Fn

P

(60)

in

1

Si son varios los montos, la fórmula a emplear es:

M1 n

P

Fk

(61)

ik

k 1

Si los montos anuales son constante e iguales a F (cuota a pagar), la ecuación anterior corresponde a una serie geométrica, por lo tanto: i)n

i (1

F

i)n

(1

P 1

(62)

Aquí, es útil definir el factor de reembolso (repayment multiplier) como: i)n

i (1

r

(1

i)n

1

(63)

por lo tanto F = r P. Para interés continuo, la fórmula es: Ejemplo. Superficie óptima de un intercambiador Queremos optimizar el área de intercambio de un generador de vapor. Se necesita enfriar una corriente de aceite, inicialmente a 400 °F, y al mismo tiempo generar vapor a 250°F y 30 psia. Las condiciones del generador son: U = 100 Btu/ (h ft2 °F) woil Cpoil = 7.5 x 104 Btu/(°F h) el precio del área es de 25 $/ft2, el beneficio proveniente del vapor generado es 2 x 10-6 $/Btu. Las horas de servicio son 8000 h/año. Suponer que se solicita un préstamo con una tasa de interés de 15% y un periodo de 10 años para su devolución en cuotas iguales. Podemos plantear la siguiente función objetivo: Max f

F

T2,A

f

2x10

6

ûH v

Wagua 8000

r I0

r 25 A [$/año]

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

(64)

(65)

15

sujeto a: w aceite Cp aceite (400

T2)

w aceite Cp aceite (400

U A (400 ln(150/(T2

T2)

ûH v

T2)

(66)

250))

w agua

(67)

donde se utilizó la diferencia de temperatura medio logarítmico (definida como la diferencia de temperaturas en el extremo caliente menos la diferencia en el extremo frío, todo dividido por el logaritmo del cociente de las dos diferencias). Utilizando estas ecuaciones y eliminando A y wagua, tenemos:

f

0.016 w aceite Cp aceite (400

T2)

25 r

waceite Cpaceite U

ln

150 T2 250

(68)

Sacando factor común w Cp, el resultado no es alterado, la solución es obtenida derivando con respecto a T2 e igualando a cero: 0

25

T2

r 1 100 T2 250

250

0.016

15.62 r

(69)

(70)

Calculando con los datos de i y n, tenemos que r es 0.2; por lo tanto, T2 = 253.1 °F, y A = 2905 ft2.

Evaluación de proyectos Todos los ejemplos considerados anteriormente fueron pequeños proyectos de inversión donde la función objetivo podía ser simplificada hasta llegar a los costos o los beneficios. Los proyectos de mayor complejidad requieren de una función objetivo acorde. En la Figura 3 se muestra el flujo de dinero de un proyecto en ejecución. Se puede apreciar como los Ingresos son derivados en Costos e Impuestos para al final producir el ingreso real en la caja denominado Cash flow.

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

16

Figura 3: Flujo de fondo de un proyecto. Dado que la DGI permite que una empresa reserve parte de sus Beneficios como fondo de ahorro para reemplazar los bienes perdurables o simplemente recuperar la inversión, se utiliza el concepto de Amortización. Existen diferentes métodos para estimar la amortización, el más común es el método de la línea recta: d

I0

Sv

(71)

n

donde Sv es el valor de rescate al cabo de la vida útil de n años del bien que se está amortizando. Generalmente, los costos de operación se estiman como un % de los costos de capital (Peters and Timmerhaus, 1980). Incluyen costos de materias primas, de mano de obra, mantenimiento, servicios, seguros, etc. Los costos de capital incluyen el precio de los equipos, su instalación, la instrumentación, aislación, etc. Una fórmula típica para estimar estos costos es: C

. S

(72)

donde S es el parámetro de tamaño, las constantes se encuentran tabuladas.

Medidas de rentabilidad Como función objetivo para un proyecto se elegirá algún índice de rentabilidad. En la siguiente tabla se muestran los índices más empleados para medir la rentabilidad de un proyecto.

Indice

Inversión inicial concentrada

Inversión distribuida

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

17

Payback period, PBP [años] Se desea el mínimo valor.

Return on investment, ROI

n 1

PBP

Se desea el máximo valor.

Valor presente neto, VAN Se desea el máximo valor.

I0

PBP

F

j

j0 n

j1

i

j

NI I0

ROI

Se desea el máximo valor.

Tasa interna de retorno, TIR = i

MI MF

F I0

VAN

M (1 F i) M (1 I i) n

i (1

F

(1

n 1

j

i)n

1

i´)n

i´ (1

j

j1

1 i´)

n

M (1 n

I0

VAN

j1

j

j0

j

M (1 I i´) n 1

Fj i´)

0

j

j0

j

j

NI: ingreso neto después de los impuestos, se parece a F. i´: es la tasa de la siguiente mejor alternativa para invertir. Tanto F (el cash flow) como I (la inversión) se toman al final del año. Los que más se utilizan son el TIR y el VAN. Cuando se optimiza un proyecto se suele maximizar el TIR o el VAN. Muchas veces se pueden hacer simplificaciones y basta con evaluar los costos o los beneficios como se hizo en los ejemplos anteriores. Ejemplo. La compra de un extrusor. Se va a comprar un extrusor a un costo de $24000, tiene una vida útil de 5 años y un valor de rescate de $4000. Los ingresos esperados son $12000 el primer año y $15000 los siguientes. Los costos de mantenimientos aumentan con la edad del equipo y son: 2000, 3000, 5000, 7000 y 8000 para cada año respectivamente. La amortización es lineal. Los impuestos son del 50%. Calcule el TIR para las siguientes situaciones: • Sin financiación. • Con el 50% de financiación con pago de capital mas interés a cinco años con una tasa del 10%. Es decir, se pagan anualmente $3170 que se reparte en intereses y capital. Caso a: Sin financiación

A. Ingresos B Costos C Beneficios (A-B)

1

2

3

4

5

12000

15000

15000

15000

15000

2000

3000

5000

7000

8000

10000

12000

10000

8000

7000

Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

18

D Amortización

4000

4000

4000

4000

4000

E Ingresos Netos antes de Impuestos (C-D)

6000

8000

6000

4000

3000

F Impuestos (05 E)

3000

4000

3000

2000

1500

G Ingresos Netos después de Impuestos (E-F)

3000

4000

3000

2000

1500

0

0

0

0

4000

7000

8000

7000

6000

9500

H Rescate I Cash Flow (D+G+H)

0

7000 1 i

8000 (1

7000

i)2

6000

i)3

(1

9500

i)4

(1

i)5

(1

24000

(73)

La solución es TIR = 0.1653 Caso b: Con 50% de financiación 1

2

3

4

5

12000

15000

15000

15000

15000

B Costos Intereses (10% de la deuda actual)

2000 1200

3000 1000

5000 790

7000 550

8000 290

C Beneficios (A-B)

8800

11000

9210

7450

6710

D Amortización

4000

4000

4000

4000

4000

E Ingresos Netos antes de Impuestos (C-D)

4800

7000

5210

3450

2710

F Impuestos (05 E)

2400

3500

2600

1720

1355

G Ingresos Netos después de Impuestos (E-F)

2400

3500

2610

7730

1355

H Principal (Cuota - Intereses)

1970

2170

2380

2620

2880

0

0

0

0

4000

4430

5330

4230

3110

6495

A. Ingresos

I Rescate J Cash Flow (D+G-H+I)

0

4430 1 i

5330 (1

4230 2

i)

(1

3110 3

i)

(1

6495 4

i)

(1

i)5

12000

(74)

El TIR es igual a 0.2708. Es mucho mejor que antes, esto se conoce como efecto palanca. Optimización y Simulación de procesos. Modelos y Función Objetivo Enrique Eduardo Tarifa - Universidad Nacional de Jujuy.

19

Optimización y Simulación de Procesos Enrique Eduardo Tarifa Facultad de Ingeniería - Universidad Nacional de Jujuy

Métodos Numéricos Introducción Los métodos analíticos imponen demasiadas restricciones a las funciones objetivos. Además, no siempre es posible resolver el sistema de ecuaciones analíticamente. Por este motivo se desarrollaron los métodos numéricos. Existen dos tipos de métodos numéricos, a saber: • Métodos directos: sólo utilizan los valores de las función objetivo. • Métodos indirectos: utilizan las condiciones necesarias, las derivadas (analíticas o numéricas) y la función objetivo. Los métodos indirectos requieren el cálculo de las derivadas primeras y segundas. Sin embargo, muchas veces obtener las derivadas es una tarea difícil, y hasta es posible que ni siquiera se conozca la forma analítica de la función objetivo. Esto plantea la necesidad de contar con métodos capaces de trabajar únicamente con los valores (experimentos) de la función objetivo. Estos son los métodos de búsqueda directa. La obtención de un valor de la función objetivo significará en algunos casos evaluar un modelo matemático, mientras que en otros significará realizar un experimento. Sea como sea, siempre será conveniente llegar al óptimo realizando la menor cantidad de evaluaciones. Esa es la misión de los métodos de búsqueda directa, a partir de los resultados de las evaluaciones realizadas, sugerirán el siguiente experimento de forma tal de aumentar la velocidad de convergencia. Es decir, que estos métodos diseñarán un adecuado plan de experiencias. El plan de experiencias puede ser secuencial o simultáneo. Cuando disponemos de un equipo por un tiempo limitado, puede ser que nos veamos obligados a realizar una serie de experimentos simultáneos. Estos experimentos son independientes, los experimentos realizados no influyen sobre la forma de realizar el siguiente. Un mejor enfoque es el plan de experiencias secuencial. Este método analiza los resultados obtenidos en un experimento para sugerir la forma de realizar el próximo. Los métodos numéricos que veremos utilizan una combinación de este tipo de plan de experiencias, pero tratan de enfatizar en el plan de experiencias secuenciales. En las siguientes secciones veremos métodos indirectos y directos para una sola variable y métodos para múltiples variables.

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

1

Métodos para una sola variable Acotación de la solución A fin de reducir el esfuerzo en la búsqueda, es conveniente limitar el espacio que el método numérico deberá explorar. En el caso de una sola variable, esto significa dar un intervalo inicial dentro del cual está el punto óptimo. Lo mejor es tener un conocimiento sobre el significado físico de las variables a fin de acotar la búsqueda a un espacio lógico; por ejemplo, las fracciones molares estarán entre 0 y 1. En el caso de no ser posible una acotación del tipo anterior, se puede emplear un método de exploración acelerada; por ejemplo, hacer variar el ln(x), o también una grilla variable como la dada por: xk

/ 2k 1

xk

1

(1)

Ejemplo (x

Min f(x)

100)2

(2)

El procedimiento de acotación utilizando la grilla variable con / = 1 da como resultado el intervalo [63, 255]. x

0

1

3

7

15

31

63

127

255

f(x)

104

9801

9409

8649

7225

4761

1369

729

2325

Métodos Indirectos Método de Newton Encuentra la raíz de la derivada de la función objetivo. Para ello, utiliza el método numérico para encontrar raíces de Newton. En este caso en particular, la ecuación de iteración es: xk

f (x k)

xk

1

(3)

f (x k)

se debe verificar continuamente si realmente el nuevo punto es mejor que el anterior. La desventaja de este método es que converge lentamente cuando la derivada segunda tiende a cero. Además, necesita conocer la derivada primera y la derivada segunda. Métodos Quasi-Newton Estos métodos imitan el método de Newton pero calculan numéricamente las derivadas, como por ejemplo: xk

1

xk

(f(xk (f(x k

h)

h)

f(xk

2 f(xk)

h))/(2 h) f(x k

h))/h 2

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

(4)

2

En este caso se utilizaron diferencias centrales para estimar las variables, pero también es posible utilizar diferencias hacia adelante ([f(x+h)-f(x)]/h) o hacia atrás ([f(x)-f(x-h)]/h). La única desventaja de este método es que requiere más evaluaciones de la función objetivo que el anterior. Método de la Secante En este caso se aplica el método para encontrar raíces de la secante. Para este caso en particular la fórmula de iteración es: x

f (a)

a (f (b)

f (a))/(b

a)

(5)

Inicialmente se elige a y b para que f’(a) tenga distinto signo de f’(b). Luego x* reemplazará a a o b según corresponda para mantener la diferencia de signos entre los extremos, y así sucesivamente. Esta variante de la secante se denomina “regula falsi”.

Métodos directos Función unimodal Se dice que una función es unimodal si la función siempre empeora cuando nos alejamos del punto óptimo. Ejemplos de funciones unimodales para el caso en que se maximiza (y*):

Figura 1: Unimodalidad en una dimensión. Criterio de eliminación de regiones Si la función es unimodal, se puede definir un criterio para eliminar regiones donde seguro el óptimo no se encuentra. Para ello necesitamos evaluar la función en dos puntos y aplicar algo de lógica. En la Figura siguiente se sobre la región eliminada para los tres casos posibles en la búsqueda de un máximo:

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

3

Figura 2: Criterio de eliminación de regiones.

Intervalo de incertitud máximo El intervalo que permanece luego de la eliminación se denomina intervalo de incertitud. Este intervalo depende de la ubicación de los puntos (x) y de la forma de la función. Para independizarlo de esta última se define el intervalo de incertitud máximo (ûMax), éste es el peor caso. Entonces, si se realizan dos experimentos x1 y x2 en el intervalo [a,b]:

û1 û2

ûMax

[a, x2] [x1, b] Max{û1,û2}

Minimax Al ser el ûMax independiente de la función evaluada, sirve para establecer un criterio de comparación entre los métodos de una sola variable denominado Minmax. Este criterio establece que cuando se comparan dos métodos, se debe determinar en cada paso el ûMax, al final se selecciona el método que tiene el mínimo ûMax. Por otra parte, podemos utilizar este criterio para diseñar un buen método. Si el criterio plantea que un buen método debe tener un ûMax mínimo, entonces podemos planear los experimentos de tal forma que ello ocurra. En otras palabras, determinamos x1 y x2 resolviendo el siguiente problema de optimización:

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

4

ûMax

Min x1,x2

En la siguiente Figura se muestra la solución de este problema en forma gráfica:

Figura 4: Criterio del Minimax. donde  = x1-x2, la solución óptima es

x

1

a b 2

0

2

En conclusión, el criterio de Minimax sugiere colocar los puntos en forma simétrica dando como resultado que los dos posibles intervalos de incertitud sean iguales al Máximo intervalo de incertitud.

Método de la dicotomía Este método plantea una secuencia de pares de experiencias simultáneas. En la siguiente Figura se muestra el correspondiente diagrama de flujo. El valor de  no puede ser muy pequeño porque puede verse afectado por los errores de redondeo de la computadora. En este trabajo se sugiere un valor igual a la mitad del intervalo final de incertitud deseado. Este procedimiento tiene una ventaja adicional y es que se puede estimar fácilmente el error de la solución; esto es, se podrá afirmar que la solución es igual x* ± .

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

5

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

6

Método del número de oro El método del número del número de oro es una mejora del anterior en el sentido que se inicia con un par de experiencias simultáneas pero luego continúa únicamente con experiencias secuenciales. Esto se logra ubicando los puntos de tal forma que el punto que queda dentro del intervalo de incertitud producido por la eliminación puede ser utilizado nuevamente. Para ello, los puntos deben cumplir con dos nuevas relaciones. La primera se obtiene de analizar la siguiente Figura:

Figura 6: Intervalos de incertidumbre.

ûj ûj 2

ûj 1

(8)

1

(9)

La otra relación es:

ûj 2 ûj 1

ûj 1 ûj

Resolviendo ambas ecuaciones se obtiene para 1 el valor número de oro.

1

5 2

1.618 que se denomina

Estas relaciones junto con el dato del intervalo inicial y la colocación en forma simétrica de las variables determinan el método por completo. En la siguiente Figura se muestra el diagrama de flujo correspondiente. Observe que esta vez el valor de  es variable y no es controlable a diferencia del método de la dicotomía. Ahora el error de la solución deberá calcularse como el Max{û*-*, *}.

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

7

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

8

Método de la serie de Fibonacci Este método, al igual que el anterior, también se inicia con dos experimentos simultáneos y luego continúa con experimentos secuenciales. La diferencia radica en que la relación de proporcionalidad empleada por el número de oro es reemplazada por un dato adicional, el  final. Se aconseja tomar este valor como igual a 0.4 / FN, donde FN es el valor N de la serie de Fibonacci definida como: Fk

Fk 1 Fk 2 F0 1 F1 1

(10)

El problema del método anterior es que debe darse como dato el valor de N que no dice mucho. Una mejora sería solicitar como dato el intervalo de incertitud final deseado y estimar a partir de allí el N. Para ello debe elegirse el menor N que hace positivo  calculado como: 

FN

û û1

FN 2

(11)

para ello basta que se cumpla: FN 

b a

û

(12)

Nuevamente el error deberá estimarse como Max{û*-*, *}. En la siguiente página se presenta el diagrama de flujo correspondiente.

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

9

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

10

Comparación Desde el punto de vista de velocidad de reducción del intervalo de incertitud, el método de la dicotomía es el mejor. Esto es porque el  se puede elegir tan pequeño que casi se elimina la mitad del intervalo por iteración. Sin embargo, requiere mayor cantidad de cálculos. El método de Fibonacci supera ligeramente en velocidad de reducción al método del número de oro, pero tiene una programación algo más compleja. En general, se recomienda el método el método del número de oro. En fórmulas, para N evaluaciones de la función objetivo, se tienen los siguientes intervalos de incertitud final:

û

û û

oro

b a FN

Fibonacci

b a

dicotomía

0

2N/2 b a 1.618N 1 FN 2 b a 0  FN FN

(13) 0

2

Interpolación cuadrática Dado tres puntos y los correspondientes valores de la función objetivo, podemos ajustar una parábola del tipo: f(x)

a x2

b x

(14)

c

Derivando e igualando a cero podemos encontrar el punto óptimo de la parábola que será una aproximación del punto óptimo buscado, la ecuación que se obtiene para el punto aproximado es: f (x)

2 a x

2

x

2

(15)

0

b 2 a

x

1 (x2 x3 ) f1 2 (x2 x3) f1

b

(16)

(x3

2

x1 ) f 2

2

(x1

2

(x3

x1) f2

(x1

2

x2 ) f 3 x2) f3

(17)

Luego, x* sustituye a uno de los puntos de la terna original de acuerdo al siguiente criterio (para minimizar): 1. Si x* está entre x2 y x3: a. Si f* < f2 y f* < f3 entonces eliminar x1. b. Si f* > f2 y f* < f3 entonces eliminar x3. c. Si f* > f2 y f* > f3 entonces hay un error numérico. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

11

2.

Si x* está entre x1 y x2: a. Si f* < f2 y f* < f1 entonces eliminar x3. b. Si f* > f2 y f* < f1 entonces eliminar x1. c. Si f* > f2 y f* > f1 entonces hay un error numérico.

Funciones de dos o más variables Introducción Ahora se verán los métodos destinados a optimizar funciones con dos o más variables. Nos concentraremos principalmente en los casos en que no existen restricciones. En tales caso sólo deberemos revisar los puntos críticos para encontrar los puntos extremos. Es decir, que primero deberemos encontrar todos los puntos que anulan todas las derivadas parciales simultáneamente. Luego estos puntos serán clasificados como mínimos o máximos de acuerdo al criterio del Hessiano. Para ello definiremos primero los Hessianos parciales de orden k como: 02f 2

0x 1

Hk

02f 0x20x1

02f 02f ... 0 x 10 x 2 0x10x k 02f 2

0x 2

...

02f 0x20x k

... 02f 0x k 0x 1

02f ... 0x k 0x 2

02f 2

0x k

Si todos los Hessianos parciales tienen determinantes positivos en un punto crítico, entonces se trata de una matriz definida positiva y el punto analizado será un mínimo. Si los determinantes tienen signos alternos, comenzando por H1 < 0, entonces el punto analizado es un máximo. Cuando algún determinante sea nulo, tenemos una matriz semidefinida y nada se puede afirmar. Por ejemplo, considere la función f(x) = x4 tiene una matriz semidefinida en el punto x=0; sin embargo, en ese punto hay un mínimo. Los problemas que enfrentaremos al tratar más de una variable son: • La propiedad de unimodalidad falla y su definición es muy restringida. • Problemas de aristas. • Imposibilidad de establecer criterios de eficiencia independientes de la función optimizada. Definición de unimodalidad Vamos a definir tres tipos de unimodalidad, ellos son: • Unimodal: Una función es unimodal si por cualquier punto del dominio pasa una trayectoria que conduce al punto óptimo, y la función siempre mejora a medida que nos acercamos al óptimo por dicha trayectoria. • Fuertemente unimodal: Si toda recta que pasa por el óptimo es unimodal. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

12

•

Linealmente unimodal: Si cualquier recta que pasa por el dominio es unimodal.

Figura 9: Unimodalidad. Un problema que surge de la gran cantidad de cálculos que debe realizarse es que el ruido provocado por el error de redondeo crece en magnitud. Esto dar origen a falsos óptimos locales cuando el paso de búsqueda se reduce mucho. Problema de aristas El problema de aristas es un problema muy común y surge cuando la función depende muy poco de alguna variable, pero la posición del óptimo depende fuertemente de todas las variables. Los métodos deben ser diseñados para superar este problema o de lo contrario reportarán puntos óptimos falsos. En la siguiente Figura se ilustra cuál es el problema.

Figura 10: Problema de arista. Criterio de eficiencia Un problema que tenemos que enfrentar ahora es la falta de un criterio de eficiencia general. Esto dificulta tanto la comparación entre métodos como el diseño de un método eficiente. Se puede afirmar que la eficiencia de un método depende de: • El problema a optimizar. • El punto de partida y los intervalos iniciales de prueba. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

13

•

La computadora y precisión utilizada.

Podemos decir que un buen método debe surgir de resolver el siguiente problema de optimización (planteado para el caso de búsqueda de máximos): ) f(x ) M f(xCálulos N 1

Max

.,dx,N

xk 1

k 1

xk

k 1

k

.k dx k

donde . es el paso que dará en la dirección dx. Restricciones Las restricciones pueden transformar un problema que era unimodal en bimodal (ver la Figura 11). Los métodos que veremos son para problemas sin restricciones, sin embargo es posible aplicarlos si se sigue la siguiente secuencia de pasos: •

• • • • •

Si el problema tiene restricciones fuertes (de igualdad), elimine cada una de ella despejando una variable y reemplazando ésta en el resto de las ecuaciones. Repita el procedimiento hasta haber eliminado todas las restricciones fuertes y sus correspondientes variables. Optimice sin considerar las restricciones débiles (desigualdades). Si todas las restricciones se cumplen, entonces ya encontró el óptimo. Identifique las restricciones que no se cumplen. Tome de a una las que no se cumplen y transfórmelas en igualdad. Resuelva nuevamente el problema. Si aún no encontró el óptimo, tome de a dos las restricciones que no se cumplen y resuelva el problema nuevamente, y así sucesivamente.

Figura 11: Región bimodal. En la siguiente Figura se ilustra un caso en que tanto h1 como h2 no se cumplen, pero basta con tomar como igualdad a h2 para encontrar el óptimo. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

14

Figura 12: Optimo con restricciones. Método de una variable a la vez Este método es funciona bien cuando las curvas de nivel son circulares, siempre llega al óptimo en dos etapas; pero fracasa cuando se encuentra con curvas de nivel elípticas que originan aristas.

Figura 13: Método de una variable a la vez. Método del máximo gradiente Para interpretar gráficamente este método necesitamos determinar la relación entre las curvas de nivel y el gradiente de una función. Sea la función f(x,y), entonces podemos definir una curva de nivel como f(x,y) = k. Diferenciando se tiene: 0f dx 0x

0f dy 0y

0

Entonces, la dirección tangente a la curva es: dy dx

fx fy

Es decir que el /f es perpendicular a la curva de nivel porque su pendiente es igual a la inversa cambiada de signo de la pendiente de la línea tangente a la curva de nivel. Otro resultado interesante de probar es que el gradiente es la dirección de mayor cambio. Para ello resolvemos el siguiente problema: Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

15

Max .

df ds

df ds

fx cos(.)

f y sin(.)

derivando e igualando a cero: df s d.

f x ( sin(.)) sin(.) cos(.)

tan(.)

fy cos(.) fy fx

Es decir, /f es la dirección de mayor pendiente positiva. El método del máximo gradiente utiliza la dirección /f para realizar la búsqueda. Si las curvas de nivel son circulares, el método llega al óptimo en un solo paso. En la Figura 14 se muestra esto junto con dos variantes del método, en la primera se optimiza a lo largo de la dirección de /f, en el otro se da un paso determinado de antemano que luego se va reduciendo. Esta última variante es la que se ilustra con el diagrama de flujo de la página siguiente.

Figura 14: Método del gradiente.

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

16

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

17

Métodos empíricos Están basados en la experiencia y son los más aconsejables. No requieren el cálculo de las derivadas parciales. En las páginas siguientes se dan los diagramas de flujo de dos métodos. El primero es el Simplex y el otro es el Nelder y Mead, el segundo es el de Hooke y Jeeves. La siguiente Figura muestra los pasos elementales que emplea el método de Nelder y Mead.

Figura 16: Pasos en el método de Nelder y Mead. En la siguiente Figura se ilustra algunos pasos del método de Hooke y Jeeves.

Figura 17: Pasos en el método de Hooke y Jeeves.

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

18

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

19

Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

20

Direcciones conjugadas La experiencia demostró que la búsqueda a través de direcciones conjugadas es una buena alternativa. Dos direcciones p y q con respeto a la matriz cuadrada definida positiva Q son conjugadas si: pT Q q

(24)

0

En general, un conjunto de n direcciones linealmente independientes serán conjugadas con respecto a la matriz cuadrada definida positiva Q si: (s i)T Q (s j)

0

0  i g j  n

1

(25)

Si Q = I tenemos el caso especial de ortogonalidad. En optimización, Q = H (el Hessiano). Si la función objetivo de n variables es cuadrática, Fletcher y Reeves (1964) demostraron que se puede alcanzar el óptimo en n etapas, si en cada etapa se optimiza en la correspondiente dirección conjugada. Las direcciones conjugadas no son únicas, se debe especificar una dirección inicial para que el resto queden determinadas. Si la función fuese cuadrática, el desarrollo de Taylor sería: f(x)

/ s k)

f(x k

f(x k)

1 (/ s k)T H(x k) (/ s k) 2

/Tf(x k) / s k

(26)

para obtener el paso que optimiza en la dirección sk, derivamos e igualamos a cero: / s k)

df(x k d/

0

/Tf(x k) s k

(s k)T H(x k) (/ s k)

(27)

despejando: /opt

/Tf(x k) s k (s k)T H(x k) (s k)

(28)

Ejemplo Minimizar: f(x)

2

2 x1

2

x2

3

(29)

partiendo del punto (1,1) y la dirección inicial (-4,-2). Primero, calculamos el Hessiano y utilizamos la definición de direcciones conjugadas para Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

21

encontrar la otra dirección: 1

[ 4 2]

4 0 s1

(30)

0

0 2 s1 2

Debido a que no hay una única solución, tomamos s11 = 1 y determinamos la componente restante: 1 [ 16 4]

0

1

s2

(31)

por lo tanto, la dirección conjugada es (1,-4). Verifiquemos ahora si se alcanza el óptimo en dos etapas. Optimizando en la dirección (-4,-2), de la ec. 25, tenemos: 4

[4 2] /0

[ 4 2]

2

4 0

4

0 2

2

20 72

(32)

Luego, el nuevo x es: x1

x0

/0 s 0

1 1

20 72

4

0.1111

2

0.4444

(33)

A continuación se optimiza en la dirección (1, -4): 1

[ 0.4444 0.8888] /1

[1 4]

4

4 0

1

0 2

4

0.111

(34)

damos el nuevo paso: x2

x1

/1 s 1

0.1111 0.4444

0.1111

1

0

4

0

(35)

que concuerda con el punto óptimo. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

22

Cálculo de la dirección conjugada sin usar derivadas Partiendo del punto x0, localice el punto xa que minimice la función en la dirección s. Luego, parta de otro punto x1 distinto del primero y localice el punto xb que minimice la función en la misma dirección s. La dirección (xb - xa) es conjugada a s. La Figura 20 muestra esta situación.

Figura 20: Determinación de direcciones conjugadas.

Métodos indirectos Los métodos directos fueron los primeros utilizados para problemas sin restricciones, y son frecuentemente utilizados en los procesos químicos debido a que son fáciles de entender y de implementar. Sin embargo, no son tan eficientes y robustos como los modernos métodos indirectos. Método del Gradiente Es el método que vimos previamente, también conocido como el método de Cauchy, se trata de un método indirecto de primer orden porque utiliza el gradiente. Este método tiene el inconveniente que a medida que nos acercamos al óptimo (o a un punto crítico) su velocidad de convergencia disminuye notoriamente (el gradiente tiende a anularse). Sin embargo, tiene una buena velocidad en las primeras etapas, cuando estamos lejos del óptimo. Siempre se debe verificar que el punto encontrado no sea un punto de ensilladura. Otra desventaja del método es su elevada sensibilidad al escalado de la función objetivo provocando una baja velocidad de convergencia y fuertes oscilaciones. Gradiente conjugado Este método fue diseñado por Fletcher y Reeves (1964). Utiliza como dirección de búsqueda una combinación lineal del gradiente actual con el gradiente de la iteración anterior. Es una notable mejora con respecto al método original, porque para funciones cuadráticas se puede demostrar que las sucesivas direcciones de búsqueda son conjugadas, por lo tanto el óptimo se encuentra en n pasos, donde n es el número de variables. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

23

En resumen, el método sigue los siguientes pasos (para minimizar):

/f(x 0)

s0 xk 1

sk 1

/k s k

xk

/f(x k 1)

sk

(36)

(37)

/Tf(x k 1) /f(x k 1) /Tf(x k) /f(x k)

(38)

donde el paso /k se hace minimizando en la dirección sk. A medida que se acerca al óptimo, se parece cada vez más al método del Gradiente, y con los mismos inconvenientes. Método de Newton Es un método indirecto de segundo orden porque utiliza una aproximación de Taylor de segundo orden: f(x)  f(x k)

1 (ûx k)T H(x k) (ûx k) 2

/Tf(x k) ûx k

(39)

Derivando e igualando a cero, tenemos:

/f(x)  /f(x k)

H(x k) ûx k

0

(40)

La ecuación anterior representa un sistema de ecuaciones donde la incógnita es el vector: ûx k

xk 1

xk

[H(x k)] 1 /f(x k)

(41)

Note que el método establece tanto la dirección como la magnitud del paso. Nos debemos asegurar que H permanece definida positiva para asegurar un mínimo. El método tiene la gran ventaja de comportarse mejor a medida que más cerca está del óptimo, al contrario del método del Gradiente.

Forzando a la matriz H Como es importante que el Hessiano se mantenga definido positivo, Marquardt (1963), Levenberg (1944) sugirieron utilizar una matriz modificada que se mantuviera definida positiva: H (x)

H(x)

 I

(42)

donde  se elige suficientemente grande como para que la nueva matriz sea definida positiva cuando el Hessiano no lo sea. Optimización y Simulación de Procesos. Métodos Numéricos. Enrique Eduardo Tarifa - Universidad Nacional de Jujuy

24

Una alternativa es: [H (x)] 1

 I

H 1(x)

(43)

donde  también es un escalar suficientemente grande para mantener la característica de definida positiva. Para estimar el valor de  se puede utilizar el valor del autovalor más negativo, esto es:

 > min(.i)

(44)

Note que si este factor es demasiado grande, entonces I predomina y el método tiende al método del Gradiente. Método de la Secante Este método sólo utiliza los valores de la función y del gradiente; luego, la matriz Hessiana es aproximada por una combinación de los otros dos. Estos métodos también se conocen como métodos quasi-Newton. Criterios de terminación Los siguientes criterios se recomiendan para evitar problemas de escalado: f(x k 1)

f(x k) < f(x k)











1

(45)

2

(46)







o si la función tiende a cero: f(x k 1)

f(x k)