UNIVERSIDAD SIMON BOLIVAR DECANATO DE ESTUDIOS DE POST-GRADO DOCTORADO INTERDISCIPLINARIO EN CIENCIAS

NUEVOS OPERADORES GENÉTICOS PARA EL DISEÑO DE CONTROLADORES CON MULTIPLES OBJETIVOS Y RESTRICCIONES NO-CONVEXAS

Tesis Doctoral presentada ante la Universidad Simón Bolívar por

Gustavo Adolfo Sánchez Hurtado

como requisito parcial para optar al grado académico de

DOCTOR EN CIENCIAS

con la asesoría de la profesora

Minaya Villasana

Mayo, 2011

v

AGRADECIMIENTOS Quisiera manifestar mi agradecimiento a las personas que han colaborado para realizar este trabajo de investigación. En primer lugar a la profesora Minaya Villasana (USB) sin cuya valiosa asesoría hubiese sido imposible llevar a término el proyecto. En segundo lugar a los profesores Miguel Strefezza y Ubaldo García Palomares (USB) por sus múltiples revisiones y comentarios. A los profesores Oliver Schüetze y Carlos Coello Coello por su cordial recibimiento y ayuda durante el tiempo en que estuvimos trabajando en conjunto en las instalaciones del CINVESTAV, México. A los profesores Carlos González y Orlando Reyes (USB) por el apoyo recibido, tanto en el ámbito personal y profesional. Por último, a mi esposa, mis hijos, y a toda mi familia por el tiempo, la paciencia y la comprensión que me han brindado.

vi

UNIVERSIDAD SIMON BOLIVAR DECANATO DE ESTUDIOS DE POST-GRADO DOCTORADO INTERDISCIPLINARIO EN CIENCIAS

NUEVOS OPERADORES GENÉTICOS PARA EL DISEÑO DE CONTROLADORES CON MULTIPLES OBJETIVOS Y RESTRICCIONES NO-CONVEXAS

Por: Gustavo Adolfo Sánchez Hurtado Carnet: 99-80737 Tutora: Prof. Minaya Villasana Fecha: Mayo, 2011

RESUMEN En determinadas circunstancias, el diseño de un controlador automático debe ser formulado como un problema de optimización con múltiples objetivos y restricciones noconvexas. Desafortunadamente, no siempre estos problemas pueden resolverse mediante algoritmos deterministas que garanticen la convergencia hacia una solución óptima. En estos casos, los algoritmos aleatorios (y en particular los genéticos) pueden representar una alternativa cuando se desea obtener mejores resultados. En efecto, la mayoría de los reportes publicados sobre el tema confirman de manera experimental la eficacia de estos algoritmos para resolver problemas de diseño de controladores automáticos. No obstante, aún persisten en esta área de investigación numerosos aspectos por mejorar e importantes problemas que considerar. En este contexto, la principal contribución del presente trabajo consiste en la propuesta y el análisis de nuevos operadores genéticos, especialmente concebidos para intentar obtener mejores aproximaciones de los Conjuntos de Pareto, en comparación con otros métodos competitivos citados frecuentemente en la literatura. Palabras Claves: operadores de variación, diseño de controladores, optimización multiobjetivo, algoritmos genéticos, restricciones no convexas.

ÍNDICE GENERAL APROBACIÓN DEL JURADO

ii

AGRADECIMIENTOS

v

RESUMEN

vi

ÍNDICE DE TABLAS

xi

ÍNDICE DE FIGURAS

xiv

ÍNDICE DE ALGORITMOS

xv

ÍNDICE DE SÍMBOLOS Y ABREVIATURAS

xvi

1 INTRODUCCIÓN 1.1 Antecedentes . . . . . . . . . . . . . . . . . . . . . . 1.2 Planteamiento del problema de investigación . . . . . 1.2.1 Consideración de las restricciones de diseño. . 1.2.2 Representación del espacio de búsqueda . . . . 1.2.3 Integración con operadores de búsqueda local 1.3 Estructura de la tesis . . . . . . . . . . . . . . . . . . 2 SISTEMAS LINEALES REALIMENTADOS 2.1 Generalidades . . . . . . . . . . . . . . . . . . 2.2 Índices de desempeño . . . . . . . . . . . . . . 2.2.1 Índices de estabilidad . . . . . . . . . . 2.2.2 Índices de respuesta temporal . . . . . 2.2.3 Rechazo de perturbaciones . . . . . . . 2.2.4 Estabilidad robusta . . . . . . . . . . . 2.3 Controladores PID . . . . . . . . . . . . . . . 2.3.1 Método de Ziegler y Nichols . . . . . . 2.3.2 Método de Hohenbichler y Abel . . . . 2.4 Reubicación de polos . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

. . . . . . . . . .

. . . . . .

1 1 8 8 9 11 12

. . . . . . . . . .

15 15 18 18 19 21 22 22 23 24 27

2.5 Diseño de controladores mediante LMIs 2.5.1 Problema H2 . . . . . . . . . . 2.5.2 Problema H∞ . . . . . . . . . . 2.6 Resumen del capítulo . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

3 OPTIMIZACIÓN MULTI-OBJETIVO 3.1 Definiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Estrategias para resolver un problema multi-objetivo . . . . . 3.2.1 Formulación de preferencias “a priori” . . . . . . . . . 3.2.2 Formulación de preferencias de manera interactiva . . . 3.2.3 Formulación de preferencias “a posteriori” . . . . . . . 3.3 Algoritmos genéticos multi-objetivo . . . . . . . . . . . . . . . 3.3.1 Operadores de variación . . . . . . . . . . . . . . . . . 3.3.2 Esquemas de selección . . . . . . . . . . . . . . . . . . 3.3.3 Métodos de archivo . . . . . . . . . . . . . . . . . . . . 3.3.4 Tratamiento de las restricciones . . . . . . . . . . . . . 3.3.5 Multi-Objective Genetic Algorithm (MOGA) . . . . . . 3.3.6 Non-dominated Sorting Genetic Algorithm (NSGA-II) 3.3.7 Strength Pareto Evolutionary Algorithm (SPEA2) . . . 3.4 Evaluación de desempeño . . . . . . . . . . . . . . . . . . . . 3.4.1 Índices de desempeño . . . . . . . . . . . . . . . . . . . 3.4.2 Pruebas estadísticas . . . . . . . . . . . . . . . . . . . 3.5 Resumen del capítulo . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . .

4 DISEÑO DE CONTROLADORES PID MULTI-OBJETIVO 4.1 Antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Método de diseño basado en el cálculo exacto de la región de estabilidad 4.3 Método de diseño basado en la aproximación la región de estabilidad . . 4.4 Resumen del capítulo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 BÚSQUEDA MULTI-OBJETIVO LOCAL 5.1 Antecedentes . . . . . . . . . . . . . . . . . 5.2 Ascenso de Colina con Paso Lateral (ACPL) 5.2.1 Caso sin información de gradientes . 5.2.2 Caso con información de gradientes . 5.2.3 Tratamiento de restricciones . . . . . 5.3 Evaluación del operador ACPL . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . .

viii 30 31 34 35

. . . . . . . . . . . . . . . . .

36 36 39 39 41 42 42 44 46 47 48 49 51 54 57 57 60 61

63 . 63 . 64 . 75 . 81

. . . . . .

82 82 84 84 89 90 91

5.3.1 Problema convexo sin restricciones . . . . . . . . . 5.3.2 Problema convexo con restricciones de intervalo . . 5.3.3 Problema no convexo con restricciones de intervalo 5.4 Sensibilidad del ACPL . . . . . . . . . . . . . . . . . . . . 5.5 Integración SPEA2-ACPL . . . . . . . . . . . . . . . . . . 5.5.1 Problema convexo sin restricciones . . . . . . . . . 5.5.2 Problema no convexo con restricciones de intervalo 5.6 Resumen del capítulo . . . . . . . . . . . . . . . . . . . . . 6 PROBLEMA H2 /H∞ 6.1 Formulación . . . . . . . 6.2 Antecedentes . . . . . . 6.3 Solución mediante LMIs 6.4 Método MOPPEA . . . 6.5 Resultados numéricos . . 6.6 Resumen del capítulo . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

. . . . . .

. . . . . . . .

ix 91 95 99 100 105 106 108 110

. . . . . .

112 112 113 115 116 120 128

7 CONCLUSIONES

129

APÉNDICE

133

REFERENCIAS

137

ÍNDICE DE TABLAS 1.1 Promedio y desviación standard del valor mínimo alcanzado por cada algoritmo, calculados sobre un total de 100 ejecuciones. Problemas P1 y P2 . .

11

2.1 Frecuencias singulares y ecuaciones de las rectas que delimitan la región de estabilidad para KP = −4.4, modelo (2.26) . . . . . . . . . . . . . . . . . .

26

3.1 Resumen de ventajas y desventajas de los algoritmos MOGA, NSGA-II y SPEA2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Resumen de indicadores de desempeño . . . . . . . . . . . . . . . . . . . . 3.3 Ejemplo de valores obtenidos para el indicador I correspondientes a 7 ejecuciones de los algoritmos A1 y A2 . . . . . . . . . . . . . . . . . . . . . . 4.1 Proporción de individuos factibles para distintos valores de L, L . . . . . . 4.2 Configuración utilizada para los algoritmos A1 y A2. Problema de diseño PID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Resultados obtenidos para el indicador de cobertura. Algoritmos A1 y A2 . 4.4 Resultados de los indicadores ESS y DMAX. Algoritmos A1 y A2 . . . . . 4.5 P-valores correspondientes a las pruebas de Wilcoxon para determinar las diferencias que son estadísticamente significativas. Algoritmos A1, A2 . . . 4.6 Configuración utilizada por SPEA2 durante la primera etapa del método de aproximación de la región factible . . . . . . . . . . . . . . . . . . . . . 4.7 Resultados del indicador de cobertura. Algoritmos A1, A2 y A3 . . . . . . 4.8 Resultados de los indicadores N, ESS y DMAX. Algoritmos A1, A2 y A3 . 4.9 P-valores correspondientes a las pruebas de Wilcoxon para determinar las diferencias que son estadísticamente significativas. Algoritmos A1, A2 y A3

56 60 60 66 72 73 73 73 78 80 80 80

5.1 Parámetros que deben ser definidos por el usuario en los casos con y sin información de gradientes . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2 Valores utilizados para los parámetros ACPL . . . . . . . . . . . . . . . . . 92 5.3 Valores utilizados para los parámetros ACPL. Problema P3 . . . . . . . . 100 5.4 Sensibilidad del operador ACPL1 ante variaciones en el vector r . . . . . . 101 5.5 Sensibilidad del operador ACPL1 ante variaciones en el vector Nnd . . . . 103

5.6 5.7 5.8 5.9 5.10

Configuración utilizada para resolver el problema P4 . Indicadores GD/IGD obtenidos para el problema P4 Configuración utilizada para las pruebas del problema Indicadores GD/IGD obtenidos para el problema P3 Indicadores GD/IGD obtenidos para el problema P5

6.1 Datos correspondientes a los problemas COMP le ib 6.2 Parámetros utilizados por SPEA2 y SPEA2-ACPL 6.3 Resultados del indicador de cobertura . . . . . . . . 6.4 Resultados de los indicadores ESS y DMAX . . . .

. . . .

. . . . P3 . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

xi 107 108 108 109 109

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

122 126 127 127

. . . .

ÍNDICE DE FIGURAS 1.1 Modelo a lazo cerrado en tiempo discreto. Controlador de dos grados de libertad. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Sistema a lazo cerrado con controlador PID . . . . . . . . . . . . . . . . . 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10

2.11 2.12

Sistema LTI a lazo abierto . . . . . . . . . . . . . . . . . . . . . . . . . . . Interconexión Planta/Controlador (lazo cerrado) . . . . . . . . . . . . . . . Diagrama de polos de un sistema con índices de estabilidad ns = 2, qa = 3 Diagrama de polos de un sistema con índices de estabilidad ns = 1, qa = 2 Indices de desempeño a partir de la respuesta ante un escalón unitario . . . Lazo de realimentación con incertidumbre no estructurada ∆ . . . . . . . . Lazo de control con PID ideal . . . . . . . . . . . . . . . . . . . . . . . . . Respuesta del sistema con los parámetros KD , KP , KI de acuerdo con las fórmulas de Ziegler y Nichols . . . . . . . . . . . . . . . . . . . . . . . . . . Límites de la región de estabilidad para KP = −4.4 para el modelo (2.26) . Respuesta al escalón del sistema a lazo cerrado para el modelo (2.26) correspondiente al controlador PID caracterizado por KP = −4.4, KI = 2.2 y KD = −10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Diagrama de polos (×) y ceros (◦) a lazo abierto para el modelo (2.26) . . Respuesta al escalón del sistema a lazo cerrado, para el modelo (2.26). Controlador diseñado mediante reubicación de polos. . . . . . . . . . . . .

3.1 Esquema general de un algoritmo genético . . . . . . . . . . . . . . . . . . 3.2 Operador discreto de cruce de un punto. El punto de cruce se encuentra en la mitad del cromosoma. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Operador discreto de cruce de dos puntos. . . . . . . . . . . . . . . . . . . 3.4 Ejemplo de cálculo de rango para el algoritmo MOGA . . . . . . . . . . . . 3.5 Ejemplo de cálculo de rango para el algoritmo NSGA-II . . . . . . . . . . . 3.6 Ejemplo de cálculo de crowding distance para el algoritmo NSGA-II . . . . 3.7 Ejemplo de cálculo de S y R para el algoritmo SPEA2 . . . . . . . . . . . 3.8 Ejemplo de cálculo de la distancia σ para el algoritmo SPEA2 . . . . . . .

6 9 15 18 19 19 20 22 23 24 27

28 30 30 43 45 45 51 52 53 55 56

xiii 3.9 Ejemplo de dos aproximaciones con distintos valores para los indicadores GD, IGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Ejemplo de cálculo para el indicador DMAX . . . . . . . . . . . . . . . . 3.11 Cálculo de los rangos para la prueba Wilcoxon . . . . . . . . . . . . . . . . 4.1 Espacio de controladores PID estabilizantes correspondiente al modelo G0(s) determinado mediante el método de Hohenbichler y Abel (2008) . . . . . . 4.2 Selección de un punto al azar en el interior de un polígono convexo . . . . 4.3 Lazo con controlador PID . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Ejemplo de las aproximaciones obtenidas por A1 y A2 . . . . . . . . . . . . 4.5 Evolución de las aproximaciones de los frentes. Algoritmo A1. Ngen es el número de generación. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Evolución de las aproximaciones de los frentes. Algoritmo A2. Ngen es el número de generación. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Perturbación y ruido que actúan sobre el lazo . . . . . . . . . . . . . . . . 4.8 Salida del lazo de control, afectada por perturbación y ruido . . . . . . . . 4.9 Ejemplo de aproximación obtenida al final de la primera etapa . . . . . . . 4.10 Ejemplo de las aproximaciones que se obtienen mediante A1, A2 y A3 . . . 5.1 Situación que se presenta cuando el punto de prueba se encuentra lejos de un COLP. Cono de descenso amplio (región sombreada). . . . . . . . . . . 5.2 Situación que se presenta cuando el punto de prueba se encuentra cerca de un COLP. Cono de descenso estrecho (región sombreada). . . . . . . . . . . 5.3 Verdadero conjunto de Pareto para el problema P1 . . . . . . . . . . . . . . 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15

Verdadero frente de Pareto para el problema P1 . . . . . . . . . . . . . . . Aproximación del conjunto de Pareto generada por ACPL1. Problema P 1 Aproximación del frente de Pareto generada por ACPL1. Problema P 1 . . Aproximación del conjunto de Pareto generada por ACPL2. Problema P 1 Aproximación del frente de Pareto generada por ACPL2. Problema P 1 . . Verdadero conjunto de Pareto para el problema P2 . . . . . . . . . . . . . . Verdadero frente de Pareto para el problema P2 . . . . . . . . . . . . . . . Aproximación del conjunto de Pareto generada por ACPL1. Problema P 2 Aproximación del frente de Pareto generada por ACPL1. Problema P 2 . . Aproximación del conjunto de Pareto generada por ACPL2. Problema P 2 Aproximación del frente de Pareto generada por ACPL2. Problema P 2 . . Ejemplo de aproximaciones del frente de Pareto generadas por ACPL1 y ACPL2. Problema P 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58 59 61

66 67 68 71 72 72 74 74 78 79

85 85 92 93 93 94 94 95 96 96 97 97 98 98 99

xiv

5.16 Aproximación obtenida mediante ACPL2 con ε1 =[0.01, 0.01] . . . . . . . . 101 5.17 Aproximación obtenida mediante ACPL2 con ε2 =[0.1, 0.1] . . . . . . . . . 102 5.18 Aproximación obtenida mediante ACPL2 con ε3 =[1, 1] . . . . . . . . . . . 102 5.19 5.20 5.21 5.22

6.1 6.2 6.3 6.4 6.5 6.6 6.7

Aproximación obtenida mediante ACPL2 con εP = 0.01 . . . . . . . . . . . Aproximación obtenida mediante ACPL2 con εP = 0.1 . . . . . . . . . . . Aproximación obtenida mediante ACPL2 con εP = 1 . . . . . . . . . . . . Aproximaciones obtenidas por ACPL1 para distintos valores del parámetro T ol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

105

Modelo lineal a lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . Representación genética propuesta. Método MOPPEA . . . . . . . . . Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema AC1 . Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema AC6 . Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema NN10 Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema WEC1 Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema AC9 .

113 118 123 124 124 125 125

. . . . . . .

. . . . . . .

103 104 104

ÍNDICE DE ALGORITMOS 3.1 3.2 3.3 3.4 3.5 5.1 5.2 6.1

Extracción del sub-conjunto de soluciones no-dominadas . . . . . . . . . . . 42 Multi-Objective Genetic Algorithm (MOGA) . . . . . . . . . . . . . . . . . 50 Asignación de rango NSGA-II . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Non-Dominated Sorting Genetic Algorithm (NSGA-II) . . . . . . . . . . . 53 Strength Pareto Evolutionary Algorithm 2 . . . . . . . . . . . . . . . . . . 54 Operador ACPL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Retorno a la región factible para el operador ACPL1 . . . . . . . . . . . . . 90 Generación de una aproximación del Frente de Pareto para el problema mixto H2 /H∞ mediante LMIs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.2 Generación de un vector de polos aleatorio . . . . . . . . . . . . . . . . . . . 121

ÍNDICE DE SÍMBOLOS Y ABREVIATURAS LQR

Linear Quadratic Regulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

x(t)

Vector de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

u(t)

Vector de señales manipuladas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

H2

Norma 2 en el espacio de funciones complejas . . . . . . . . . . . . . . . . . . . . . . 3

H∞

Norma infinita en el espacio de funciones complejas . . . . . . . . . . . . . . . . 3

LMI

Linear Matrix Inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

PID

Proportional Integral Derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

MOEA

Multi-Objective Evolutionary Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

MOGA

Multi-Objective Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

G(s)

Matriz o función de transferencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

SPEA2

Strength Pareto Evolutionary Algorithm 2 . . . . . . . . . . . . . . . . . . . . . . . . . 5

MRCD

Multi-objective Robust Control Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

NSGA − II Non-Dominated Sorting Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 7 λ

Vector de auto-valores de un sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

MOPSO

Multiple Objective Particle Swarm Optimization . . . . . . . . . . . . . . . . . . . 7

MOPPEA

Multi-Objective Pole Placement with Evolutionary Algorithms . . . . 13

COMP leib COnstraint Matrix-optimization Problem library . . . . . . . . . . . . . . . . . .13 y(t)

Vector de señales medidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

w(t)

Vector de señales exógenas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

z(t)

Vector de señales controladas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Ln2 [0, +∞)

Espacio lineal de señales de energía finita . . . . . . . . . . . . . . . . . . . . . . . . . 15

R

Conjunto de los números reales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

C

Conjunto de los números complejos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

nu

Dimensión del vector de señales manipuladas . . . . . . . . . . . . . . . . . . . . . 15

nw

Dimensión del vector de señales exógenas . . . . . . . . . . . . . . . . . . . . . . . . . 15

xvii ny

Dimensión del vector de señales medidas . . . . . . . . . . . . . . . . . . . . . . . . . 15

nz

Dimensión del vector de señales controladas . . . . . . . . . . . . . . . . . . . . . . 15

n

Dimensión del vector de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

|X|

Número de elementos del conjunto X . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

ncl

Grado del sistema a lazo cerrado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

ns

Número de polos estables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

qa

Abscisa del polo que se encuentra más a la derecha . . . . . . . . . . . . . . . 18

ess

Error en estado estacionario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

tr

Tiempo de subida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

ts

Tiempo de asentamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

Mp

Porcentaje de sobre-elongación máxima . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

R(s)

Conjunto de funciones complejas racionales de coeficientes reales . . 21

KP

Constante proporcional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

KI

Constante integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

KD

Constante derivativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Tu

Período de las oscilaciones críticas. Método de Ziegler-Nichols . . . . . 24

Ku

Ganancia de la cadena directa. Método de Ziegler-Nichols . . . . . . . . 24

POM

Problema de optimización Multi-Objetivo . . . . . . . . . . . . . . . . . . . . . . . . 36

X

Espacio factible de variables de decisión . . . . . . . . . . . . . . . . . . . . . . . . . . 36

n

Número de variables de decisión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

M

Número de funciones objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

x

Vector de variables de decisión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

F

Espacio Objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

z1 z2

z 1 domina a z 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

C(x,b)

Cono de diversidad de tipo b en x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

COGP

Conjunto Óptimo Global de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

COLP

Conjunto Óptimo Local de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

FOGP

Frente Óptimo Global de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

FOLP

Frente Óptimo Local de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

∇f

Gradiente de f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

xviii IWTP

Interactive Weighted Tchebycheff Procedure . . . . . . . . . . . . . . . . . . . . . . 41

µs

Número total de individuos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

λs

Número de hijos generados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

GD

Distancia generacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

IGD

Distancia generacional inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

DMAX

Distancia máxima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .58

C

Cobertura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

ESS

Espaciado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

H0

Hipótesis nula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

H1

Hipótesis alternativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

α

Nivel de significancia de la prueba de Wilcoxon . . . . . . . . . . . . . . . . . . . 61

L, L

Intervalo de variación para cada parámetro PID . . . . . . . . . . . . . . . . . . 65

Ns (L, L)

Número de individuos factibles muestreados en L, L . . . . . . . . . . . . . . 65

K

Conjunto de controladores estabilizantes . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Vmax

Volumen máximo del espacio de búsqueda . . . . . . . . . . . . . . . . . . . . . . . . 77

B(x0 , r)

Vecindario de definido alrededor de x0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Nnd

Vector de límites para el cálculo de direcciones de diversidad . . . . . . 88

ε

Parámetros para calcular la longitud del paso lateral . . . . . . . . . . . . . . 88

ab

Dirección acumulada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Niter εP

Número de puntos muestreados en el vecindario del punto de origen 88 Tolerancia para decidir si un punto pertenece a P∗ . . . . . . . . . . . . . . . . . 89

T ol

Tolerancia para decidir si un punto es factible . . . . . . . . . . . . . . . . . . . . 90

rmax

Cota para delimitar un semi-plano complejo . . . . . . . . . . . . . . . . . . . . . 120

imax

Cota para delimitar un semi-plano complejo . . . . . . . . . . . . . . . . . . . . . 120

H2 /H∞

Número de iteraciones del algoritmo basado en LMIs . . . . . . . . . . . . 122

∆γ 2

Define el incremento en la cota de la inecuación H2 . . . . . . . . . . . . . . 122

∆γ inf

Define el incremento en la cota de la inecuación Hinf . . . . . . . . . . . . 122

γ min

Define la mínima cota de la inecuación H2 . . . . . . . . . . . . . . . . . . . . . . . 122

γ max

Define la máxima cota de la inecuación Hinf . . . . . . . . . . . . . . . . . . . . 122

CAPÍTULO I INTRODUCCIÓN

El concepto de sistema es tal vez uno de los más fértiles en la historia de la ciencia y la tecnología (Bertalanffy, 1950). Es posible definir un sistema real (compuesto por equipos, maquinaria, objetos, etc) como una interconexión de componentes que intercambian información, energía y materia, con la finalidad de lograr uno o varios objetivos. Un controlador automático es un sistema que se interconecta con otro denominado planta, para intentar que el conjunto se comporte de acuerdo a ciertas especificaciones establecidas previamente, sin que sea necesario ningún tipo de intervención humana. Los objetivos técnológicos que justifican la inversión económica que implica el proceso de diseño, instalación y mantenimiento de un controlador automático pueden ser múltiples: mejorar el funcionamiento, reducir costos, aumentar la seguridad, disminuir el consumo de energía y/o la emisión de desechos, etc. Prácticamente todos los sistemas necesarios para dar soporte a la vida humana moderna necesitan controladores automáticos para su correcto funcionamiento (fábricas, máquinas, reactores, plantas de tratamiento de agua, vehículos, turbinas, sistemas de bombeo, computadoras, equipos bio-médicos, etc). En algunos casos resulta incluso imposible realizar las operaciones en modo manual (sin la intervención de un controlador automático) debido a la complejidad de los procesos y la necesidad de precisión y rapidez en la toma de decisiones.

1.1. Antecedentes Hasta la segunda mitad del siglo XX, los métodos de análisis y diseño de controladores automáticos se basaban en la interpretación de diagramas y gráficos propuestos por investigadores pioneros en el área tales como Bode, Evans, Nyquist y Nichols. De manera sorprendente por su sencillez, estos métodos permitieron controlar sistemas de bajo orden, de una sola entrada y una sola salida, obteniendo un desempeño aceptable para la época.

2 A partir de 1960, la disponibilidad de equipos de computación cada vez más rápidos y con mayor capacidad de memoria, hizo posible resolver problemas de mayor complejidad, obteniendo al mismo tiempo mejor desempeño de los lazos de control. El diseño de un controlador automático comienza a partir de entonces a formularse formalmente como un problema de optimización con múltiples objetivos y restricciones. Es decir, se considera que el proceso de diseño termina cuando se ha conseguido una solución que optimiza un conjunto de funciones y satisface una serie de restricciones expresadas en forma de desigualdades. En este sentido, uno de los problemas de diseño multi-objetivo más conocidos es el Linear Quadratic Regulator (LQR), el cual consiste en diseñar un controlador que estabilice el sistema a lazo cerrado y que minimice al mismo tiempo dos cantidades: a) la norma cuadrática de la desviación de los estados con respecto a su estado de reposo y b) la norma cuadrática del esfuerzo de control (Kalman, 1960). Considere por ejemplo un sistema lineal descrito por la ecuación de estado •x(t) = Ax(t) + Bu(t) , x(0) = 0

(1.1)

que se supone controlable y observable. El problema LQR consiste en determinar el vector de entrada u(t) en función del vector de estados x(t) de forma que la suma ponderada
∞
∞ T x(t) Qx(t)dt + u(t)T Ru(t)dt (1.2)

0

0

tenga el mínimo valor posible. En esta formulación, las matrices Q y R juegan un papel importante ya que ayudan a ponderar la influencia de cada término en el diseño final del controlador. De hecho, considerar uno solo de estos términos a efectos de la determinación de una estrategia de control, podría conducir a situaciones catastróficas. Si se diseña una solución intentando minimizar la energía del vector de estados x(t), es posible que aumente la energía de la señal de entrada u(t) de manera exagerada, lo cual pudiera producir daños o al menos saturación de los actuadores. Si por el contrario diseñamos el controlador intentando minimizar la energía de u(t), entonces el tiempo de respuesta del sistema pudiera resultar demasiado lento y la acción del controlador sería inocua. Es importante constatar que las dos cantidades que se desea

3 minimizar se encuentran realmente en conflicto puesto que si se logra que una de ellas disminuya, es probable que la otra aumente al mismo tiempo. El ejemplo analizado en el párrafo anterior ilustra de manera sencilla la necesidad de conseguir un compromiso entre diferentes criterios en conflicto, lo cual es un fenómeno ampliamente estudiado en ingeniería de control (rápidez vs precisión, desempeño vs robustez) y en otras áreas tales como la economía (calidad vs costo, costo vs tiempo, etc). De hecho, durante las últimas décadas se ha dedicado un enorme esfuerzo en el desarrollo de métodos que permitan resolver problemas que involucran objetivos en conflicto, en particular aquellos formulados en función de la norma H2 o H∞ de las matrices de transferencia características de un lazo (Apkarian et al, 2008; Boyd et al, 1994; Doyle et al, 1989; Scherer, 1997). En este contexto, una de las estrategias que cuenta actualmente con mayor aceptación consiste en reducir el problema (algunas veces al costo de obtener soluciones conservadoras) a la determinación de la factibilidad de una inecuación matricial lineal (en inglés, LMI). De acuerdo con el investigador Stephen Boyd (1994), la conexión entre el estudio de las LMIs y los problemas de control fue establecida por primera vez aproximadamente en 1892, año en que el matemático ruso Lyapunov demuestra en su tesis doctoral que el sistema descrito por: •x(t) = Ax(t) (1.3) es estable, si y solo si existe una matriz S positiva definida tal que AT S + SA < 0

(1.4)

Una de las principales ventajas de la formulación del problema de diseño en términos de LMIs, es que existen herramientas de cómputo muy eficientes para decidir con respecto a su factibilidad. De hecho, Boyd (1994) es catégorico al afirmar que “gracias a estas herramientas, el problema teórico de control de sistemas lineales, tal como se planteaba durante la etapa clásica, se encuentra definitivamente resuelto”. No obstante, debido a diversas razones técnicas o económicas a menudo es necesario adoptar una estructura pre-establecida para el controlador: Proporcional-IntegralDerivativo (PID), compensador por adelanto/atraso, control decentralizado, etc. También puede ser necesario considerar índices de desempeño de muy variada naturaleza en

4 diferentes canales de la planta (ver capítulo 2). Lo cual finalmente trae como consecuencia que las funciones que intervienen en el problema de optimización sean no-convexas y que la aproximación convexa (de tipo LMI por ejemplo) pueda resultar demasiado conservadora. En estas circunstancias, una alternativa posible consiste en el empleo de algoritmos aleatorios los cuales pueden ser aplicados para resolver el problema de diseño original, sin necesidad de introducir aproximaciones conservadoras. De hecho, estos algoritmos no exigen propiedades matemáticas tales como continuidad o convexidad, ni requieren información sobre el gradiente de las funciones objetivo para su funcionamiento. Se basan exclusivamente en el muestreo de un número determinado de puntos en el espacio de búsqueda, acompañado de algún mecanismo que intenta guiar el proceso de búsqueda hacia las regiones con mayor potencial para producir buenas soluciones. En contrapartida, una importante desventaja de los algoritmos aleatorios es que en general no garantizan convergencia hacia una solución óptima. De hecho, pueden producir resultados diferentes (aunque bastante aproximados entre sí) cada vez que se ejecutan a partir de las mismas condiciones iniciales: incertidumbre que genera serias dudas con respecto a la confiabilidad en un entorno industrial, particularmente cuando se consideran aplicaciones críticas y decisiones que deben tomarse en tiempo real. A pesar de estos cuestionamientos, numerosos investigadores coinciden en que los algoritmos aleatorios pueden representar una alternativa válida para intentar obtener mejores resultados, cuando las soluciones obtenidas previamente por otros métodos no son completamente satisfactorias. Dentro de este contexto, un algoritmo genético se define como un tipo de algoritmo aleatorio que intenta simular el proceso de evolución biológica con el objeto de resolver un problema de optimización. Desde hace más de una década se han utilizado para resolver problemas en el área del control de procesos que involucran múltiples objetivos en contradicción (Fleming, 2001). A continuación se presentan algunas referencias, con el fin de dar una idea del trabajo desarrollado en esta área. Uno de los primeros artículos fue presentado por Fonseca y Fleming (1994). Estos autores describen una aplicación del algoritmo Multi-Objective Genetic Algorithm (MOGA) para el diseño de un regulador lineal de segundo orden, parametrizado por su ganancia,

5 sus polos y sus ceros, de la forma C(s) = K

(s + z1 ) (s + z2 ) (s + p1 ) (s + p2)

(1.5)

La planta a lazo abierto estudiada tiene como función de transferencia −s + 4 G(s) = 2 s (s + 4)

(1.6)

y consideran los siguientes objetivos de diseño: 1. Estabilidad del sistema a lazo cerrado 2. Valor medio cuadrático de la señal de control. 3. Valor medio cuadrático de la señal de salida. 4. Sensibilidad ante perturbaciones. 5. Sensibilidad ante incertidumbres en el modelo. 6. Estabilidad del controlador. Los autores afirman que su método de diseño permite encontrar soluciones muy cercanas a las producidas resolviendo el problema LQR con distintos valores de las matrices R y Q, pero considerando una estructura de control más sencilla. Posteriormente, Fonseca (1995) publica su tesis de doctorado, en la cual presenta una aplicación de MOGA para el control de una turbina de gas caracterizada por un sistema de ecuaciones diferenciales no-lineales. Para este problema se propuso un controlador parametrizado por 8 variables de decisión (ganancias) y 7 objetivos de diseño. Este autor afirma que la presencia de múltiples discontinuidades y amplias mesetas en cada una de las funciones objetivo, hace que los métodos que utilizan información de gradiente presenten dificultades cuando se intentan aplicar a esta clase de problemas y justifica por ende el enfoque genético. Liu et al (1995) resuelven un problema de diseño con dos grados de libertad, tal como es posible apreciar en la figura 1.1. El modelo de la planta a controlar viene dado en tiempo discreto por b+ (z)b− (z) (1.7) G(z) = a+ (z)a− (z)

6

Figura 1.1: Modelo a lazo cerrado en tiempo discreto. Controlador de dos grados de libertad. donde los polinomios b+ (z) y a+ (z) contienen los ceros y polos que el controlador puede cancelar. El problema consiste en diseñar las dos funciones de transferencia g(z)a+ (z) (z − 1)d(z)b+ (z) a 1 a2 C2(z) = a0 + + 2 + ... + z z C1(z) =

(1.8) an zn

(1.9)

donde los parametros ai , i ∈ {1, 2, . . . , n} y los coeficientes de los polinomios g(z) y d(z) son las variables del problema de optimización. El algoritmo genético que estos autores proponen es similar al Strength Pareto Evolutionary Algorithm 2 (SPEA2) propuesto por Zitzler et al (2001). Comparan sus resultados con respecto a los obtenidos mediante un método de optimización de tipo quasi-Newton y en sus conclusiones afirman que logran disminuir los tiempos de subida, establecimiento y sobre-elongación, utilizando para ello mayor energía en la señal de control. En su tesis doctoral, Herreros (2000) presenta un estudio sobre el diseño de controladores robustos multi-objetivo por medio de algoritmos genéticos. Formula un problema de tipo min-max que tiene como variables de decisión los coeficientes característicos de la función de transferencia del controlador y del modelo de incertidumbre. Propone un nuevo algoritmo llamado MRCD (Multi-Objective Robust Control Design) junto con sus respectivos operadores. Para validar sus resultados, los compara con respecto a los generados mediante LMIs. De acuerdo con el autor, el algoritmo genético produce mejores aproximaciones (lo cual se ilustra gráficamente) aunque no presenta pruebas estadísticas que sustenten esta afirmación.

7 En 2002, Liu et al publican una monografía que cubre buena parte de los conceptos básicos en el área de control multi-objetivo. Como idea novedosa proponen una codificación genética de los individuos de la forma p = (λ, R, L)

(1.10)

donde λ es el vector de auto-valores, R es la matriz de auto-vectores derechos y L es la matriz de auto-vectores izquierdos de la matriz de estados a lazo cerrado. Estos autores consideran como objetivos de diseño las sensibilidades individuales de cada auto-valor. En (Takahashi et al., 2004) se propone un algoritmo genético que asegura la generación de un conjunto de soluciones que cumple con las condiciones necesarias para ser un Conjunto de Pareto. Este algoritmo se aplica al problema de diseño mixto H2 /H∞, consiguiendo, de acuerdo con los autores, mejores soluciones que las obtenidas mediante LMIs. Sin embargo estos autores tampoco presentan pruebas estadísticas que sustenten este hecho. En (Tan et al., 2005) se presenta un estudio sobre optimización multi-objetivo evolutiva, con aplicaciones para problemas de ingeniería tales como control de procesos y R programación de horarios. Se presenta además una interfaz de diseño en MATLAB llamada MOEA-Toolbox, la cual permite la formulación y ejecución de los programas de manera amigable (Joos, 2002). En (Pedersen, 2005) se proponen algunas variaciones del algoritmo Non-Dominated Sorting Genetic Algorithm (NSGA-II) para adaptarlo al diseño de controladores de bajo orden en el dominio frecuencial. Se analiza con detalle el problema del escalamiento de las funciones objetivo y su efecto en la evolución del proceso de búsqueda y en la distribución final de las soluciones. En los últimos años se han reportado aplicaciones de diversos conceptos novedosos relacionados con dominancia y evolución. Por ejemplo, Martinez et al (2006) proponen una estrategia de optimización multi-objetivo llamada Programación Física (PF) para resolver problemas de control robusto. Heo et al (2006) utilizan el algoritmo Multiple Objective Particle Swarm Optimization (MOPSO) para diseñar el controlador de una planta generadora de energía eléctrica. En (Gambier et al., 2006) se propone un nuevo método de diseño basado en juegos dinámicos cooperativos. En (Dellino et al., 2007) se propone una metodología híbrida para resolver problemas de alta complejidad en un

8 contexto multi-disciplinario. Finalmente, en (Wang et al., 2009) se presenta un algoritmo que intenta optimizar la posición de los polos del sistema a lazo cerrado utilizando MOPSO.

1.2. Planteamiento del problema de investigación A partir del análisis de las referencias presentadas en la sección anterior, es posible constatar que se ha desarrollado un importante trabajo en el área de la aplicación de Multi-Objective Evolutionary Algorithms (MOEAs) para resolver problemas de control automático. En líneas generales, se ha puesto mucho énfasis en ilustrar las ventajas de esta metodología de forma experimental, mediante casos de estudio y aplicaciones prácticas. El gran desafío en esta área de investigación, tal como lo plantean en su trabajo Fleming y Purshouse (2002), consiste en “convencer a los escépticos de que estos nuevos métodos merecen plena confianza en un entorno industrial ”. No obstante, todavía persisten numerosos aspectos prácticamente inexplorados en esta dirección, algunos de los cuales se introducen brevemente en los siguientes párrafos.

1.2.1

Consideración de las restricciones de diseño.

La restricción de estabilidad, por sólo nombrar un ejemplo, es fundamental para cualquier diseño. Un sistema inestable no debe funcionar en la práctica, puesto que pueden producirse daños o incluso destrucción de sus componentes. El enfoque tradicional para atacar los problemas de control automático mediante algoritmos genéticos consiste en definir una función para penalizar fuertemente a los individuos inestables, de forma que en el transcurso del proceso evolutivo se logre que la población esté compuesta por una mayoría de individuos estables. Como ejemplo, considere la función de penalidad propuesta por Herreros (2000)   0 si max [Re(λ(K))] < 0 f (K) = (1.11)  max [Re(λ(K))] en caso contrario donde K representa la matriz de transferencia del controlador evaluado y λ es el vector de autovalores del sistema a lazo cerrado. Con esta función se penalizan los individuos inestables (no factibles) con una cantidad igual a la mayor abscisa de estabilidad del sistema a lazo cerrado. Aunque esta estrategia tiene la ventaja de ser muy fácil de implementar, en la práctica puede producir fenómenos indeseados tales como convergencia

9 prematura de la población. El reto consiste en explorar métodos alternativos que permitan un manejo más eficiente de las restricciones tales como la generación de individuos estables por construcción y/o la reparación mediante operadores especiales (ver capítulo 4 y 6).

1.2.2

Representación del espacio de búsqueda

Es posible representar matemáticamente un controlador lineal e invariante en el tiempo, de una entrada y una salida, de dos formas: • Mediante una función de transferencia K(s) =

amsm + am−1sm−1 + . . . + a1s + a0 sn + bn−1sn−1 + . . . + b1s + b0

• Mediante un par de ecuaciones matriciales de estado   x˙ (t) = A x (t) + B y(t) K K K K  u(t) = C x (t) + D y(t) K K

(1.12)

(1.13)

K

La ventaja de estas dos representaciones es que son generales y por lo tanto permiten imponer cualquier estructura deseada para el controlador. Su desventaja es que, desde el punto de vista de la búsqueda genética, no es sencillo generar y mantener individuos que satisfagan restricciones de estabilidad. Para ilustrar este fenómeno, considere el modelo de planta propuesto por Fonseca y Fleming (1993): 4−s G(s) = 2 s (4 + s)

(1.14)

Suponga que se desea diseñar un controlador PID ideal mediante un algoritmo genético y se decide adoptar una representación de los individuos que consiste simplemente en un vector que contiene los parámetros KD , KP y KI característicos del controlador (ver figura 1.2). Suponga además que en alguna etapa del proceso de búsqueda se seleccionan los siguientes padres factibles (estabilizantes)     p1 = (1.15) KD1 KP 1 KI 1 = 3.5 1 2     p2 = (1.16) KD2 KP 2 KI 2 = 3.5 2 6

10

Figura 1.2: Sistema a lazo cerrado con controlador PID los cuales generan los siguientes vectores de polos a lazo cerrado 

−0.1781 + 3.5024i

   −0.1781 − 3.5024i λp1 =    −0.0719 + 0.8033i  −0.0719 − 0.8033i 

−0.1972 + 3.0607i    −0.1972 − 3.0607i λp2 =    −0.0528 + 1.5964i  −0.0528 − 1.5964i

               

(1.17)

(1.18)

Mediante el operador de cruce uniforme, estos dos individuos generan el par de hijos   h1 = (1.19) 3.5 2 2   (1.20) h2 = 3.5 1 6

los cuales a su vez generan los siguientes vectores de polos a lazo cerrado 

0.0163 + 3.3630i    0.0163 − 3.3630i λh1 =    −0.2663 + 0.7978i  −0.2663 − 0.7978i

       

(1.21)



−0.4214 + 3.3099i    −0.4214 + 3.3099i λh2 =    0.1714 + 1.4582i  0.1714 − 1.4582i

       

11

(1.22)

Puesto que ninguno de los hijos generados produce un sistema a lazo cerrado estable, el objetivo que se plantea en este caso consiste en desarrollar representaciones alternas del espacio de controladores que permitan mejorar el proceso de generación de individuos factibles y la exploración del espacio de búsqueda mediante operadores de variación genética especialmente adaptados (ver capítulo 6).

1.2.3

Integración con operadores de búsqueda local

Una de las principales desventajas de los algoritmos genéticos es que no existen garantías con respecto a la precisión de las soluciones obtenidas. Considere por ejemplo que se desea resolver el problema de minimización de la función de Rosenbrock (P1) P1 : x1 , x2 ∈ [−10, 10]min100(x2 − x21)2 + (1 − x1 )2

(1.23)

En la tabla 1.1 se muestra un ejemplo de los resultados obtenidos luego de 100 ejecuciones de un algoritmo genético simple (GS) y 100 ejecuciones del algoritmo determinista BFGS, de tipo quasi-Newton (Broyden, 1970). Se presenta el promedio y la desviación estandard del menor valor de la función objetivo encontrado durante la ejecución de cada algoritmo. Para las pruebas se utilizaron poblaciones iniciales aleatorias (en el caso genético) y punto inicial aleatorio (en el caso BFGS). Note que el promedio de la mejor solución obtenida por el algoritmo genético es aproximadamente 1000 veces mayor que el promedio logrado por BFGS. Tabla 1.1: Promedio y desviación standard del valor mínimo alcanzado por cada algoritmo, calculados sobre un total de 100 ejecuciones. Problemas P1 y P2 P1 P2 Genético 5.4271(7.3176) 1.1762(0.6551) BFGS 0.0042(0.0151) 42.0898(41.7667)

12 Considere ahora el siguiente problema (P2)

P2 : x1 , x2 ∈ [−10, 10]min2 + x21 − cos(20πx1 ) + x22 − cos(20πx2)

(1.24)

A partir de la misma configuración utilizada en el caso anterior, se obtienen los resultados mostrados en la tabla 1.1. Note que ahora el promedio BFGS es aproximadamente 35 veces mayor que el logrado por el genético, lo cual se explica considerando que el algoritmo determinista queda “atrapado” en cualquiera de los mínimos locales característicos del problema P2. Los anteriores ejemplos permiten ilustrar un fenómeno sobre el cual existe bastante consenso dentro de la comunidad de investigadores en el área: los algoritmos genéticos facilitan la exploración (búsqueda gruesa), pero presentan serias limitaciones para la explotación de las regiones (búsqueda fina). Surge entonces la idea de combinar la capacidad de exploración del algoritmo genético con la capacidad de explotación de un operador de búsqueda local, para intentar que las zonas más prometedoras sean correctamente aprovechadas. En la literatura, este tipo de algoritmos reciben el nombre de meméticos (Moscato, 1989), haciendo referencia a las unidades de difusión cultural llamadas memes (Dawkins, 1976). En la actualidad son pocas las publicaciones en las que se reportan aplicaciones de algoritmos meméticos para problemas de control automático (Caponio et al, 2007; Shyr et al 2002; Wanner et al, 2008). Un objetivo que se plantea en este trabajo consiste en producir de manera eficiente mejores soluciones considerando un algoritmo memético (ver capítulos 5 y 6) en comparación con el genético puro.

1.3. Estructura de la tesis En los próximos capítulos se presenta un conjunto de operadores, representaciones genéticas y adaptaciones algorítmicas, especialmente concebidos para intentar obtener mejores aproximaciones de los Conjuntos de Pareto, en comparación con otros métodos competitivos citados comúnmente en la literatura. En el capítulo 2 se introducen conceptos básicos relacionados con el control de sistemas lineales a lazo cerrado. Se presenta el problema general de control, haciendo énfasis en su formulación como un problema multi-objetivo. Se describen algunos índices que usualmente se utilizan para especificar los requerimientos que debe cumplir un controlador.

13 También se repasan algunos métodos clásicos de diseño incluyendo PIDs, reubicación de polos y LMIs. En el capítulo 3 se introducen conceptos generales relacionados con la optimización de múltiples objetivos, conjuntos y frentes de Pareto. Se hace un breve recuento de las ideas que pueden utilizarse para resolver problemas de esta naturaleza, tomando en cuenta el archivo de las soluciones no-dominadas y el tratamiento de las restricciones. Seguidamente se describe el funcionamiento de los algoritmos MOGA, NSGA-II y SPEA2, junto con los índices de desempeño y las pruebas estadísticas que se han propuesto para su evaluación. En el capítulo 4 se proponen nuevos operadores de variación para el diseño de controladores PID, basados en la generación y mantenimiento de individuos estables durante del proceso de búsqueda. Se utiliza el método propuesto en (Hohenbichler & Abel, 2008) para generar los vértices de la región factible correspondiente al sistema a lazo cerrado. Adicionalmente se propone un nuevo método para aproximar la geometría de la región factible, basado en muestreo aleatorio uniforme (Calafiore & Dabbene, 2006). Los resultados de este capítulo fueron previamente publicados en (Sánchez et al., 2008a) y (Sánchez et al., 2009). En el capítulo 5 se presenta un nuevo operador de búsqueda local llamado Ascenso de Colina con Paso Lateral (ACPL) especialmente diseñado para problemas multi-objetivo. En caso de no encontrar soluciones dominantes (como en el ascenso de colina tradicional), el operador utiliza la información recolectada para intentar moverse hacia algún cono de diversidad (paso lateral). También se presenta un análisis de la sensibilidad del operador con respecto a sus parámetros característicos. Por último se propone un esquema para la integración del ACPL con el algoritmo SPEA2, dando como resultado un algoritmo memético denominado SPEA2-ACPL, para el cual se presentan los resultados obtenidos durante las pruebas que se realizaron en el caso de un problema convexo y uno no-convexo. Estos resultados fueron previamente publicados en (Schütze et al., 2008) y (Lara et al., 2010). En el capítulo 6 se propone un nuevo método de diseño de controladores llamado MultiObjective Pole Placement with Evolutionary Algorithms (MOPPEA), utilizando como base del proceso de búsqueda a los algoritmos SPEA2 y SPEA2-ACPL. Se presentan los resultados obtenidos para el diseño mixto H2 /H∞ en el caso de cinco modelos extraidos de COMP leib (Leibfritz, 2004) y se comparan con respecto a los obtenidos mediante LMIs.

14 Los resultados del capítulo fueron publicados en (Sánchez et al., 2007) y (Sánchez et al., 2008b). Por último, se presentan las conclusiones del trabajo y se describen algunas direcciones para futuras investigaciones.

CAPÍTULO II SISTEMAS LINEALES REALIMENTADOS

En este capítulo se introducen algunos conceptos básicos relacionados con el control de sistemas lineales a lazo cerrado. Se presenta el problema general de control, haciendo énfasis en su formulación como un problema multi-objetivo. Se definen los índices que usualmente se incluyen como parte de las especificaciones que debe cumplir un controlador automático. Por último se realiza un breve repaso de tres estrategias bien conocidas que se han sido propuestas para resolver problemas de diseño de controladores lineales: sintonización PID, reubicación de polos mediante realimentación de estados y LMIs.

2.1. Generalidades Un sistema lineal e invariante (ver figura 2.1) puede definirse mediante las ecuaciones de estado    x(t) ˙ = Ax(t) + Bw w(t) + Bu u(t)   (2.1) z(t) = Cz x(t) + Dzw w(t) + Dzu u(t)     y(t) = C x(t) + D w(t) + D u(t) y

yw

yu

donde u ∈ Ln2 u [0, +∞) es el vector de señales manipuladas, w ∈ Ln2 w [0, +∞) es el vector n de entradas exógenas (perturbaciones), y ∈ L2 y [0, +∞) es el vector de señales medidas, z ∈ Ln2 z [0, +∞) es el vector de señales controladas y x ∈ Ln2 [0, +∞) es el vector de estado (el espacio L2[0, +∞) es definido en el Apéndice).

Figura 2.1: Sistema LTI a lazo abierto

16

Las matrices A ∈ Rn×n , Bw ∈ Rn×nw , Bu ∈ Rn×nu , Cz ∈ Rnz ×n , Cy ∈ Rny ×n , Dzw ∈ Rnz ×nw , Dzu ∈ Rnz ×nu , Dyw ∈ Rny ×nw y Dyu ∈ Rny ×nu se denominan matrices de estado del sistema a lazo abierto. Un sistema se dice estable en el sentido BIBO (Bounded Input Bounded Output) si a cualquier entrada acotada, corresponde una salida acotada. Es posible demostrar que el sistema (2.1) es estable en el sentido BIBO si y sólo si los autovalores de la matriz A pertenecen a C− (semi-plano abierto complejo izquierdo). En este caso se dice que la matriz de estado A es de tipo Hurwitz. Los conceptos de controlabilidad y observabilidad (Kalman, 1963) juegan un papel importante para el diseño de controladores. Se dice que un sistema es controlable si es posible transferir un estado inicial a cualquier otro estado finito en un intervalo de tiempo finito. Si todos los estados son controlables, se dice que el sistema es completamente controlable. Por otro lado, un sistema es completamente observable si es posible reconstruir la información de cada variable de estado en función de las mediciones de las salidas y las entradas. De esta forma, el sistema (2.1) se dice completamente controlable y observable si y sólo si las matrices   n − 1 (2.2) Bu ABu · · · A Bu son de rango completo.



Cy Cy A · · · Cy An−1

T

(2.3)

La Transformada de Laplace (también definida en el Apéndice) de los vectores y(t) y z(t) se escribe Z(s) = Gzw (s)W (s) + Gzu (s)U (s)

(2.4)

Y (s) = Gyw (s)W (s) + Gyu (s)U (s) con Gzw (s) = Cz (sI − A)−1Bw + Dzw G (s) = C (sI − A)−1B + D zu

z

u

zu

(2.5)

17

Gyw (s) = Cy (sI − A)−1 Bw + Dyw Gyu (s) = Cy (sI − A)−1 Bu + Dyu

donde Gzw (s), Gzu (s), Gyw (s) y Gyu (s) se denominan matrices de transferencia a lazo abierto. La necesidad de interconectar el sistema a lazo abierto con un controlador automático surge de los siguientes hechos: • Existen perturbaciones (cuya acción es imposible evitar) que desplazan al sistema del régimen de operación deseado. • Es imposible conocer con exactitud la estructura y/o los parámetros del sistema real (incertidumbre). • Ante las perturbaciones y en presencia de incertidumbre, el sistema debe reaccionar de una manera que sería imposible lograr mediante un sistema a lazo abierto (sin realimentación). Cuando se interconecta el sistema (2.1) con un controlador lineal (ver figura 2.2) representado por   x˙ = A x + B y K K K K (2.6)  u=C x K K

el sistema a lazo cerrado resultante se escribe   x˙ = A x + B w cl cl cl cl  z=C x +D w cl cl

(2.7)

cl

donde



Acl = 

A

Bu CK



,

BK Cy AK + BK Dyu CK   Ccl = Cz Dzu CK , Dcl = Dzw



Bcl = 

Bw BK Dyw

 

(2.8)

Por su parte la matriz de transferencia a lazo cerrado puede escribirse como Gzwcl (s) = Gzw (s) + Gzu (s)K(s)(I − Gyu (s)K(s))−1 Gyw (s)

(2.9)

18

Figura 2.2: Interconexión Planta/Controlador (lazo cerrado) con

K(s) = CK (sI − AK )−1BK

(2.10)

Es evidente, a partir de estas ecuaciones, que no es trivial la selección de matrices AK , BK y CK tales que el sistema a lazo cerrado sea estable y cumpla con una lista predeterminada de requerimientos.

2.2. Índices de desempeño A menudo se requiere que las salidas del sistema a lazo cerrado sigan cierta dinámica en respuesta a las perturbaciones y en presencia de incertidumbre. En otras palabras, se requiere dotar a la conexión (figura 2.2) de propiedades que el sistema a lazo abierto (figura 2.1) no posee de manera aislada. En el transcurso de las últimas décadas, se han propuesto numerosos criterios o índices para medir el desempeño de un controlador automático, es decir, para cuantificar que tan bien cumple con los objetivos para los cuales fue diseñado. A continuación se definen algunos de los más comunes.

2.2.1

Índices de estabilidad

Sea ncl el grado y λ el vector de polos del sistema a lazo cerrado. Es posible considerar como índice de estabilidad el número de polos estables ns = |{λi / Re(λi ) < 0}|

(2.11)

19 o a la abscisa del polo que se encuentre más a la derecha en el plano complejo, sea qa tal que (2.12) qa = max Re(λi ) i

En las figuras 2.3 y 2.4 se muestra un ejemplo del cálculo de estos índices para dos configuraciones de polos distintas. Note que de acuerdo al índice ns el primer sistema sería mejor que el segundo. Con respecto al índice qa la conclusión sería exactamente la opuesta.

Figura 2.3: Diagrama de polos de un sistema con índices de estabilidad ns = 2, qa = 3

Figura 2.4: Diagrama de polos de un sistema con índices de estabilidad ns = 1, qa = 2

2.2.2

Índices de respuesta temporal

Cuando el sistema a lazo cerrado es estable, es posible definir índices en función de su respuesta temporal ante señales de prueba en distintos canales de entrada. Considere el caso de un sistema de tipo Single-Input Single-Output (SISO).

20 Para una entrada de referencia r(t) en forma de escalón, es decir   0 para t < 0 r(t) =  R para t ≥ 0

(2.13)

el error en estado estacionario ess se define como:

(2.14)

ess = lim |r(t) − y(t)| t→∞

donde y(t) es la respuesta del sistema. Note que de manera similar es posible definir el error en estado estacionario para entradas más complejas tales como rampas, senoidales, triangulares, etc. El tiempo de subida tr al 90%, el tiempo de asentamiento ts al 5% y el porcentaje de sobre-elongación máxima Mp pueden definirse como se muestra a continuación (ver figura 2.5). Mp(%)

y(t)

1.05 1 0.95 0.90

ess

t tr

ts

Figura 2.5: Indices de desempeño a partir de la respuesta ante un escalón unitario   tr (90%) =arg min t : y(t) = 0.9× lim y(t) t→∞ t      ts (5%) =arg max y(t)− lim y(t) ≥ 0.05 t→∞ t

 y(t) − lim y(t) ∞ t→∞ × 100 si y(t) > lim y(t)  ∞ t→∞  y (t)∞ Mp =  0 en caso contrario

(2.15)

(2.16)

(2.17)

21 Uno de los conflictos más comunes cuando se analiza la respuesta de un sistema a lazo cerrado es el que existe entre el tiempo de subida tr y el porcentaje de sobre-elongación máxima Mp. En efecto, si se diseña el controlador para obtener el menor tiempo de subida, es probable que se obtenga al mismo tiempo una mayor sobre-elongación. Un problema de diseño multi-objetivo que hiciera intervenir los cuatro índices que se acaban de definir más la restricción de estabilidad pudiera plantearse como   min ess(K) tr (K) ts (K) Mp(K) K ∈R(s)nu ×ny (2.18) sujeto a: qa(K) < 0

2.2.3

Rechazo de perturbaciones

Sea w(t) una perturbación cuyo efecto sobre la salida z(t) es necesario rechazar y Gzwcl (s) la matriz de transferencia definida en (2.9). El rechazo de la señal de perturbación puede medirse en función de una norma (ver Apéndice A) tal como W (s)Gzwcl (s, K)i

(2.19)

donde i ∈ {1, 2, ∞} y la matriz W (s) se incluye de manera arbitraria para ponderar la importancia de las bandas de frecuencia en las cuales w(t) concentra su mayor energía. Considere por ejemplo una perturbación en forma de impulso unitario w(t) = δ(t). Para minimizar la energía de la señal de salida, se debe plantear un problema de la forma min

K ∈R(s)nu ×ny

W (s)Gzwcl (s, K)2 (2.20)

sujeto a qa (K) < 0

Si en lugar de minimizar la energía de la señal de salida se necesita minimizar su máxima amplitud (por ejemplo para prevenir posibles situaciones peligrosas), se debe entonces plantear el problema min

K ∈R(s)nu ×ny

W (s)Gzwcl (s, K)∞

sujeto a qa(K) < 0

(2.21)

22

2.2.4

Estabilidad robusta

Ningún modelo matemático puede predecir con absoluta certeza la respuesta de un sistema en condiciones reales de funcionamiento. En consecuencia, cualquier método de diseño de controladores debe considerar este hecho, introduciendo incertidumbre en el modelo. Considere por ejemplo el lazo mostrado en la figura 2.6. La matriz de transferencia ∆(s) representa una incertidumbre no estructurada que se adiciona al modelo de la planta a lazo abierto G(s). Suponga que ∆(s) esté acotada de forma que se cumpla ∆(s)∞ ≤ 1

(2.22)

Entonces es posible demostrar que el sistema a lazo cerrado es estable si y sólo si se cumple   K(s) (I + G(s)K(s))−1  < 1 (2.23) ∞

Figura 2.6: Lazo de realimentación con incertidumbre no estructurada ∆ De manera general, se considera que mientras menor sea la norma que aparece en lado izquierdo de la desigualdad (2.23) más robusto será el lazo de control con respecto las incertidumbres sobre el modelo de la planta. Sin embargo, a medida que se mejora la robustez del diseño, es posible que se degrade al mismo tiempo el desempeño del lazo medido en términos de la norma de otro canal del modelo. De esta forma se constata nuevamente la necesidad de adoptar un enfoque multi-objetivo en el tratamiento de los problemas de diseño.

2.3. Controladores PID De acuerdo con Åström y Hägglund (1995) cerca del 95% de los controladores instalados en la industria son del tipo Proporcional-Integral-Derivativo (PID). Su principal

23 ventaja en la práctica consiste en que vienen pre-programados en los equipos electrónicos comerciales y por ende son bien conocidos tanto por los ingenieros como por los operadores y responsables del mantenimiento de las plantas. Desde un punto de vista teórico, un PID asegura un error en estado estacionario nulo para una referencia en forma de escalón gracias a su efecto integrador. Por otro lado, la acción derivativa interviene cuando hay cambios bruscos en el error (si el error cambia lentamente, la contribución principal se debe a los modos proporcional e integral). Considerando las ventajas expuestas anteriormente, se entiende la necesidad de verificar si es posible resolver un problema de diseño mediante una estructura PID (o alguna adaptación de la misma) antes de considerar soluciones más complejas. En la literatura es posible encontrar varias formas para expresar la función de transferencia de un controlador PID. La más simple es la conocida como forma ideal KPID (s) =

KD s2 + KP s + KI s

(2.24)

caracterizada por tres parámetros reales KD , KP , KI que deben ser ajustados por el diseñador del sistema de control (ver figura 2.7).

Figura 2.7: Lazo de control con PID ideal

2.3.1

Método de Ziegler y Nichols

El investigador Aidan O’Dwyer (2006) ha contabilizado 408 reglas para el ajuste de los parámetros PID, de las cuales 293 fueron propuestas después de 1992. Una de las más importantes históricamente fue publicada por Ziegler y Nichols en 1942 y consiste básicamente en las siguientes fórmulas KP = 0.6Ku Ku KI = 1.2 Tu KD = 0.075Ku Tu

(2.25)

24 donde Tu es el período de las oscilaciones que se obtienen a lazo cerrado cuando se incrementa progresivamente la ganancia de la cadena directa desde cero hasta el valor Ku . Considere el siguiente modelo a lazo abierto G(s) =

−0.5s4 − 7s3 − 2s + 1 (s + 1)(s + 2)(s + 3)(s + 4)(s2 + s + 1)

(2.26)

Para este modelo, las oscilaciones críticas se obtienen para Ku = 6.2 y el período de las oscilaciones es Tu = 4.2 lo cual, según las fórmulas de Ziegler y Nichols corresponde a KP = 3.7, KI = 1.8 y KD = 1.95. La respuesta del sistema a lazo cerrado ante una entrada de referencia en forma de escalón se muestra en la figura 2.8. Para este controlador, los índices de respuesta temporal tienen los siguientes valores ess = 0, tr (90%) = 24.5, ts (5%) = 30.9 y Mp = 0. 1

y(t)

0.5

0

−0.5

0

20

40

60

80

100

t

Figura 2.8: Respuesta del sistema con los parámetros KD , KP , KI de acuerdo con las fórmulas de Ziegler y Nichols Aunque goza de mucha aceptación debido a su sencillez, el ámbito de aplicación del método de Ziegler y Nichols, y de la mayoría de los métodos similares, es bastante limitado. Por ejemplo, ninguno garantiza el logro de un compromiso satisfactorio entre un conjunto de índices de desempeño, en relación con los problemas de diseño multi-objetivo formulados en la sección anterior.

2.3.2

Método de Hohenbichler y Abel

Cuando se decide adoptar una estrategia de control PID, surge la duda acerca de cuáles son los valores de los parámetros KP , KD , KI (si es que éstos existen) que garantizan la estabilidad del sistema a lazo cerrado.

25 En este sentido, Hohenbichler y Abel (2008) propusieron un método para calcular los límites de la región de controladores PID estabilizantes en el espacio (KP , KD , KI ). Dado un modelo de planta G(s) a lazo abierto, estos autores demostraron que es posible determinar los intervalos de la ganancia proporcional KP tales que exista al menos un plano (KD , KI ) estabilizante. La única limitación para este resultado es que el sistema a lazo abierto no posea ceros en el eje imaginario. La idea básica de dicho método se fundamenta en el teorema propuesto en (Ackermann & Kaesbauer, 2003), el cual afirma que las regiones estabilizantes en el plano KI −KD para un KP fijo, están formadas por polígonos convexos cuyos vértices pueden ser determinados. Considere una función de transferencia a lazo abierto expresada como: G(s) =

A(s) , B(s) = sR(s) R(s)

(2.27)

con A(s) = am sm + am−1sm−1 + . . . + a1s + a0, am = 0 B(s) = bn sn + bn−1 sn−1 + . . . + b1 s, bn = 0

(2.28)

y para un valor s = jω, denote: RA (ω) = Re [A(jω)] , RB (ω) = Re [B(jω)]

(2.29)

IA (ω) = Im [A(jω)] , IB (ω) = Im [B(jω)] Suponga que se interconecta G(s) con un controlador PID ideal de manera que la ecuación característica del lazo cerrado se escriba como (ver figura 2.7):   P (s) = KD s2 + KP s + KI A(s) + B(s) = 0

(2.30)

Ackermann y Kaesbauer (2003) demostraron el siguiente resultado: Para un valor fijo de KP , la región de estabilidad del polinomio P (s) en el plano KI − KD está compuesta por polígonos convexos, cuyos vertices se definen mediante la intersección de las siguientes rectas: (2.31) KI = 0    0 si n < m + 2   KD = (2.32) − abmn , si n = m + 2     ∞, si n > m + 2

KI = KD ω 2 −

26 RA (ω)RB (ω) + IA (ω)IB (ω) 2 (ω) + I 2 (ω) RA A

(2.33)

donde ω > 0 (denominada frecuencia singular) es solución real de la ecuación ω=

IA (ω)RB (ω) − RA (ω)IB (ω) 2 (ω) + I 2 (ω)) KP (RA A

(2.34)

Considere la función de transferencia propuesta anteriormente en (2.26). Para KP = −4.4 las rectas que delimitan la región de estabilidad en el plano KI − KD , junto con sus respectivas frecuencias singulares, se muestran en la tabla 2.1. En la figura 2.9 se grafican dichas rectas, de manera que es posible observar los límites de la región de estabilidad para este valor particular de KP . La respuesta al escalón del sistema a lazo cerrado correspondiente a los valores KP = −4.4, KI = 2.2 y KD = −10 se muestra en la figura 2.10. En este caso los índices de respuesta temporal tienen los siguientes valores ess = 0, tr (90%) = 15.9, ts(5%) = 32.6 y Mp = 17%. Note que a diferencia del lazo diseñado mediante Ziegler-Nichols, esta respuesta presenta mayor sobre-elongación pero en contrapartida presenta un menor tiempo de subida. Queda claro entonces que el valor de los índices de desempeño del sistema a lazo cerrado dependen de la posición del vector de parámetros KP , KI , KD en el espacio de controladores estabilizantes.

Tabla 2.1: Frecuencias singulares y ecuaciones de las rectas que delimitan la región de estabilidad para KP = −4.4, modelo (2.26) ω KI − KD KI = 0 KD = ∞ 0.32 KI = 0.1KD + 6.8 3.65 KI = 24KD − 88.8

27 8 6

KI

4 2 0 −2 −4 −70

−60

−50

−40

−30 K

−20

−10

0

D

Figura 2.9: Límites de la región de estabilidad para KP = −4.4 para el modelo (2.26)

2.4. Reubicación de polos En esta sección suponga que se desea diseñar un controlador para el sistema lineal   x(t) ˙ = Ax(t) + Bu(t) (2.35)  y(t) = Cx(t) + Du(t)

que se supone controlable y observable. Considere el cambio de variables de estado definido mediante x = T x, u = Ru + F x (2.36) donde T y R son matrices invertibles. Se obtiene de esta forma una representación equivalente del sistema a lazo abierto:  ·  x= (T −1AT + T −1BF )x + (T −1 BR)u (2.37)  y = (CT + DF )x + DRu Por otra parte, si el vector de estados se encuentra disponible y se realimenta mediante una matriz estática K ∈ Rnu ×n de la forma u = Kx, u = R−1(KT − F )x = Kx tendremos dos representaciones equivalentes a lazo cerrado  ·   x= (A + BK)x  x˙ = (A + BK)x ⇔  y = (C + DK)x  y = (C + DK)x

(2.38)

(2.39)

28 1.2 1 0.8 0.6

y(t)

0.4 0.2 0 −0.2 −0.4 −0.6 10

20

30

40

50 t

60

70

80

90

100

Figura 2.10: Respuesta al escalón del sistema a lazo cerrado para el modelo (2.26) correspondiente al controlador PID caracterizado por KP = −4.4, KI = 2.2 y KD = −10 con

A = T −1 AT + T −1 BF, B = T −1 BR

(2.40)

C = (CT + DF ), D = DR

(2.41)

A + BK = T −1 (A + BK)T

(2.42)

y además

En el caso de sistemas de una sola entrada y una sola salida, la transformación definida por T =



V

V A . . . V An−1

donde V es un vector tal que

−1

, R = 1, F = −V An T

V Ak−1B = 0, ∀k < n − 1, V An−1 B = 1   produce el par A, B en la forma  0 1    0 0   . . A =  .. ..    0 0  0 0

(2.43)

(2.44)

controlable 0 0

... 0



  1 0 ... 0   .. . .  . ... 0 , .   0 0 1 0   0 0 0 0



     B=    

0



        0   1 0 .. .

(2.45)

29 En particular, si se desea que el sistema a lazo cerrado esté caracterizado por el polinomio p(s) = sn + pn−1sn−1 + pn−2 sn−2 + . . . + p1 s + p0 (2.46) basta con seleccionar K=− de donde es posible deducir



p0 p1 . . . pn−1

K = (RK + F )T −1



(2.47) (2.48)

Si el vector de estados no se encuentra disponible, es posible realimentar un observador de la forma · x = A x + Bu + L(C x  − y) (2.49)

donde (C x − y) es el término de corrección y L ∈ Rn×ny es la matriz de ganancias del observador, la cual debe ser calculada de forma que los autovalores de la matriz A + LC pertenezcan todos a C− . La representación completa en espacio de estados del controlador diseñado se expresa como  ·  x = (A + BK + LC) x − Ly (2.50)  u = K x Considere el sistema a lazo abierto definido por la función de transferencia (2.26). Para asegurar un error en estado estacionario nulo ante una entrada de referencia en forma de escalón unitario, se incluye un polo en el origen en el modelo a lazo abierto. El diagrama de polos y ceros correspondiente se muestra en la figura 2.11. Suponga que se desea aplicar el método explicado en esta sección de manera que a lazo cerrado se tenga el siguiente vector de polos   (2.51) λK = −1 −2 −3 −4 −5 −0.5 + 0.866i −0.5 − 0.866i

Es posible verificar que este objetivo se logra mediante el vector de realimentación de estados   K = − 10 13.75 14.37 14.84 17.03 11.56 15

Para completar el diseño es necesario calcular el vector de ganancias del observador L. Suponga que se desea obtener el vector de polos para el observador λL = λK , lo cual se logra mediante el vector de ganancias  T (2.52) L = − 0 0 0 0 0 0 40

30 La respuesta al escalón del sistema a lazo cerrado se muestra en la figura 2.12. En este caso los índices de respuesta temporal tienen los siguientes valores ess = 0, tr (90%) = 2.85, ts (5%) = 10.65 y Mp = 13%, los cuales son bastante mejores que los obtenidos mediante los controladores PID presentados en la sección anterior. En contrapartida, la complejidad de este controlador es mucho mayor que la de un PID.

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −14

−12

−10

−8

−6

−4

−2

0

Figura 2.11: Diagrama de polos (×) y ceros (◦) a lazo abierto para el modelo (2.26) 1.5

1

y(t)

0.5

0

−0.5

−1

−1.5

0

5

10

15

t

Figura 2.12: Respuesta al escalón del sistema a lazo cerrado, para el modelo (2.26). Controlador diseñado mediante reubicación de polos.

2.5. Diseño de controladores mediante LMIs Una LMI es una inecuación matricial lineal de la forma: m  xi Fi > 0 F (x) = F0 + i=1

(2.53)

31

donde cada xi ∈ R es una variable escalar real y cada Fi ∈ Rn×n es una matriz simétrica

constante. El símbolo de desigualdad en la ecuación (2.53) significa que la matriz F debe ser positiva definida. En numerosos problemas de diseño de controladores, las variables de decisión que deben ser ajustadas son matrices simétricas. Considere por ejemplo la inecuación de Lyapunov AT X + XA < 0 donde



X =

x1 x2





(2.54) 

a a  , A =  11 12  x2 x3 a21 a22

(2.55)

la cual puede reescribirse como en (2.53) definiendo F0 = 0 y −x1F1 − x2 F2 − x3 F3 > 0   2a11 a12  F1 =  a12 0   2a21 a11 + a22  F2 =  a11 + a22 2a12   0 a21  F3 =  a21 2a22

(2.56) (2.57)

(2.58)

(2.59)

De esta manera, resolver la inecuación (2.54) equivale a determinar una matriz X tal que la matriz AT X + XA sea negativa semi-definida. En las próximas secciones se formulan los problemas de diseño H2 y H∞ , de amplia importancia en el área de sistemas automáticos, los cuales pueden ser resueltos determinando la factibilidad de una LMI.

2.5.1

Problema H2 En esta sección se considera el siguiente    · x A Bj        zj  =  Cj Dj    y C Fj

modelo lineal a lazo abierto:   B x     E j   wj    0 u

(2.60)

32 y un controlador lineal K representado mediante:      · xK xc AK BK   =  u CK DK y

(2.61)

de forma que a lazo cerrado se tiene



zj = Tj wj = 

con

Aclj Bclj Cclj Dclj





 wj

(2.62)

 A + BDK C BCK  Acl =  AK BK C   Bj + DK Fj  Bclj =  BK F j   Cclj = Cj + Ej DK C E j CK

(2.63)

(2.64) (2.65) (2.66)

Dclj = Dj + Ej DK Fj

El problema H2 consiste en diseñar un controlador tal que (2.67)

Tj (s)2 ≤ γ 2j

donde γ 2j es el valor deseado para la norma de la matriz de transferencia a lazo cerrado, el cual debe ser fijado por el diseñador. Es posible demostrar (Scherer et al., 1997) que la desigualdad (2.67) se cumple si y  B,  C,  D  tales que: sólo si existen matrices X > 0, Y > 0, A, 



 + BC  AX + XAT + B C   T   A T + A + B DC    

 j Bj + B DF

T



T

  Bj + B DFj   T   + BC   j 0 

(2.69)

33

traza(Q) < γ 2j

(2.70)

 j=0 Dj + Ej DF

(2.71)

Si estas inecuaciones se cumplen, las matrices del controlador correspondiente se calculan mediante las fórmulas: MN T = I − XY  DK = D

(2.72)

−T  − DCX)M  CK = ( C  − Y B D)  BK = N −1(B

AK = N −1(A − N BcCX − Y BCc M T − Y (A + BDc C)X)M −T

Considere el siguiente modelo a lazo abierto, propuesto en (Scherer et al., 1997):   1 0 0 10 2       0 1   −1 1 0   A B1 B     1 0 0 2 −5       (2.73)  C1 D1 E1  =      0 1 0 0 0    C F1 0    0 0 1 0 0    2 0 0 1 0

Para verificar si el anterior conjunto de inecuaciones es o no factible es posible utilizar R . En este caso, fijando γ = 10 las LMIs resultan la instrucción feasp de MATLAB 21 factibles y se obtienen los siguientes resultados   473.7130 −24.7516 −6.9934     (2.74) X =  −24.7516 3.6762 −0.0953    −6.9934 −0.0953 4.2502 



0.9540 1.5255 −4.2185     Y =  1.5255 11.8585 −20.8031    −4.2185 −20.8031 115.0885

(2.75)

34 T1(s)2 = 3.1748 ≤ γ 21

(2.76)

y la función de transferencia del controlador correspondiente es K2 (s) = 1206.314

2.5.2

(s + 5.008)(s + 1.365) (s + 460.4)(s + 5.16)(s + 0.6827)

(2.77)

Problema H∞

De manera análoga a lo planteado en la sección anterior, el problema H∞ consiste en diseñar K tal que (2.78) Tj (s)∞ ≤ γ ∞j  B,  C,  D  tales que: lo cual se cumple si y sólo si existen matrices X, Y, A, 



 + BC  AX + XAT + B C

  T  T  (A + A + B DC)    j )T  (Bj + B DF   Cj X + Ej C

T

∗

∗  T  + BC  AY + Y AT + B C ∗  j )T (Y Bj + BF

 Cj + Ej DC

Con γ ∞1 = 10, las LMIs resultan factibles y se obtienen los  391.3011 −36.2150 −8.0631   X =  −36.2150 56.4724 −4.3847  −8.0631 −4.3847 48.2515

−γ ∞j I

∗

    0, el número de variables de decisión (Roy, 1971). En general, la reducción del valor de una de las funciones objetivo conduce al incremento de otra, por lo cual no existe un x∗ ∈ X que sea mínimo de todas simultaneamente. Para Osyczka (1985) el símbolo “min” debe ser interpretado en el caso multi-objetivo como la búsqueda de un conjunto de soluciones que produzcan valores aceptables para todos los objetivos. Muchos métodos que se han propuesto para resolver este problema se basan en una transformación escalar del vector objetivo, de la forma U : RM → R

(3.2)

37 con el fin de replantear el problema de optimización en forma escalar, como min

z1 ,z2 ,...,zM

sujeto a :

U (z1, z2 , . . . , zM )

(3.3)

zi = fi (x) , 1 ≤ i ≤ M , x ∈ X

Kalianmoy Deb (2001) considera que reducir el problema vectorial a uno escalar, implica en general alguna pérdida en la calidad de las soluciones que pueden obtenerse. Los algoritmos que producen una sola solución deben ejecutarse en repetidas ocasiones para identificar los mejores compromisos, utilizando conocimiento previo sobre el problema para definir prioridades entre los objetivos. Un enfoque alternativo consiste en la búsqueda de soluciones especiales llamadas nodominadas, de acuerdo con los conceptos que se definen a continuación . Dado el problema planteado en (3.1), se denomina Espacio Objetivo y se denota F al conjunto definido por  T (3.4) z ∈F ⇐⇒ ∃x ∈ X : z = F(x) = f1 (x) f2(x) ... fM (x)

Sean z1 y z2 dos vectores pertenecientes a F. Se dice que z2 es dominado por z1 y se denota z1 z2 si y solo si z1i ≤ z2i , i = 1, 2, ..., M

∃i ∈ {1, .., M} : z1i < z2i

(3.5)

Si todas las componentes de z1 son estrictamente menores que las de z2 , es decir, si se cumple (3.6) z1i < z2i , i = 1, 2, ..., M

se dice entonces que z2 es estrictamente dominado por z1 y se denota z1 ≺ z2. También puede suceder que los vectores z1 y z2 sean no dominados mutuamente, es decir, tales que z1  z2 ∧ z2  z1

(3.7)

En este caso se dice que z1 y z2 son incomparables y se denota z1  z2 Sea BM el conjunto de números binarios de M dígitos tales que   b = bM bM −1 . . . b2 b1, bi ∈ {0, 1}, i = 1, 2, . . . , M b ∈ BM ⇐⇒  b = 00 . . . 000, b = 11 . . . 111

(3.8)

(3.9)

38 Se denomina cono de diversidad de tipo b definido en el punto x al conjunto C(x, b) definido por   fi (y) ≤ fi (x) si bi = 0 y ∈ C(x, b) ⇐⇒ para i = 1, 2, . . . , M (3.10)  f (y) > f (x) si b = 1 i i i

Note que el número de conos de diversidad es igual a 2M − 2. Como ejemplo, x1 ∈ C(x0 , 110) si y sólo si f1(x1) ≤ f1(x0) f2(x1) > f2(x0)

(3.11)

f3(x1) > f3(x0) Se denomina Conjunto Óptimo Global de Pareto (COGP) y se denota P ∗ al conjunto para el cual se cumple x∗ ∈ P ∗ ⇔

∄x ∈ X , F(x) F(x∗)

(3.12)

De manera similar, se denomina Conjunto Óptimo Local de Pareto (COLP) y se denota PL∗ a cualquier conjunto para el cual se cumple x∗ ∈ PL∗ ⇔

∄x ∈ XL , F(x) F(x∗)

(3.13)

donde XL es un sub-conjunto de X . Se denomina Frente Óptimo Global de Pareto (FOGP) y se denota PF ∗ al conjunto para el cual se cumple z∗ ∈ PF ∗ ⇔

∃x∗ ∈ X , z∗ = F(x∗)

(3.14)

Se denomina Frente Óptimo Local de Pareto (FOLP) y se denota PF ∗L al conjunto para el cual se cumple (3.15) z∗ ∈ PF ∗L ⇔ ∃x∗ ∈ XL , z∗ = F(x∗) Note que si todas las funciones objetivo son continuas y derivables, una condición necesaria para x∗ ∈ P ∗ es que se cumpla ∇fi (x∗ )T d ≥ 0, ∀i ∈ {1, .., M}, ∀d ∈ RM

(3.16)

De hecho, es posible calcular una dirección de descenso para todas las funciones objetivo (3.1) a partir de la información sobre los gradientes ∇fi (x), i ∈ {1, .., M}.

En efecto, sea x ∈ X y sea d ∈ Rn , d = 0, definida mediante d=−

M  i=1

α∗i ∇fi (x)

donde α∗ ∈ RM es solución del problema  2   min  M α ∇f (x)  i i i =1 M 2 α∈R M  sujeto a αi = 1, αi ≥ 0

39

(3.17)

(3.18)

i=1

entonces Schäffler et al (2002) demostraron que si d es distinto de 0, entonces d es una dirección de descenso para todas las funciones f1, f2, . . . , fM

3.2. Estrategias para resolver un problema multi-objetivo Una gran cantidad de estrategias se han propuesto para resolver problemas que involucran objetivos en conflicto. Las preferencias del diseñador pueden ser consideradas antes, durante o después del proceso de búsqueda. De esta forma es posible establecer la siguiente clasificación.

3.2.1

Formulación de preferencias “a priori”

El rasgo común a este tipo de métodos es que el diseñador define sus preferencias antes de ejecutar el algoritmo, ya sea mediante coeficientes de ponderación, metas mínimas a alcanzar, orden lexicográfico, etc. Método min-max Consiste en plantear el problema de la forma:   c (f (x) − γ 1)  1 1     c2 (f2 (x) − γ 2)      max min  ∗  x∈X i∈{1,2,...,M }       ∗   cM (fM (x) − γ M )

(3.19)

En esta formulación los ci , γ i se eligen para normalizar los valores de las funciones objetivo y en general se utilizan las fórmulas:

40 1 fi max − fi min γ i = fi min ci =

(3.20)

donde fi max =

max fi (x), i ∈ {1, 2, ..., M} (3.21) x∈X fi min = min (x), i ∈ {1, 2, ..., M} x∈X Note que de igual forma podría plantearse la minimización de cualquier otra norma para medir la magnitud del vector de funciones objetivo. Método de suma ponderada El problema se plantea de la forma min x∈X

M  i=1

wi ci (fi (x) − γ i )

(3.22)

donde los wi son coeficientes que se eligen para ponderar la importancia relativa de cada función objetivo. La desventaja de este enfoque es que no es sencillo seleccionar los valores de wi , por lo cual en general se requieren varios intentos de ensayo y error para ajustarlos. Los parámetros ci , γ i se calculan de la misma forma que en el método anterior. Método de restricción ε Este método se conoce igualmente como método de las desigualdades (Zakian & Al-Naib, 1973). Consiste en resolver el siguiente problema escalar min f1(x) x∈X sujeto a:

(3.23)

fi (x) ≤ εi , i ∈ {2, 3, ..., M} En este caso es necesario experimentar sucesivamente con distintos valores del vector ε, relajando las restricciones que sean más dificiles de satisfacer, hasta conseguir una solución satisfactoria. Programación de metas Consiste en resolver el siguiente problema escalar M  min w i pi x∈X i=1

sujeto a:

fi (x) − pi = ti , i ∈ {1, 2, 3, ..., M} pi ≥ 0

(3.24)

41 donde los wi , ti , i ∈ {1, 2, 3, ..., M} son parámetros que deben ser definidos por el diseñador. Método de optimización lexicográfica En este método se establece un orden de preferencia (lexicográfico) entre los objetivos. Se resuelve en primer lugar el problema escalar min f1(x) (3.25) x∈X con α1 valor mínimo de f1. En una segunda etapa se intenta resolver el problema min f2 (x) x∈X sujeto a:

(3.26)

f1(x) = α1 y así sucesivamente. Por ejemplo, en una tercera etapa se intentaría resolver min f3 (x) x∈X sujeto a: f1(x) = α1

(3.27)

f2(x) = α2

3.2.2

Formulación de preferencias de manera interactiva

Un ejemplo de este tipo de métodos es el Interactive Weighted Tchebycheff Procedure (IWTP), el cual consiste en resolver una secuencia de problemas escalares del tipo min α α

sujeto a: α ≥ wT F(x)

(3.28)

x∈X Durante el proceso, el diseñador debe interactuar con el algoritmo, proporcionando el vector de pesos w más adecuado, con el fin de explorar las zonas del espacio objetivo que sean más interesantes a su criterio.

42

3.2.3

Formulación de preferencias “a posteriori”

Estos métodos generan un conjunto de potenciales soluciones al problema hasta que se cumpla alguna condición de parada. Al finalizar la ejecución, el diseñador debe seleccionar la solución definitiva entre las no-dominadas que se generaron durante el proceso de búsqueda. Un ejemplo sencillo consiste en generar vectores al azar, uniformemente distribuidos en el espacio de variables y luego filtrarlos para extraer el sub-conjunto de soluciones no-dominadas (ver algoritmo 3.1). Otro ejemplo que puede ser incluido dentro de esta categoría son los llamados algoritmos genéticos multi-objetivo, los cuales se describen en la siguiente sección.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

Entradas : P : Población generada por un algoritmo de optimización Salidas : A : Conjunto de soluciones no-dominadas extraídas de P i = 1; j = 1; Mientras i ≤ |P | hacer p1 ← P i ; f lag ← false; Mientras (j ≤ |P |) ∧ (f lag = f alse) hacer p2 ← P j ; Si p2 p1 entonces f lag ← true; finsi j = j + 1; Fin Si f lag = false entonces Incluir p1 en A; finsi i = i + 1; Fin Algoritmo 3.1 : Extracción del sub-conjunto de soluciones no-dominadas

3.3. Algoritmos genéticos multi-objetivo Los Algoritmos Genéticos (AG) ejecutan operaciones aleatorias sobre una población de individuos, semejantes a las que se realizan durante el proceso de la evolución biológica (mutaciones y recombinaciones). De manera similar a lo que ocurre en la naturaleza, los individuos más adaptados tienen mayor probabilidad para sobrevivir y reproducirse con respecto a los menos aptos. Jimenez y Sánchez (2002) identifican los siguientes elementos

43

Figura 3.1: Esquema general de un algoritmo genético característicos de un AG (ver figura 3.1): • Una representación del espacio de soluciones. • Un método para crear nuevos individuos. • Un criterio para comparar individuos entre sí, usualmente expresado mediante una función de adaptación. • Un método para seleccionar los individuos más adaptados. • Un conjunto de operadores de variación. • Un conjunto de parámetros que guían la evolución del algoritmo: tamaño de la población, número de iteraciones, probabilidades de aplicación de los operadores, etc. • Un criterio de parada. Tradicionalmente se considera que los AG realizan de manera adecuada la búsqueda global (exploración) pero no son eficientes para la explotación de las soluciones (Hart, 1994). Otra importante desventaja de estos algoritmos es que no existen reglas precisas para ajustar los parámetros que controlan la dinámica del algoritmo, ni para decidir el tipo de operador que conviene utilizar (Eiben & Smith, 2003). No obstante, numerosos autores han publicado reportes en los que afirman que este tipo de herramientas han permitido resolver problemas de optimización de alta complejidad, especialmente en aplicaciones de control automático (Fleming & Purshouse, 2002).

44

3.3.1

Operadores de variación

A continuación se describen algunos operadores de variación que transforman uno o más vectores reales en otro u otros vectores reales. Cada operador se aplica sobre uno o varios individuos con una probabilidad, denotada pc en el caso de los operadores de cruce y pm en el caso de los operadores de mutación. La terminología que se adopta en este capítulo fue propuesta en (Eiben & Smith, 2003). En particular, se denomina cromosoma a la estructura que contiene la información genética de un individuo, la cual puede venir expresada en diversos formatos alfa-numéricos. Cada cromosoma está compuesto por diversos genes, los cuales se asocian con las variables de decisión del problema de optimización. Operadores de cruce Los operadores de cruce generan nuevos individuos (hijos) a partir de la información de dos o más padres. De acuerdo con la naturaleza de la operación que llevan a cabo, es posible identificar dos tipos:

• Operadores discretos, caracterizados por uno, dos o más puntos de cruce, tal como se muestra en las figuras 3.2 y 3.3. Cada punto de cruce define una sección de cromosoma que será heredada por uno de los hijos. En el caso extremo, si cada gen se elige aleatoriamente a partir del cromosoma de uno de los padres, el operador de cruce se denomina uniforme.

• Operadores continuos, entre los cuales es posible mencionar el operador aritmético. Sean x1 y x2 los cromosomas correspondientes a dos individuos pertenecientes a la población. El operador de cruce aritmético genera dos nuevos individuos x1 y x2 mediante x1 = αx1 + (1 − α)x2 x2 = (1 − α)x1 + αx2

(3.29)

donde α es un escalar aleatorio uniformente distribudo en el intervalo [0, 1]. Note que si la región factible es convexa y los padres x1 , x2 son factibles, los hijos x1, x2 generados por este operador son igualmente factibles.

45

Figura 3.2: Operador discreto de cruce de un punto. El punto de cruce se encuentra en la mitad del cromosoma.

Figura 3.3: Operador discreto de cruce de dos puntos. Un operador continuo citado frecuentemente en la literatura es el Simulated Binary Cross-Over (SBX) propuesto en (Deb & Agrawal, 1995). Este operador produce dos hijos a partir de dos padres mediante las ecuaciones : x1i =

1 (1 + β i )x1i + 2 1 (1 − β i )x1i + x2i = 2 con

1 (1 − β i )x2i 2 1 (1 + β i )x2i 2

 1  (2ri ) ηc +1 , si ri ≤ 1 2   η 1+1 βi = c  1 2(1−ri )

(3.30)

(3.31)

donde ri es una variable aleatoria distribuida uniformemente en el intervalo [0, 1] y η c es un coeficiente que regula la distancia entre los cromosomas de los padres y los hijos. Un valor η c ≫ 1 produce hijos que con alta probabilidad se encuentran muy cerca de sus padres en el espacio de variables.

46 Operadores de mutación En el mundo natural, una mutación es un fenómeno poco común. En mamíferos la tasa de mutación es de 1 en 2.2 × 109 bases nucleotídicas (Kumar & Subramanian, 2002) mientras que en los virus es del orden de 1 en 100 (Malpica et al., 2002). En la mayoría de los casos las mutaciones son letales para los organismos (los mutantes no suelen alcanzar la edad reproductiva) pero a lo largo del proceso evolutivo, contribuyen a incrementar la diversidad genética de la especie. De cualquier forma, el operador de mutación se aplica únicamente a un porcentaje reducido de los genes totales de la población. Este porcentaje se puede seleccionar de dos maneras: de forma fija, especificando el mismo porcentaje de mutación en todas las generaciones del algoritmo genético y de forma variable, por ejemplo reduciendo el porcentaje de una generación a otra. De esta manera, se intenta realizar una búsqueda más amplia al principio del algoritmo y más restringida en las generaciones finales. Un operador de mutación puede ser representado de la forma #x1 , x2 , . . . , xn $ → #x1 , x2 , . . . , xn $

(3.32)

donde el vector x representa al individuo original y x representa al mutante. Un tipo común de mutación es la perturbación aletoria, definida mediante xi = xi + δ i , ∀i ∈ {1, 2, . . . , n}

(3.33)

donde cada δ es un vector aleatorio de distribución uniforme, gaussiana u otra.

3.3.2

Esquemas de selección

Por convención, se denota µs al número de padres que conforman la población actual, a partir de los cuales se generan λs hijos. En un algoritmo de tipo (µs , λs ) los hijos son directamente seleccionados para ser parte de la próxima generación (sin competencia en contra de sus padres). Si λs < µs es necesario completar a partir de los antiguos padres. Si λs > µs es posible eliminar los peores hijos hasta que queden exactamente µs . Por el contrario, en un algoritmo de tipo (µs + λs) se establece una competencia entre padres e hijos para seleccionar los µs sobrevivientes.

47 Se han propuesto numerosos métodos para la selección de un sub-conjunto de individuos a partir de una población inicial. La mayoría se basa en la definición de una función de adaptación (fitness) F it(i) ≥ 0 para cada individuo i, lo cual permite realizar comparaciones. Por ejemplo, el torneo binario consiste en seleccionar al azar (con probabilidad uniforme) dos individuos pertenecientes a la población. Sean F it(i) y F it(j) los valores de la función de adaptación correspondientes a los individuos i, j. El torneo binario selecciona el individuo i si y sólo si se cumple F it(i) ≥ F it(j) (se dice en este caso que i está más adaptado que j). En caso contrario se selecciona el individuo j. Otro método consiste en asignar a cada individuo i una probablidad de selección proporcional a su función de adaptación de la forma F it(i) pi = µs i=1 F it(i)

(3.34)

Sin embargo, el método descrito en el párrafo anterior presenta desventajas debido a que suele conducir a dos situaciones extremas: convergencia prematura o convergencia lenta. Por esta razón, Eiben y Smith (2003) recomiendan una alternativa basada en el ordenamiento de la población en valores decrecientes de la función de adaptación y la asignación de la probabilidad de selección mediante la fórmula pi =

2 − s 2(µs − r(i) + 1)(s − 1) + µs µs (µs + 1)

(3.35)

donde r(i) es el rango del individuo en la lista ordenada y s ∈ [1, 2] es un parámetro que controla la presión de selección. Cuando s se acerca a 2, mayor presión se ejerce para seleccionar los mejores individuos. Por el contrario, cuando s se acerca a 1, todos los individuos tienen idéntica probabilidad de selección.

3.3.3

Métodos de archivo

El problema de la gestión de archivo se presenta de manera natural cuando se intenta generar una aproximación finita, lo más cercana y bien distribuida posible con respecto al verdadero frente de Pareto. La idea básica para la gestión del archivo consiste en comparar la solución actual con respecto a las que ya existen en el archivo, para posteriormente decidir si es conveniente o no incluir la nueva, eliminando las soluciones dominadas por esta última.

48 Para comprobar si una solución pertenece o no al archivo es posible considerar un criterio de similitud, ya sea en el espacio objetivo o en el espacio de variables. El objetivo consiste en evitar que se archiven soluciones duplicadas o muy similares, de tal manera que se preserve la diversidad genética y se mejore la distribución de las soluciones. También es posible considerar cualquier otro criterio que se estime conveniente para favorecer el tipo de soluciones que se desea archivar, o utilizar conceptos especiales de dominancia, que incluyan las preferencias del diseñador. Cuando el archivo sobrepasa su capacidad máxima es necesario eliminar algunas soluciones, tomando en cuenta por ejemplo una medida de densidad, para favorecer a las soluciones pertenecientes a las regiones menos pobladas.

3.3.4

Tratamiento de las restricciones

Un problema de optimización multi-objetivo con restricciones puede plantearse como sujeto a  T min f1(x) f2(x) . . . fM (x) x∈X sujeto a:

g1 (x) ≤ 0 g2 (x) ≤ 0 .. .

(3.36)

gR (x) ≤ 0 Una estrategia frecuentemente considerada consiste en transformar el problema en uno sin restricciones mediante la incorporación de una función de penalidad por cada objetivo, de la forma min f1 (x) + p1(x) x∈X min f2 (x) + p2(x) x∈X .. .

(3.37)

min fM (x) + pM (x) x∈X donde usualmente se toma pi (x) =

R  j =1

wij [max (0, gj (x))] , i = 1, 2, . . . M

(3.38)

49 y los parámetros wij son pesos que deben ser ajustados para ponderar el efecto de cada restricción sobre cada función objetivo. Precisamente, una desventaja de este enfoque es que puede ser difícil seleccionar estos pesos de manera que el algoritmo explore la región factible de manera eficiente (Coello et al., 2007). Algunos investigadores han propuesto transformar las restricciones del problema en R nuevos objetivos, de manera que el problema se reformula como min f1(x) x∈X min f2(x) x∈X .. . min fM (x) x∈X

min g1(x) x∈X min g2(x) x∈X .. .

(3.39)

min gR (x) x∈X

Otra forma de atacar el problema de las restricciones consiste en modificar los esquemas de selección de individuos, para dar prioridad a los individuos factibles con respecto a los no factibles. También es posible utilizar operadores de reparación los cuales utilizan algún mecanismo para transformar un individuo no factible en factible. En general se recomienda no eliminar por completo a todos los individuos no factibles, con el objeto de permitir la captura de las regiones del conjunto de Pareto que pertenezcan a las fronteras de la región de factibilidad. A continuación se presentan tres algoritmos frecuentemente citados en la literatura y que, con numerosas variantes y mejoras, conforman el estado del arte en el área de la optimización genética multi-objetivo.

3.3.5

Multi-Objective Genetic Algorithm (MOGA)

Propuesto originalmente por Fonseca (1995), este algoritmo ha sido utilizado en numerosas aplicaciones de control (Chipperfield & Fleming, 1996; Chipperfield et al, 1998; Duarte et al, 2000). Por ejemplo, en (Griffin et al., 2000) se utilizó para resolver un problema de diseño caracterizado por 10 objetivos y 12 variables de decisión. En (Adra et al., 2005a) se propuso una versión híbrida para optimizar el sistema de control de una aeronave. También en (Silva et al., 2008) se utilizó para diseñar el controlador H∞ de una turbina de gas. Estos autores redujeron los tiempos de simulación mediante aproximaciones del proceso no-lineal.

50 El funcionamiento de MOGA se explica a continuación (ver algoritmo 3.2). A cada individuo se le asigna un rango rMOGA de acuerdo al número de otros individuos que lo dominan.

1 2 3 4 5 6 7 8 9 10

Entradas : P0 : Población inicial N : Número máximo de generaciones µs número de padres que conforman la población λs número de hijos a ser generados pc probabilidad de cruce pm probabilidad de mutación σ share : Tamaño deseado para cada sub-población Salidas : P : Aproximación del Conjunto de Pareto k = 1; Mientras k ≤ N hacer Pk = Pk−1 ; Calcular la función de adaptación para cada individuo de Pk de acuerdo con la ecuación 3.41; Seleccionar los padres; Aplicar los operadores de variación para generar los hijos; Calcular la función de adaptación para cada individuo de Pk incluyendo a los hijos generados en el paso anterior; Seleccionar los sobrevivientes y almacenarlos en Pk ; k = k + 1; Fin Algoritmo 3.2 : Multi-Objective Genetic Algorithm (MOGA) De esta forma, el rango rMOGA de un individuo i ∈ {1, 2, . . . µs } viene dado por: rMOGA (i) = 1 + psup (i)

(3.40)

donde psup(i) es el número de individuos que dominan a i (los no-dominados tendrán rango 1). En la figura 3.4 se muestra un ejemplo del rango correspondiente a una población de 7 individuos, en un problema de dos objetivos. Para seleccionar los mejores individuos, MOGA adopta un mecanismo llamado fitness sharing que consiste en calcular la función de adaptación de cada individuo mediante la fórmula µ − rMOGA (i) + 1 F it(i) = s (3.41) µs  µs Φ(dij ) j =1

 2   1 − dij σshare , si dij ≤ σ share Φ(dij ) =  0, si d > σ ij

share

(3.42)

51

Figura 3.4: Ejemplo de cálculo de rango para el algoritmo MOGA donde

  dij = Fi − Fj 2

(3.43)

F es el vector de funciones objetivo y σ share representa la distancia mínima deseada (en el espacio objetivo) entre dos soluciones pertenecientes a la aproximación final del Frente de Pareto. Este mecanismo penaliza a los individuos que habitan en zonas densamente pobladas, intentando conservar la diversidad genética y mejorar la distribución.

3.3.6

Non-dominated Sorting Genetic Algorithm (NSGA-II)

Este algoritmo (Deb et al., 2000) es uno de los más ampliamente citados en la literatura. En cuanto a las aplicaciones de control, en su tesis de doctorado Lagunas (2004) lo utiliza para sintonizar los parámetros de un controlador PID de dos grados de libertad. Por su parte, Dupuis et al (2004) presentan el diseño de un controlador PID cuyos objetivos se plantean en función de las diferencias entre los índices deseados y los índices simulados. También, Reyes et al (2009) presentan el diseño multi-objetivo de un controlador difuso caracterizado por 50 variables de decisión. El NSGA-II asigna un rango rNSGA−II (i) a cada individuo i ∈ {1, 2, . . . µs}, calculado como se explica a continuación. En cada iteración j (ver algoritmo 3.4) se determinan los individuos no-dominados y se les asigna rango j. Seguidamente se eliminan los individuos no-dominados de la sub-población actual y se comienza una nueva iteración (ver algoritmo 3.3). En la figura 3.5 se muestra un ejemplo de los rangos correspondientes a la misma población presentada en la figura 3.4.

52

1 2 3 4 5 6 7

Entradas : P0 : Población inicial Salidas : rNSGA−II : Vector de rangos correspondiente a P0 j = 1; Mientras |Pj | > 0 hacer Generar el conjunto Qj de vectores no-dominados pertenecientes a Pj ; Actualizar el vector rNSGA−II asignando rango j a los individuos pertenecientes a Qj ; Pj+1 = Pj /Qj ; j = j + 1; Fin Algoritmo 3.3 : Asignación de rango NSGA-II

Rango 2

Rango 3

Rango 1

Figura 3.5: Ejemplo de cálculo de rango para el algoritmo NSGA-II

La preservación de diversidad de NSGA-II se basa en calcular para cada individuo una cantidad llamada crowding distance, calculada como se explica a continuación (ver figura 3.6). En primer lugar se seleccionan los individuos que tengan el mismo rango. Luego se ordenan estos individuos en valores crecientes con respecto al objetivo i ∈ {1, 2, . . . M}. Sea ki el índice que corresponde a un individuo cualquiera en la lista ordenada con respecto al objetivo i, de forma que sus vecinos más cercanos tengan índices ki − 1 y ki + 1.

53

1 2 3 4 5 6 7 8 9 10

Entradas : P0 : Población inicial N : Número máximo de generaciones µs número de padres que conforman la población λs número de hijos a ser generados pc probabilidad de cruce pm probabilidad de mutación Salidas : P : Aproximación del Conjunto de Pareto k = 1; Mientras k ≤ N hacer Pk = Pk−1 ; Calcular el vector de rangos y crowding distance para cada individuo en Pk ; Seleccionar los padres mediante torneo binario modificado; Aplicar operadores de variación; Calcular el vector de rangos y crowding distance para cada individuo en Pk incluyendo los hijos generados en el paso anterior; Seleccionar los sobrevivientes mediante torneo binario modificado y almacenarlos en Pk ; k = k + 1; Fin Algoritmo 3.4 : Non-Dominated Sorting Genetic Algorithm (NSGA-II)

crowding distance = ∞

d1

d2

crowding distance = d1 + d2

crowding distance = ∞

Figura 3.6: Ejemplo de cálculo de crowding distance para el algoritmo NSGA-II

54 El valor del crowding distance se calcula mediante la fórmula: crowding distance =

M  i=1

Fiki +1 − Fiki −1

(3.44)

El mecanismo de selección utilizado por NSGA-II es un torneo binario modificado: entre dos individuos se prefiere al que tenga menor rango. Si los rangos son iguales, entonces se prefiere al que tenga el mayor valor de crowding distance.

3.3.7

Strength Pareto Evolutionary Algorithm (SPEA2)

Este algoritmo fue propuesto en (Zitzler et al., 2001) como una versión mejorada del algoritmo SPEA (Zitzler, 1999). Pocas aplicaciones se han publicado en el área de sistemas de control. Una de ellas fue presentada por Wang et al (2008) en la cual describen el diseño de un controlador robusto para un sistema de calentamiento a base de vapor. El SPEA2 utiliza un archivo externo A para almacenar los individuos no dominados (ver algoritmo 3.5) generados durante el proceso de búsqueda.

1 2 3 4 5 6 7 8 9 10

Entradas : P0 : Población inicial N : Número máximo de generaciones µs número de padres que conforman la población λs número de hijos a ser generados pc probabilidad de cruce pm probabilidad de mutación Salidas : A k = 1; Mientras k ≤ N hacer Calcular la función de adaptación F it(i) para los individuos pertenecientes a Pk y A; Seleccionar los padres mediante torneo binario modificado; Aplicar operadores de variación; Calcular la función de adaptación F it(i) para los individuos pertenecientes a Pk , A y los hijos generados en el paso anterior; Actualizar el archivo A; Seleccionar los sobrevivientes mediante torneo binario modificado y almacenarlos en Pk ; k = k + 1; Fin Algoritmo 3.5 : Strength Pareto Evolutionary Algorithm 2

55 A cada individuo i se le asigna un valor S(i) llamado fuerza (del inglés, strength) que se calcula como:   (3.45) S(i) =  j : j ∈ P ∪ A ∧ Fi Fj 

es decir es el número de individuos que son dominados por i en el conjunto unión de P y A. Posteriormente se calcula la suma de las fuerzas de los individuos que dominan a i obteniendo  S(j) (3.46) R(i) = j j ∈P ∪P,z zi En la figura 3.7 se presenta un ejemplo del cálculo de estos valores. S=1 R=2

S=0 R=13

S=2 R=0 S=1 R=6 S=3 R=0 S=3 R=0

S=1 R=8 S=2 R=0

Figura 3.7: Ejemplo de cálculo de S y R para el algoritmo SPEA2 Por otra parte se calcula un valor de densidad D(i) para cada individuo, definida como: D(i) =

1 2 + σi

(3.47)

donde σ i corresponde a la distancia entre el individuo i y su k-ésimo vecino más cercano  (ver figura 3.8). El valor de k se toma como el entero más cercano a |P | + |A|. Finalmente, la función de adaptación F it(i) de cada individuo se obtiene como F it(i) =

1 R(i) + D(i)

(3.48)

Para seleccionar los padres selecionan todos los individuos cuyo F it(i) sea mayor que 1 (nodominados). Si este número es menor que λs , se completa con individuos cuyo F it(i) sea el mayor posible. En caso contrario se utiliza un mecanismo para eliminar los individuos que habiten en las zonas de más alta densidad.

56

k-ésimo vecino más cercano distancia

3 1

σ

2

Población de 9 individuos → k = 3

Figura 3.8: Ejemplo de cálculo de la distancia σ para el algoritmo SPEA2 Para finalizar esta sección, en la tabla 3.1 se presenta un resumen de las ventajas y desventajas de los tres algoritmos descritos anteriormente.

Tabla 3.1: Resumen de ventajas y desventajas SPEA2 Algoritmo Ventajas MOGA Utiliza un mecanismo de asignación de rango sencillo NSGA-II Utiliza un mecanismo rápido para la asignación de rango. No necesita especificar parámetros adicionales SPEA2 Utiliza un archivo externo para almacenar las soluciones no dominadas. No necesita especificar parámetros adicionales

de los algoritmos MOGA, NSGA-II y Desventajas El usuario necesita especificar el parámetro σ share No integra la gestión de un archivo externo para las soluciones nodominadas Es más lento comparado con MOGA Y NSGA-II

57

3.4. Evaluación de desempeño El tema de la evaluación del desempeño de algoritmos multi-objetivo se encuentra todavía en desarrollo. En esta sección se presentan en primer lugar los indicadores que serán utilizados en los próximos capítulos. Luego se describe brevemente la prueba de Wilcoxon, la cual permite la comprobación de hipótesis relacionadas con las medianas de dos conjuntos de datos.

3.4.1

Índices de desempeño

Básicamente, son tres los aspectos que se miden mediante los índices que se han propuesto (Coello et al., 2007): • Cantidad de soluciones no-dominadas. • Distancia entre las soluciones no-dominadas y las soluciones del COGP. • Distribución de las soluciones no-dominadas: se desea obtener una distribución lo más uniforme y completa que sea posible (que cubra la mayor parte del COGP). Es interesante notar que la evaluación de las aproximaciones es también un problema de naturaleza multi-objetivo. Por ello, los expertos recomiendan usar varios indicadores para evaluar los distintos aspectos del desempeño de un algoritmo. En lo que sigue se denota P a una aproximación cualquiera de P ∗ y PF a la aproximación correspondiente de PF ∗ . Distancia generacional Este indicador, denotado GD, mide la distancia promedio de las soluciones pertenecientes a la aproximación PF con respecto a PF ∗ .  |PF|  1   GD = |di |2 (3.49) |PF | i=1   (3.50) di = min ∗ zi − zk  zk ∈PF De igual forma es posible evaluar la Distancia Generacional Inversa (IGD) mediante la fórmula  |PF ∗ |  1  IGD = |di |2 (3.51) |PF ∗| i=1

58 (3.52)

  di = min zi − zk  zk ∈PF

Mientras menores sean ambas medidas, mejor será la aproximación correspondiente. Sin embargo, estos indicadores no son necesariamente compatibles: si uno tiene un valor cercano a cero, esto no implica que el otro también lo tenga. En la figura 3.9 se muestran dos ejemplos para ilustrar esta situación. La gráfica de la izquierda corresponde a un caso en cual se tendría un valor bajo para GD y alto para IGD. La gráfica de la izquierda muestra un ejemplo de la situación inversa. f2

f2

Caso: GD ↑ IGD ↓

Caso: GD ↓ IGD ↑

f1

f1

Figura 3.9: Ejemplo de dos aproximaciones con distintos valores para los indicadores GD, IGD Distancia máxima Es la distancia máxima medida entre dos soluciones cualesquiera pertenecientes a la aproximación. En otras palabras, este indicador mide la extensión de la aproximación del frente de Pareto (ver figura 3.10).   DMAX = max zi − zk  zi ,zk ∈PF

(3.53)

Cobertura Este indicador binario fue propuesto en (Zitzler, 1999), para comparar dos aproximaciones PF 1 y PF 2, sin necesidad de conocer el COGP. Se define como sigue a continuación: C (PF 1, PF 2 ) =

|{z2 ∈ PF 2 tales que ∃z1 ∈ PF 1 : z1 z2}| |PF 2|

(3.54)

59

f2

DMAX

f1 Figura 3.10: Ejemplo de cálculo para el indicador DMAX Mide la proporción de elementos de PF 2 que son dominados por al menos un elemento de PF 1 . Por ejemplo, C(PF 1, PF 2) = 1 implica que todos los elementos de PF 2 son dominados por al menos un elemento de PF 1 . Espaciado Propuesto en (Schott, 1995), este indicador mide la desviación de las distancias entre los vecinos más cercanos en el espacio objetivo. Se define como sigue a continuación:

di =

 |PF|   1   d i − d 2 ESS = |PF | i=1

(3.55)

min

(3.56)

zk ∈PF ,k=i

 i  z − zk 

|PF| 1  di ,d= |PF| i=1

En la tabla 3.2 se muestra un resumen de los diferentes indicadores junto con el valor ideal para cada uno de ellos. El símbolo ∞ significa que mientras más grande es el valor, mejor es la aproximación con respecto al indicador correspondiente.

60 Tabla 3.2: Resumen de indicadores de desempeño Indicador Símbolo Valor Ideal Distancia generacional GD 0 Distancia generacional inversa IGD 0 Distancia máxima DMAX ∞ Cobertura C 1 Espaciado ESS 0

3.4.2

Pruebas estadísticas

Los algoritmos genéticos pueden ser considerados fundamentalmente como procesos estocásticos. El mismo algoritmo, partiendo de las mismas condiciones iniciales, produce resultados distintos cada vez que se ejecuta, lo cual justifica el empleo de medidas y pruebas estadísticas para comparar los resultados producidos por dos algoritmos cualesquiera. En particular, la prueba de Wilcoxon utiliza como estadístico de decisión W la suma de los rangos de cada observación, siendo el rango el número que le corresponde a cada observación cuando se ordenan todas en orden ascendente. Suponga que se realizan siete ejecuciones de dos algoritmos, A1 y A2, comenzando a partir de la misma población inicial. Posteriormente se calculan los valores correspondientes a un indicador de desempeño I obteniendo los resultados de la tabla 3.3. Tabla 3.3: Ejemplo de valores obtenidos para el ciones de los algoritmos A1 y A2 I A1 1 0.44 2 0.55 3 0.61 4 0.53 5 0.56 6 0.45 7 0.54 Promedio 0.52

indicador I correspondientes a 7 ejecuA2 0.50 0.48 0.58 0.52 0.59 0.65 0.70 0.57

La prueba consiste en formular una hipótesis nula de la forma H0 : I A1 = I A2

(3.57)

61 y una hipótesis alternativa H1 : I A1 = I A2

(3.58)

con un nivel de significancia α, donde I A1 y I A2 son las medianas de cada proceso estocástico. En la figura 3.11 se muestra el cálculo de los estadísticos WA1, WA2 para el ejemplo de la tabla 3.3.

Figura 3.11: Cálculo de los rangos para la prueba Wilcoxon

Si el p-valor correspondiente a la prueba es menor que α entonces se rechaza la hipótesis R y es igual H0. Este valor puede obtenerse mediante la instrucción ranksum de MATLAB a 0.3829 para este caso, por lo cual no existe evidencia estadística para rechazar la hipótesis H0 con un nivel de significancia α = 0.05.

3.5. Resumen del capítulo En este capítulo se definieron algunos conceptos básicos relacionados con el área de optimización multi-objetivo y algoritmos genéticos, tales como conjunto y frente de Pareto, dominancia, direcciones de descenso, conos de diversidad, operadores de variación y esquemas de selección. Se realizó un breve repaso de las estrategias que pueden ser utilizadas para resolver un problema de esta naturaleza, las cuales pueden ser agrupadas en 3 categorías: • Formulación de preferencias a priori • Formulación de preferencias de manera interactiva • Formulación de preferencias a posteriori

62 Con respecto a esta última categoría, se presentó un breve análisis de los algoritmos MOGA, NSGA-II y SPEA2, los cuales son los más citados en la literatura y constituyen el estado del arte en esta área. Un punto importante para la evaluación del desempeño de estos algoritmos consiste en la definición de índices para medir la calidad de las aproximaciones generadas. En este capítulo se definieron índices tales como la distancia generacional, distancia generacional inversa, cobertura, espaciamiento, etc. Puesto que el proceso de búsqueda es fundamentalmente estocástico, se hace imprescindible el uso de pruebas estadísticas para la correcta evaluación de los resultados. La prueba de Wilcoxon presentada permite la comprobación de hipótesis relacionadas con las medianas de dos conjuntos de datos, sin necesidad de conocer las distribuciones probabilísticas asociadas. Esta prueba será utilizada en los próximos capítulos para evaluar el desempeño de los algoritmos propuestos.

CAPÍTULO IV DISEÑO DE CONTROLADORES PID MULTI-OBJETIVO

En este capítulo se proponen nuevos operadores de variación, especialmente concebidos para el diseño de controladores PID mediante algoritmos genéticos. La idea principal consiste en garantizar la generación y el mantenimiento de una alta proporción de individuos factibles durante todo el proceso de búsqueda. En primer lugar se utiliza el método propuesto en (Hohenbichler & Abel, 2008) para generar los vértices de la región factible correspondiente al sistema a lazo cerrado. Posteriormente se propone un nuevo método para aproximar la geometría de la región factible, basado en muestreo aleatorio uniforme (Calafiore & Dabbene, 2006). Parte de los resultados presentados en este capítulo fueron previamente incluidos como parte de las publicaciones (Sánchez et al., 2008a) y (Sánchez et al., 2009).

4.1. Antecedentes Es posible encontrar en la literatura numerosas referencias de trabajos previos relacionados con el problema que se intenta resolver en este capítulo, basados tanto en algoritmos deterministas como aleatorios. Con respecto a los primeros, Åström et al (1998) utilizaron el algoritmo de NewtonRaphson para minimizar la energía de la señal de error en un lazo de control PID sujeto a la acción de una perturbación de salida. Estos autores señalan que el resultado de la optimización depende fundamentalmente del punto de inicio seleccionado. Hwang & Hsiao (2002) resuelven el problema no-convexo de la siguiente forma: parten de un valor arbitrario para la constante derivativa y a partir de allí resuelven el problema de dos variables (proporcional/integral) incluyendo restricciones para asegurar estabilidad.

64 Granado et al (2007) presentan el diseño de un controlador PID multivariable robusto, el cual garantiza la estabilidad a lazo cerrado de sistemas lineales sujetos a incertidumbre poliédrica. El algoritmo de diseño se basa en resolver un conjunto de BMIs (Bilinear Matrix Inequalities) de manera iterativa. Tal como se mencionó en la introducción, la desventaja de este tipo de enfoques es que las soluciones encontradas pueden ser en ocasiones demasiado conservadoras y por otra parte únicamente consiguen una solución en cada ejecución (ver capítulo 6). Con respecto al empleo de algoritmos aleatorios, Herreros (2000) propuso dos nuevos algoritmos en su tesis de doctorado, a saber MRCD y MRCD-min-max. Este autor diseñó controladores PIDs para los sistemas de ensayo de bajo órden propuestos en (Astrom et al., 1998), consiguiendo mejores soluciones que las obtenidas mediante LMIs. En (Dupuis et al., 2004) se reporta el diseño de un PID para controlar la velocidad de giro de un motor de corriente continua, utilizando el algoritmo NSGA-II. Los objetivos fueron planteados en función de la diferencia entre la dinámica obtenida durante la simulación y la dinámica deseada. Se tomó en cuenta el criterio de estabilidad de Routh-Hurwitz a efectos de generar individuos estables. Por último, en su tesis doctoral Lagunas (2004) utilizó el NSGA-II para diseñar controladores PIDs en el caso de los sistemas de ensayo SISO propuestos en (Astrom et al., 1998). Se consideraron los siguientes objetivos: seguimiento de la señal de referencia, rechazo de perturbación en la carga, atenuación de ruido en la medición, robustez ante variaciones de los parámetros de la planta y acotamiento de la señal de control.

4.2. Método de diseño basado en el cálculo exacto de la región de estabilidad Los algoritmos genéticos operan en el interior de un espacio de búsqueda finito, generalmente definido por el diseñador, el cual debe incluir la región factible para que el proceso de búsqueda sea capaz de conseguir soluciones válidas. Típicamente, la población inicial es generada mediante muestreo aleatorio uniforme en el espacio de búsqueda. Esto significa que cada cromosoma xi , i = 1, 2, . . . µs, es generado mediante la fórmula

  xij = Lj + rj Lj − Lj , j = 1, 2, . . . n

(4.1)

65

donde L, L ∈ Rn son vectores que definen el intervalo de búsqueda para cada variable de decisión y cada ri es una variable aleatoria distribuida uniformemente en [0, 1]. Sin embargo, este método puede resultar bastante ineficiente para producir una alta proporción de individuos factibles, en aplicaciones de diseño de controladores. Para ilustrar este hecho, suponga que se desea generar de manera aleatoria un conjunto de controladores PID, caracterizados por su función de transferencia KPID (s) =

KD s2 + KP s + KI s

(4.2)

para estabilizar el modelo (Ackermann & Kaesbauer, 2003): G0 (s) =

−0.5s4 − 7s3 − 2s + 1 A(s) = B(s) (s + 1)(s + 2)(s + 3)(s + 4)(s2 + s + 1)

(4.3)

el cual posee una complejidad superior a los que típicamente se han analizado en trabajos previos (órden 6 y fase no-mínima). Sean L, L ∈ R3 los vectores que definen el intervalo de variación para cada parámetro PID, de forma que KP ∈ [L1 , L1]

(4.4)

KI ∈ [L2 , L2] KD ∈ [L3 , L3] Suponga que se toman N muestras en el interior del volumen definido por L, L y sea Ns (L, L) el número de individuos factibles correspondiente, es decir: Ns (L, L) = con

N  i=1

St(KP IDi )

  0 si max {Re [λ(KPID )]} ≥ −ε St (KPID ) =  1 en caso contrario

(4.5)

(4.6)

donde ε > 0 es un parámetro que define la tolerancia considerada para la región de estabilidad y λ(KP ID ) denota el vector de raices de la ecuación B(s) + KP ID (s)A(s) = 0

(4.7)

La tabla 4.1 muestra el promedio de la proporción de individuos factibles NNs para distintos valores de L, L, calculada para N = 1000, ε = 0.01 y repitiendo el experimento un total

66 de 100 veces. Note que la mejor proporción (≈ 40%) se logra en el caso L = [−1, −1, −1] y L = [1, 1, 1] (el valor mostrado entre paréntesis corresponde a la desviación standard del cociente NNs ). Tabla 4.1: Proporción de individuos factibles para distintos valores de L, L Ns L L N [−0.1, −0.1, −0.1] [0.1, 0.1, 0.1] 0(0) [−1, −1, −1] [1, 1, 1] 0.3871(0.0155) [−10, −10, −10] [10, 10, 10] 0.1440(0.0130) [−100, −100, −100] [100, 100, 100] 0(0) [−1000, −1000, −1000] [1000, 1000, 1000] 0(0)

Es posible concluir a partir de estos resultados que definir el espacio de búsqueda de manera arbitraria puede ser muy ineficiente. Precisamente con el objeto de evitar esta situación, se propone a continuación un método alternativo que permite resolver el problema de diseño en dos etapas consecutivas.

Figura 4.1: Espacio de controladores PID estabilizantes correspondiente al modelo G0(s) determinado mediante el método de Hohenbichler y Abel (2008) Durante la primera se determinan los límites de la región de controladores PID estabilizantes (ver figura 4.1) utilizando el método de Hohenbichler y Abel (2008) presentado en el capítulo 2.

67 Durante la segunda etapa se genera una distribución uniforme de individuos en el interior de los polígonos determinados por los límites calculados durante la etapa anterior, utilizando el método que se explica a continuación. La información de base son los vértices de los polígonos de estabilidad. Suponga que para un polígono de V vértices se ordena esta información de la forma        

v1 v2 .. . vV

       

(4.8)

donde v1 está conectado con v2 y vV , v2 está conectado con v1 ,v3 y así sucesivamente. En primer lugar se seleccionan dos aristas del polígono al azar. A su vez, en cada arista se selecciona un punto al azar, sean xa y xb . Posteriormente se selecciona un punto xc perteneciente al segmento que une xa con xb (ver figura 4.2) mediante la fórmula xc = αxa + (1 − α)xb

(4.9)

donde α es un escalar aleatorio perteneciente al intervalo [0, 1]. En el caso de la región estabilizante correspondiente a G0 (s) se determinaron 25 polígonos de estabilidad y se generaron 2 individuos al interior de cada uno de ellos, para obtener una población inicial de 50 individuos.

Figura 4.2: Selección de un punto al azar en el interior de un polígono convexo En cuanto a los operadores de variación, note que aunque para un mismo valor de KP la región factible es convexa, es posible que al cruzar dos individuos con distinto valor de

68 KP se produzcan hijos no factibles. En este caso se propone utilizar el siguiente mecanismo de reparación: • Suponga que se conocen los dos padres x1y x2, a partir de los cuales se generan los descendientes x1 y x2 mediante el operador de cruce aritmético de la forma: x1 = αx1 + (1 − α)x2 x2 = (1 − α)x1 + αx2

(4.10)

donde α es un escalar aleatorio perteneciente al intervalo [0,1]. • Si cualquiera de los descendientes no es factible, entonces se disminuye α a la mitad y se genera un nuevo par de hijos. • La mutación se realiza de forma similar, fijando α = 1 y a partir de un individuo factible x se genera una mutación mediante x = x + αδ

(4.11)

donde δ es un vector aleatorio de distribución gaussiana. Si el individuo resultante no es factible entonces se disminuye α a la mitad y se genera una nueva mutación.

Considere el lazo de control mostrado en la figura 4.3. Los bloques Wn (s) y Wd(s) corresponden a funciones arbitrarias de ponderación en frecuencia. Las variables que intervienen en el lazo son r(s), d(s) y n(s) como señales de referencia, perturbación y ruido respectivamente.

Figura 4.3: Lazo con controlador PID

69 Tradicionalmente, para el problema de regulación se considera r(s) = 0 y por lo tanto a lazo cerrado se tiene: KP ID (s)G(s) Wn (s)n(s) 1 + KPID (s)G(s) 1 + Wd (s)d(s) 1 + KP ID (s)G(s) = −T (s)Wn (s)n(s) + S(s)Wd (s)d(s)

y(s) = −

con T (s) = S(s) =

KPID (s)G(s) 1 + KP ID (s)G(s) 1 1 + KPID (s)G(s)

(4.12) (4.13) (4.14)

(4.15) (4.16)

La función de transferencia S(s) se conoce como función de sensibilidad, mientras que T (s) se conoce como función de sensibilidad complementaria. Note que siempre se cumple S(s) + T (s) = 1 (4.17) A partir de la ecuación (4.12) se concluye que para compensar los efectos de la perturbación y del ruido, es necesario diseñar el controlador para que las normas de S(s) y T (s) tengan el menor valor posible. Sin embargo, la ecuación (4.17) implica que no se puede disminuir ambas normas simultaneamente. El modelo de planta a lazo abierto es G0 (ver ecuación 4.3) y los filtros Wn y Wd se definen como: 1 Wn(s) = 1, Wd(s) = (4.18) s+1 El problema que se desea resolver se plantea como min f1 (KPID )

(4.19)

min f2 (KPID )

(4.20)

KP ID ∈K KP ID ∈K

con f1(KP ID ) = Wd (s)S(s)2 f2(KP ID ) = Wn (s)T (s)2 donde K es el conjunto de controladores estabilizantes.

(4.21)

70 El problema de diseño que involucra exclusivamente normas H2 en los distintos canales de la función de transferencia a lazo cerrado, puede ser resuelto de manera eficiente mediante métodos de programación convexa, en caso de que no se imponga ninguna estructura especial del controlador. Por ejemplo, Khargonekar y Rotea (1991) resolvieron este problema asumiendo una estructura de orden completo para los controladores y considerando un único objetivo en términos de la suma convexa de las normas cuadráticas, de la forma min α f1(K) + (1 − α)f2 (K)

K ∈R1×n

(4.22)

donde para cada valor α ∈ [0, 1] se obtiene una solución perteneciente al frente de Pareto del problema multi-objetivo. En el presente caso, la dificultad radica en que se desea obtener controladores con una estructura seleccionada a priori de tipo PID, lo cual dificulta el planteamiento en términos de programación convexa. Como base del proceso evolutivo se adopta el algoritmo SPEA2 (ver capítulo 3) considerando que: • Permite seleccionar la cantidad máxima y las características de las soluciones de la aproximación, gracias a un manejo adecuado del archivo externo. • Optimiza la distribución de las soluciones en el frente, sin necesidad de ajustar parámetros dependientes de su geometría. • Su desventaja con respecto al NSGA-II y al MOGA es que su tiempo de ejecución, para un mismo problema, es relativamente más lento. Sin embargo, para el tipo de problemas que se consideran en este trabajo, la velocidad de cómputo no se considera una variable determinante. En concreto, se presentan a continuación los resultados obtenidos a partir de dos versiones de SPEA2: • A1: El cual opera en un espacio de búsqueda definido de manera arbitraria como −10 ≤ KP , KI , KD ≤ 10

(4.23)

y utiliza una función de penalidad en cada objetivo para manejar la restricción de estabilidad.

71 • A2: El cual opera en el espacio de controladores estabilizantes determinado mediante el método de Hohenbichler y Abel, utilizando operadores especiales de cruce y mutación (ver ecuaciones 4.10 y 4.11). Se realizaron 30 ejecuciones de cada algoritmo y se almacenaron los resultados obtenidos. En la figura 4.4 se muestra un ejemplo de las aproximaciones del frente de Pareto obtenidas por cada algoritmo. En dicha figura se puede apreciar que A2 produce individuos (marcados con
) que claramente dominan a las obtenidos por A1 (marcados con ×). 0.65 Aproximación obtenida por A1 Aproximación obtenida por A2

0.6 0.55 0.5

Controlador K1

f2

0.45 0.4 Controlador K2 0.35 0.3 0.25 0.2 0.66

0.67

0.68

0.69 f1

0.7

0.71

0.72

Figura 4.4: Ejemplo de las aproximaciones obtenidas por A1 y A2 La configuración utilizada para las pruebas de A1 y A2 se muestra en la tabla 4.2. Note que se permitió un máximo de 2500 evaluaciones, distribuidas durante 50 generaciones, para cada algoritmo. Esto motivado a que previamente se verificó de manera experimental que estos números son suficientes para alcanzar una etapa del proceso en la cual las aproximaciones de los frentes evolucionan lentamente con respecto al inicio. En las figuras 4.5 y 4.6 es posible apreciar un ejemplo de la evolución de las aproximaciones producidas por A1 y A2 después de haber transcurrido 1,25,40,48 y 50 generaciones. Note que en ambos casos, las tres últimas generaciones corresponden a aproximaciones muy parecidas entre sí. Por otra parte, si se permite una magnitud mayor en el número de evaluaciones (p.ej. 10.000 o más) entonces en la práctica no es posible detectar diferencias significativas entre los indicadores de desempeño producidos por cada algoritmo y no sería posible evaluar el efecto de los nuevos operadores propuestos.

72 1.7 Ngen = 1 Ngen = 25 Ngen = 40 Ngen = 48 Ngen = 50

1.6 1.5 1.4

f2

1.3 1.2 1.1 1 0.9 0.8 0.68

0.7

0.72 f

0.74

0.76

0.78

1

Figura 4.5: Evolución de las aproximaciones de los frentes. Algoritmo A1. Ngen es el número de generación. Tabla 4.2: Configuración utilizada para los algoritmos A1 y A2. Problema de diseño PID Parámetro Valor Representación Números reales Operador de cruce Aritmético (A1) y Aritmético con retorno (A2) Probabilidad de cruce 0.9 Operador de mutación Perturbación gaussiana con retorno a la región factible Probabilidad de mutación 0.1 Condición de parada 50 generaciones Número máximo de evaluaciones 2500 Tamaño máximo del archivo externo 50 Esquema de selección (50+50)

Ngen = 1 Ngen = 25 Ngen = 40 Ngen = 48 Ngen = 50

0.6 0.55 0.5

f

2

0.45 0.4 0.35 0.3 0.25 0.65

0.7

0.75 f

1

Figura 4.6: Evolución de las aproximaciones de los frentes. Algoritmo A2. Ngen es el número de generación.

73 Los resultados obtenidos para los indicadores C, N, DMAX y ESS (definidos en el capítulo 3) se muestran en las tablas 4.3 y 4.4. Tabla 4.3: Resultados obtenidos para el indicador de cobertura. Algoritmos A1 y A2 C(Ai , Aj ) A1 A2 A1 − 0.2043(0.3103) ∗ A2 0.5946(0.4552) − Tabla 4.4: Resultados de los indicadores ESS y DMAX. Algoritmos A1 y A2 N DMAX ESS A1 15.4667(15.1128) 0.1767(0.1559) 0.0118(0.0224) A2 36.5333(14.0804)∗ 0.3692(0.1629)∗ 0.0080(0.0082) Tabla 4.5: P-valores correspondientes a las pruebas de Wilcoxon para determinar las diferencias que son estadísticamente significativas. Algoritmos A1, A2 Hipótesis Nula p − valor Conclusión C(A1, A2) = C(A2, A1) 7.7540e − 004 Rechazar N (A1) = N (A2) 1.9066e − 006 Rechazar DMAX(A1) = DMAX(A2) 2.4321e − 005 Rechazar ESS(A1) = ESS(A2) 0.4118 No Rechazar El indicador N representa simplemente la cantidad de individuos almacenados en el archivo de soluciones no-dominadas al finalizar la ejecución del algoritmo. El asterisco corresponde al indicador que es significativamente mejor de acuerdo a la prueba de Wilcoxon, para un nivel de significancia α = 0.05. En la tabla 4.5 se presentan los p-valores correspondientes a la prueba de Wilcoxon, los cuales nos permiten identificar las diferencias entre los promedios que son significativas desde el punto de vista estadístico. Note que A2 obtiene mejores valores que A1 con respecto a los indicadores C, N y DMAX. Estos resultados podrían haberse anticipado, puesto que A2 utiliza información especial para mantener a la población dentro de la región factible durante todo el proceso de búsqueda. Con respecto al espaciamiento de las soluciones, el indicador ESS no permite identificar mayores diferencias, lo cual puede constatarse de manera gráfica observando la figura 4.4. Para ilustrar lo que representa la diferencia entre las aproximaciones obtenidas por A1 y A2 en términos prácticos, en la figura 4.8 se presenta una simulación del efecto de dos controladores K1 y K2 (ver figura 4.4) ante la acción de las señales de perturbación y ruido mostradas en la figura 4.7. En particular, es posible observar una reducción significativa de los picos de la señal de salida.

74

1

Perturbación p(t) Ruido n(t)

0.8

y(t)

0.6 0.4 0.2 0 −0.2 0

10

20

30

40

50 t

60

70

80

90

Figura 4.7: Perturbación y ruido que actúan sobre el lazo

1.4 1.2

Respuesta del lazo con controlador K1 Respuesta del lazo con controlador K2

1 0.8

y(t)

0.6 0.4 0.2 0 −0.2 −0.4

0

20

40

60

80

100

t

Figura 4.8: Salida del lazo de control, afectada por perturbación y ruido

75

4.3. Método de diseño basado en la aproximación la región de estabilidad Tal como hemos visto en las sección anterior, la restricción de estabilidad es fundamental para el diseño de controladores automáticos. Cuando se intenta resolver el problema mediante algoritmos genéticos, es posible utilizar una función auxiliar para penalizar los individuos inestables (no factibles) y tratar de guiar a la población hacia la región estabilizante. Sin embargo, esta estrategia puede conducir a fenómenos no deseados tales como convergencia prematura en caso de que sea difícil encontrar soluciones factibles. Ante esta situación, los investigadores han planteado dos tipos de enfoque: • Utilizar un algoritmo que opere en un espacio de búsqueda con límites dinámicos, que puedan adaptarse para no quedar atrapado dentro de una región no factible. • Utilizar algun método para obtener información acerca de la geometría de las regiones factibles. Varios investigadores han propuesto espacios de búsqueda con límites dinámicos. De hecho el algoritmo Multi-Objective Robust Control Design (MRCD) propuesto en (Herreros, 2000) utiliza límites dinámicos para el espacio de búsqueda y una división de la población en islotes distribuidos. En (Arakawa et al., 1998) se propone el algoritmo Adaptive Range Genetic Algorithm (ARange GA) en el cual se adopta un espacio de búsqueda adaptivo que no necesita especificar el rango para la población inicial. En (Khor et al., 2002) se utiliza un enfoque de aprendizaje inductivo/deductivo para dirigir a la población hacia zonas que se encuentran fuera de los límites iniciales de búsqueda. Los autores afirman que su método es capaz de adaptar la precisión de la representación genética con el fin de mejorar las soluciones. Otro algoritmo que puede citarse dentro de esta categoría es el Adaptive Range MultiObjective Genetic Algorithm (ARMOGA) propuesto en (Sasaki et al., 2002). Este algoritmo calcula estadísticas en función de las muestras que se han tomado para decidir los límites del espacio de búsqueda que se utilizarán en la siguiente iteración. Los autores

76 afirman que su método permite reducir el número de evaluaciones necesarias para obtener la aproximación de Pareto. En otro orden de ideas, los métodos aleatorios han ganado aceptación para el estudio de sistemas inciertos y una cantidad considerable de resultados han sido publicados en la literatura (Calafiore & Dabbene, 2006). La idea básica consiste en considerar los parámetros inciertos como variables aleatorias, lo que conduce a evaluar su desempeño en términos de probabilidades. De manera recíproca, la síntesis probabilística consiste en determinar los valores de los parámetros para que la probabilidad de obtener un buen desempeño sea alta. A continuación se propone un método de diseño que intenta determinar una aproximación de la geometría de la región factible mediante muestreo aleatorio puro. Es necesario aclarar que el objetivo no es el cálculo de la región estabilizante como tal, ya que en este caso el algoritmo de Hohenbichler y Abel permite determinarla de manera eficiente como se mostró anteriormente. El objetivo que se persigue es mostrar que en algunos casos es posible obtener aproximaciones que permiten mejorar la calidad de los diseños finales obtenidos por el algoritmo genético. Durante la primera etapa del método propuesto se resuelve un sub-problema que consiste en estimar los vértices del espacio de búsqueda que maximizan simultaneamente el volumen y la proporción de individuos factibles, muestreados aleatoria y uniformemente. Durante la segunda etapa se selecciona un conjunto de vértices en particular con el fin de continuar la búsqueda evolutiva considerando los objetivos de diseño del problema original y enfocándose en la región seleccionada. Sea n el número de variables de decisión del problema (4.19) y L, L ∈ Rn vectores reales que definen el intervalo de variación para cada una de las variables. En el caso del diseño PID hemos visto que se tiene n = 3 y además KP ∈ [L1 , L1]

(4.24)

KI ∈ [L2 , L2] KD ∈ [L3 , L3] Consideremos el siguiente sub-problema:





g1(L, L)  min  L,L∈Rn g2(L, L)

(4.25)

77 con

Ns (L, L) N n  Vmax − (Li − Li )

g1(L, L) = 1 −

(4.26)

i=1 (4.27) Vmax y donde Ns es el número de individuos factibles obtenidos mediante muestreo aleatorio en el interior del volumen definido mediante L, L, ver ecuación (4.5). El número N de elementos del conjunto de muestreo puede ser estimada mediante la fórmula de Chernoff:

g2 (L, L) =

N= y entonces se cumple:

donde

2 1 ln 2 2ǫ δ

    Ns   Prob p − ≤ǫ ≥1−δ N p = Prob {St (KPID ) = 1}

(4.28)

(4.29)

(4.30)

con ǫ, δ cantidades positivas que deben ser ajustadas en función de la precisión deseada. Por otra parte, en la ecuación (4.27) la cantidad Vmax es el volumen máximo del espacio de búsqueda, calculado como n

(Li − Li ) Vmax = max L,L∈Rn i=1

(4.31)

Para la pruebas se seleccionaron los valores ǫ = δ = 0.1 lo cual conlleva a un muestreo de 100 individuos por cada sub-región. Se consideró el mismo modelo a lazo abierto estudiado en la sección anterior (ver figura 4.3) con Vmax = 20 × 20 × 20. En la tabla 4.6 se muestra la configuración utilizada para la ejecución del SPEA2 durante esta primera etapa. Note que para este problema se permitió un total de 1000 evaluaciones de las funciones g1 y g2, lo cual es suficiente para seleccionar una aproximación de la región factible. La figura 4.9 permite apreciar un ejemplo del frente obtenido al finalizar esta etapa. En dicha figura también es posible apreciar los intervalos PID correspondientes a las dos soluciones extremas del frente de Pareto. Note el conflicto agudo entre los objetivos g1 y g2, en el cual practicamente no existe solución de compromiso.

78 Tabla 4.6: Configuración utilizada por SPEA2 durante la primera etapa del método de aproximación de la región factible Parámetro Valor Representación Números reales Operador de cruce Aritmético con reparación Probabilidad de cruce 0.9 Operador de mutación Perturbación gaussiana con reparación Probabilidad de mutación 0.1 Condición de parada 20 generaciones Número máximo de evaluaciones 1000 Tamaño máximo del archivo externo 50 Esquema de selección (50+50) 1.2 1

g2

0.8 KP∈ [−4.9,2.4] KI∈ [0.19,6.24] K ∈ [−5.9,3.2]

0.6

I

0.4 KP∈ [−20,20] KI∈ [−20,20] KI∈ [−20,20]

0.2 0 0

0.1

0.2

0.3

0.4

0.5 g1

0.6

0.7

0.8

0.9

Figura 4.9: Ejemplo de aproximación obtenida al final de la primera etapa Con respecto a los operadores de variación, a lo largo de esta primera etapa es necesario utilizar un mecanismo de reparación para corregir los individuos que no correspondan a un conjunto válido de vértices. Esta reparación consiste simplemente en ordenar la información contenida en los cromosomas para asegurar que aparezcan en forma ordenada. Una vez culminada la primera etapa, se selecciona un conjunto de vértices en particular con el fin de continuar el proceso de búsqueda considerando el problema de diseño original y enfocándose en la región delimitada por los vértices seleccionados. A continuación se presentan los resultados obtenidos mediante los siguientes algoritmos: • A1: que opera en un espacio de búsqueda definido por KP , KI , KD ∈ [−10, 10]

(4.32)

79 • A2: que opera en el espacio de controladores estabilizantes determinado mediante el método de Hohenbichler y Abel, junto con operadores especiales de variación. • A3: correspondiente al método de diseño en dos etapas presentado en esta sección, utilizando los límites (4.33)

KP ∈ [−4.95, 2.43] KI ∈ [0.19, 6.24] KD ∈ [−5.97, 3.25]

que corresponden a g1 = 0 y g2 = 0.99 (ver figura 4.9). Estos valores indican que se le da más importancia al objetivo g1 , relacionado con la proporción de individuos factibles muestreados. La configuración utilizada para A1, A2 y A3 es idéntica a la mostrada anteriormente en la tabla 4.2. En la figura 4.10 se puede observar un ejemplo de las aproximaciones que se obtiene mediante cada uno de los algoritmos. En las tablas 4.7 y 4.8 se presentan los valores promedios y las desviaciones correspondientes a los indicadores C, N, DMAX y ESS. 0.65 0.6 0.55 0.5 Aproximación obtenida por A1 Aproximación obtenida por A2 Aproximación obtenida por A3

f2

0.45 0.4 0.35 0.3 0.25 0.2 0.66

0.67

0.68

0.69 f1

0.7

0.71

0.72

Figura 4.10: Ejemplo de las aproximaciones que se obtienen mediante A1, A2 y A3 Como era de esperarse, el algoritmo A2 produce los mejores indicadores en cuanto al indicador de cobertura, para el cual A1 y A3 son estadísticamente equivalentes. Sin embargo, A3 produce mejores indicadores que A1, lo cual ilustra la eficacia del método de aproximación de la región factible para el problema propuesto. En la tabla 4.9 se presentan los p-valores correspondientes a las pruebas de Wilcoxon para determinar las diferencias que son estadísticamente significativas.

80

Tabla 4.7: Resultados del indicador de cobertura. Algoritmos A1, A2 y A3 C(Ai , Aj ) A1 A2 A3 A1 − 0.2043(0.3103) 0.0580(0.1790) ∗ A2 0.5946(0.4552) − 0.6074(0.2898)∗ A3 0.6483(0.4228) 0.0723(0.1680) −

Tabla 4.8: Resultados de los indicadores N, ESS y DMAX. Algoritmos A1, A2 y A3 N DMAX ESS A1 15.4667(15.1128) 0.1767(0.1559) 0.0118(0.0224) A2 36.5333(14.0804)∗ 0.3692(0.1629)∗ 0.0080(0.0082) A3 26.5333(14.7712) 0.2423(0.1050) 0.0074(0.0076)

Tabla 4.9: P-valores correspondientes a las pruebas de Wilcoxon para determinar las diferencias que son estadísticamente significativas. Algoritmos A1, A2 y A3 Hipótesis Nula p − valor Conclusión C(A1, A2) = C(A2, A1) 7.7540e − 004 Rechazar C(A1, A3) = C(A3, A1) 6.0177e − 008 Rechazar C(A2, A3) = C(A3, A2) 2.4955e − 009 Rechazar N (A1) = N (A2) 1.9066e − 006 Rechazar N (A1) = N (A3) 9.5896e − 004 Rechazar N (A2) = N (A3) 0.0096 Rechazar DMAX(A1) = DMAX(A2) 2.4321e − 005 Rechazar DMAX(A1) = DMAX(A3) 0.0199 Rechazar DMAX(A2) = DMAX(A3) 0.0025 Rechazar ESS(A1) = ESS(A2) 0.4118 No Rechazar ESS(A1) = ESS(A3) 0.4917 No Rechazar ESS(A3) = ESS(A2) 0.9352 No Rechazar

81

4.4. Resumen del capítulo En este capítulo se presentaron dos nuevos métodos de diseño PID multi-objetivo que utilizan resultados de la teoría de control y algoritmos genéticos para resolver el problema de optimización. El primer método determina los límites de la región de controladores PID estabilizantes utilizando el método propuesto por Hohenbichler y Abel (2008). Posteriormente genera la población inicial y utiliza operadores especiales de variación para mantenerse en el interior de la región factible. El segundo método determina una aproximación de los límites resolviendo previamente un primer sub-problema: se intenta maximizar simultaneamente el volumen y el número de individuos factibles muestreados aleatoria y uniformemente. A continuación se selecciona un conjunto de vértices en particular con el fin de continuar el proceso de búsqueda considerando los objetivos originales y enfocándose en la región delimitada por el conjunto de vértices seleccionados. Las pruebas estadísticas muestran que, en el caso del problema de diseño bajo estudio, ambos métodos permiten obtener mejores valores de los indicadores de desempeño con respecto al método tradicional.

CAPÍTULO V BÚSQUEDA MULTI-OBJETIVO LOCAL

En este capítulo se presenta un nuevo operador genético específicamente diseñado para la búsqueda local multi-objetivo, llamado Ascenso de Colina con Paso Lateral (ACPL). En caso de no encontrar soluciones dominantes el operador utiliza la información recolectada para intentar moverse hacia algún cono de diversidad (paso lateral). Se propone un esquema para la integración del ACPL con el algoritmo SPEA2, generando de esta forma un algoritmo híbrido (memético) denominado SPEA2-ACPL. Se presentan los resultados obtenidos cuando se utiliza el operador propuesto para resolver problemas que involucran funciones convexas y no-convexas. Parte de los resultados presentados en este capítulo fueron previamente incluidos como parte de las publicaciones (Schütze et al., 2008) y (Lara et al., 2010).

5.1. Antecedentes Algunos algoritmos de optimización se caracterizan por lograr buenas propiedades de convergencia local, lo cual viene asociado al costo de la información sobre las derivadas de la función objetivo. Por otra parte, los algoritmos genéticos (y muchos otros similares) permiten explorar el espacio de búsqueda para intentar aproximarse a los mínimos globales. Los algoritmos meméticos (Moscato, 1989) intentan combinar las ventajas de un algoritmo genético y un algoritmo de búsqueda local. Una de las primeras referencias en el caso multi-objetivo fue presentada por Ishibuchi & Murata (1996) con el nombre de Multi-Objective Genetic Local Search (MOGLS). Este autor propuso utilizar una función escalar de ponderación con pesos aleatorios para seleccionar los padres a partir de los cuales se debe realizar la búsqueda local. Otro memético importante desde el punto de vista histórico es el M-PAES propuesto en (Knowles & Corne, 2000). Este algoritmo utiliza dos archivos de almacenamiento:

83 uno para las soluciones no-dominadas globales (que funciona de la misma manera que en el SPEA2 por ejemplo) y el otro para soluciones no-dominadas locales que se emplea únicamente para decidir si se acepta el individuo producido por la búsqueda local. En (Murata et al., 2003) se propuso un operador que utiliza una regla generalizada de dominancia, basada en el número de objetivos que se ven mejorados al aplicar el operador. Típicamente la estrategia más utilizada consiste en reemplazar el individuo actual por otro que lo domina en todos los objetivos. Una de las primeras referencias en el caso de funciones definidas para un rango continuo es el artículo (Goel & Deb, 2002), en el que un operador de búsqueda local fue combinado con el NSGA-II. En las primeras pruebas el operador se aplicó al finalizar la ejecución del NSGA-II (de manera que realmente no hay interacción entre el operador de búsqueda local y los operadores genéticos). Posteriormente se aplicó el operador de búsqueda local en cada generación. Los autores encontraron que el operador local mejoraba la calidad de las soluciones obtenidas, pero aumentando considerablemente el número de evaluaciones del vector objetivo. En (Sharma et al., 2007) se combinó una búsqueda local de tipo Sequential Quadratic Programming (SQP) con los algoritmos NSGA-II y SPEA2 para resolver un conjunto de problemas benchmark. Los autores afirman que si no hay frentes de Pareto locales, el algoritmo híbrido converge más rápido hacia un COGP sin perder diversidad en las soluciones. En (Adra et al., 2005b) tres métodos de búsqueda local fueron hibridizados con MOGA: Recocido Simulado, Ascenso de Colina y Búsqueda Tabú. Los tres híbridos fueron comparados contra el MOGA tomando en cuenta el mismo número de evaluaciones. Los autores afirman que el híbrido MOGA - Ascenso de Colina produjo los mejores resultados. También observaron que el proceso de búsqueda local puede provocar el fenómeno de convergencia prematura si no se incorpora un mecanismo de preservación de la diversidad. En (Wanner et al., 2006) el operador de búsqueda local utiliza aproximaciones cuadráticas de todas las funciones objetivo, utilizando la información obtenida a partir del punto seleccionado. Se presentan los resultados obtenidos mediante la hibridización con el algoritmo SPEA2.

84

5.2. Ascenso de Colina con Paso Lateral (ACPL) Sea xk ∈ X a un punto factible a partir del cual se desea llevar a cabo el proceso de búsqueda local y xk+1 el nuevo punto generado por el operador. En el caso multi-objetivo sería deseable que esta operación verifique las siguientes propiedades: 1. Si xk se encuentra “lejos” de un COLP, entonces xk+1 debería ser tal que xk+1 xk (descenso en todos los objetivos). 2. En caso contrario, xk+1 debería ser tal que xk+1  xk , lo cual puede contribuir a la exploración del COLP. 3. Por otra parte el operador debe ser capaz de utilizar información de gradiente, tomar en cuenta las restricciones del problema e interactuar fácilmente con un algoritmo de búsqueda global. En esta sección se presenta el operador ACPL, especialmente diseñado para problemas multi-objetivo intentando cumplir en la medida de lo posible con las propiedades antes enumeradas. En caso de no encontrar soluciones dominantes, el operador utiliza la información recolectada para intentar generar soluciones que pertenezcan a los conos de diversidad (paso lateral). La idea principal se inspira en el análisis geométrico realizado por Brown y Smith (2005). Cuando un punto se encuentra “lejos” de un COLP, los conos de ascenso, descenso y diversidad presentan aproximadamente el mismo ángulo (ver figura 5.1). Si por el contrario, el punto se encuentra “cerca” de un COLP, los ángulos de ascenso y descenso se hacen menores (ver figura 5.2) y como consecuencia es menos probable seleccionar al azar una dirección en estos conos.

5.2.1

Caso sin información de gradientes

A continuación se describe el funcionamiento del operador, asumiendo que no se dispone de la información del gradiente de ninguna función objetivo (ver algoritmo 5.1). Sea x0 ∈ X (con X ⊂ Rn ) un punto factible a partir del cual se comienza la búsqueda local con el objetivo de generar un conjunto Xnd de nuevos puntos factibles no dominados con respecto a x0.

85

Figura 5.1: Situación que se presenta cuando el punto de prueba se encuentra lejos de un COLP. Cono de descenso amplio (región sombreada).

Figura 5.2: Situación que se presenta cuando el punto de prueba se encuentra cerca de un COLP. Cono de descenso estrecho (región sombreada).

86

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

Entradas : x0, r, Nnd, ε,Niter Salidas : Xn k = 1; ab = 0, ∀b ∈ BM ; nb = 0, ∀b ∈ BM ; Mientras k ≤ Niter hacer Seleccionar x1 ∈B(x0 , r); Si f (x1) f (x0 ) entonces d = x1 − x0 ; Calcular t de acuerdo a (5.4); x2 = x0 + td; Xnd = Archivar(x2 ); ab = 0, ∀b ∈ BM ; nb = 0, ∀b ∈ BM ; x0 = x2 ; sino si f (x0 ) f (x1 ) entonces d = x0 − x1 ; Calcular t de acuerdo a (5.4); x2 = x0 + td; Xnd = Archivar(x2 ); ab = 0, ∀b ∈ BM ; nb = 0, ∀b ∈ BM ; x0 = x2 ; sino b = tipo de cono de diversidad de x1 con respecto a x0 ; nb = nb + 1; 1 x0 ab = ab + xx1 − −x0  ; finsi k = k + 1; Fin k = 1; Mientras k ≤ M hacer Si nb = r, Nnd b entonces Calcular tb de acuerdo a (5.10); x2 = x0 + tb ab ); Xn = Archivar(x2) finsi k = k + 1; Fin Algoritmo 5.1 : Operador ACPL

87

En primer lugar se selecciona de manera aleatoria un punto de prueba x1 perteneciente al vecindario B(x0, r) definido por   B(x0 , r) = x ∈ X : x0i − ri ≤ xi ≤ x0i + ri , i = 1, 2, . . . , n

(5.1)

En este punto existen varias posibilidades. Si se cumple f (x1) f (x0) entonces se define la dirección de búsqueda como (5.2) d = x1 − x0 Si por el contrario se cumple la condición f (x0 ) f (x1 ) entonces es posible realizar la búsqueda en la dirección inversa d = x0 − x1. El siguiente paso consiste en llevar a cabo una búsqueda lineal multi-objetivo en la dirección d, resolviendo M sub-problemas de la forma min fj (x0 + tj d), tj ∈R

j = 1, 2, . . . , M

(5.3)

para obtener un conjunto de soluciones {t∗1, t∗2 , . . . , t∗M }, a partir de las cuales se debe seleccionar la longitud del paso t (Dennis & Schnabel, 1983). Para intentar garantizar el descenso de todas las funciones objetivo generalmente la estrategia más conservadora consiste en seleccionar t = min {t∗1, t∗2 , . . . , t∗M } (5.4) A continuación se describe un método aproximado con el fin de resolver los subproblemas escalares (5.3). Sea f : R → R la función que se desea minimizar. Como datos de entrada se tienen dos valores t0 < t1 tales que f (t0) < f (t1 ). Defina ∆ = t1 − t0

(5.5)

tl = t0 t m = t1 tr = t0 + 2∆ En primer lugar se comprueba la condición de curvatura f (tr ) − f (tm) f (tm ) − f (tl ) < tm − t l tr − tm

(5.6)

Si esta condición se cumple, entonces es posible aproximar a f por un polinomio de la forma p(t) = at2 + bt + c (5.7)

88 cuyos coeficientes a, b, y c pueden ser calculados a partir de los valores p(tl ) = f (tl ), p(tm ) = f (tm) y p(tr ) = f (tr ). En particular, el mínimo de p(t) se alcanza para b (5.8) 2a cantidad que puede servir como estimado de la solución del problema (5.3) para cada función objetivo. t∗ = −

Si la condición de curvatura (5.6) no se cumple, entonces se puede conjeturar que todavía es posible dar un paso más grande fijando tl = tm, tm = tr y tr = 2tr . Este procedimiento puede realizarse tantas veces como sea necesario, aumentando progresivamente la longitud del paso hasta que se cumpla la condición (5.6) o hasta que se llegue a un punto fuera de la región factible (ver sección de tratamiento de las restricciones más adelante). Suponga ahora que x0 y x1 resulten incomparables. El operador ACPL identifica el tipo de cono de diversidad al cual pertenece x1 con respecto a x0 y almacena esta información. Seguidamente se selecciona un nuevo punto de prueba x2 en B(x0 , r) y se repite el procedimiento descrito anteriormente para el procesamiento de x1 .

0 Si durante este proceso se recolectan nb ≥ Nnd b puntos pertenecientes a C(x , b) entonces se intenta realizar movimientos laterales utilizando las direcciones acumuladas: nb  1 xj − x0 b a = (5.9) nb xj − x0 0 j ∈{índices de los puntos pertenecientes a C (x ,b)}

donde los xj que intervienen en la expresión anterior corresponden a la secuencia de puntos pertenecientes a C(x0 , b). El vector Nnd debe ser definido por el usuario, en función del presupuesto disponible para el número de evaluaciones de las funciones objetivo.

La longitud del paso lateral tb en la dirección ab puede ser calculada mediante la ecuación εb (5.10) tb = b a ∞ donde ε es un vector de parámetros que regulan la distancia entre dos soluciones vecinas del conjunto de Pareto, para cada cono de diversidad. Todo el proceso descrito en los párrafos anteriores se ejecuta la cantidad de Niter veces: mientras mayor sea este parámetro, mayor número de puntos serán muestreados en el vecindario del punto de origen.

89

5.2.2

Caso con información de gradientes

Cuando se tiene acceso a todos los gradientes de las funciones objetivo, es posible calcular la dirección de descenso propuesta en el capítulo 3 (ver ecuación 3.17). Posteriormente se puede calcular la longitud del paso mediante el método que se describió en la sección anterior o preferiblemente mediante una búsqueda inexacta intentando que se satisfagan las condiciones de Wolfe para todas las funciones objetivo (Dennis & Schnabel, 1983). Con respecto al movimiento lateral, en cuanto se detecta la condición 2 M     ∗ ∗ αi ∇fi (x ) ≤ εP    i=1

2

(5.11)

donde εP > 0 es una cantidad pequeña a definir, se puede conjeturar que el punto actual se encuentra cerca de un COLP. En este sentido, un método para intentar obtener puntos incomparables pertenecientes a un vecindario de x∗ fue propuesto en (Hillermeier, 2001). En primer lugar se determinan los vectores qi ∈ Rn , i = 1, 2, ..., M − 1 extraídos de las primeras n filas de las M − 1 columnas que conforman la matriz QK ∈ R(M +n)×(M −1), obtenida a su vez a partir de la factorización ortogonal QR " ! ′ ∗ ∗ T (5.12) QN QK ∗ R = (M (x , α ))

donde M ′ ∈ R(M +n)×(M +n) es la matriz jacobiana de    M α ∇f (x) i=1 i i  M(x, α) =   M i=1 αi − 1

(5.13)

Note que los vectores qi forman una base ortonormal de un hiperplano de dimensión M − 1 tangente al COLP en x∗ y por lo tanto pueden ser utilizados (así como cualquier combinación lineal) como potencial dirección de diversidad. La longitud del paso lateral se determina como εb (5.14) tb = b q ∞ donde ε es el vector que se introdujo previamente en la ecuación (5.10).

90

5.2.3

Tratamiento de restricciones

Para el manejo de las restricciones, en el caso en que no se maneja información de gradientes, se propone adoptar un mecanismo de retorno a la región factible. Suponiendo que el punto de origen es factible, se realiza una búsqueda lineal multi-objetivo en la dirección de descenso, tal como se explicó en las secciones anteriores. Si el resultado no es factible, se ejecuta el algoritmo 5.2, basado en el método de bisección hasta alcanzar la precisión deseada, especificada mediante el parámetro T ol.

1 2 3 4 5 6 7 8 9 10 11 12 13

Entradas : x0 ∈ X , x1 ∈ / X , T ol Salidas : x # ∈ x0x1 ∩ X con miny∈∂ X y − x # ≤ T ol in0 = x0; out0 = x1 ; i = 0; Mientras outi − ini  ≥ T ol hacer mi = ini + 12 (outi − ini ); Si mi ∈ X entonces ini+1 = mi ; outi+1 = outi ; sino ini+1 = ini ; outi+1 = mi ; i = i + 1; Fin Algoritmo 5.2 : Retorno a la región factible para el operador ACPL1

Con respecto al caso en que se maneja información de gradiente, al momento en que se redacta este documento todavía no se cuenta con una versión del operador que tome en cuenta el gradiente de las funciones de restricción de manera eficiente. No obstante, para las restricciones de intervalo sobre las variables, del tipo xi min ≤ xi ≤ xi max

(5.15)

se propone utilizar un mecanismo sencillo de proyección a la región factible de la forma    x si xi min ≤ xi ≤ xi max   i xpi = (5.16) xi min si xi min > xi     x i max si xi max < xi donde xp es el vector proyectado. Este método tiene como ventaja que no se necesita especificar ningún parámetro adicional desde el punto de vista del funcionamiento del operador.

91 En la tabla 5.1 se muestra un resumen de los parámetros a ser fijados para el correcto funcionamiento del operador, en los casos con o sin información de gradiente. Tabla 5.1: Parámetros que deben ser definidos por el usuario en los casos con y sin información de gradientes Parámetro Sin ∇fi (x) Con ∇fi (x) r ∈Rn+ Si No M −2 2 nd Si No N ∈N+ M −2 2 ε ∈R+ Si Si Niter ∈N+ Si Si εP ∈R+ No Si T ol ∈ R+ Si Si

5.3. Evaluación del operador ACPL 5.3.1

Problema convexo sin restricciones En esta sección considere el siguiente problema   4 2 (x1 − 1) + (x2 − 1)  P1 :min2  x∈R (x1 + 1)2 + (x2 + 1)2

(5.17)

Las figuras 5.3 y 5.4 muestran el verdadero conjunto y frente de Pareto, obtenidos resolviendo el problema escalar $ % $ % P1′ :min2 α (x1 − 1)4 + (x2 − 1)2 + (1 − α) (x1 + 1)2 + (x2 + 1)2 x∈R

para 1000 valores de α distribuidos uniformemente entre 0 y 1.

A continuación se presentan los resultados obtenidos al aplicar el operador ACPL sucesivamente a partir de un punto de origen. Para verificar si el operador es capaz de generar todo el frente, la prueba se divide en tres etapas. En primer lugar se aplica el operador hasta que se detecta el primer paso lateral. Inmediatamente se modifica el funcionamiento del operador para aceptar únicamente las soluciones que mejoren un objetivo, por ejemplo f1. Durante la tercera etapa se aceptan únicamente las soluciones que mejoren el otro objetivo, por ejemplo f2 . Se realizaron pruebas utilizando las dos modalidades del operador: sin información de gradiente (denotado ACPL1) y con información de gradiente (denotado ACPL2). Para la

92 1 0.8 0.6 0.4

x

2

0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1

−0.5

0 x

0.5

1

1

Figura 5.3: Verdadero conjunto de Pareto para el problema P1 prueba de ambas versiones se utilizó el mismo punto de inicio x0 = [5, 5]. Los valores de los parámetros característicos del operador se muestran en la tabla 5.2. Tabla 5.2: Valores utilizados para los parámetros ACPL Parámetro Valor ε [0.1, 0.1] εP 0.01 r [0.01, 0.01] Niter 10 nd N [5, 5] T ol 10−4

En las figuras 5.5 y 5.6 se muestra un ejemplo de apróximación generada por ACPL1, en el espacio de objetivos y en el espacio de variables respectivamente. En total, ACPL1 requirió 925 evaluaciones de las funciones objetivo y ACPL2 únicamente 143. Sin embargo ACPL2 empleó 50 evaluaciones de la matriz jacobiana y 50 evaluaciones de las matrices hessianas. Note que la aproximación producida por ACPL2 es prácticamente idéntica al frente de Pareto generado mediante escalarización, tal como es posible apreciar en las figuras 5.7 y 5.8.

93

8 7 6

f

2

5 4 3 2 1 0

0

5

10 f1

15

20

Figura 5.4: Verdadero frente de Pareto para el problema P1

1.5

1

Verdadero Conjunto de Pareto Aproximación obtenida mediante ACPL1

x2

0.5

0

−0.5

−1

−1.5 −1

−0.5

0 x1

0.5

1

Figura 5.5: Aproximación del conjunto de Pareto generada por ACPL1. Problema P 1

94

7

Verdadero frente de Pareto Aproximación obtenida mediante ACPL1

6

f

2

5 4 3 2 1 0 0

2

4

6

8

10

12

14

16

18

f

1

Figura 5.6: Aproximación del frente de Pareto generada por ACPL1. Problema P 1

1 0.8 0.6

Verdadero Conjunto de Pareto

0.4

Aproximación obtenida mediante ACPL2

x2

0.2 0 −0.2 −0.4 −0.6 −0.8 −1

−0.5

0 x1

0.5

1

Figura 5.7: Aproximación del conjunto de Pareto generada por ACPL2. Problema P 1

95 8 Verdadero frente de Pareto Aproximación obtenida mediante ACPL2

7 6

f

2

5 4 3 2 1 0

0

5

10

15

20

25

f1

Figura 5.8: Aproximación del frente de Pareto generada por ACPL2. Problema P 1

5.3.2

Problema convexo con restricciones de intervalo En esta sección considere el siguiente problema   4 2 (x − 1) + (x2 − 1)  1  P2 : min x∈[0.5,1]×[1,2] (x1 + 1)2 + (x2 + 1)2

(5.18)

Las figuras 5.9 y 5.10 muestran el verdadero conjunto y frente de Pareto, obtenidos resolviendo el problema escalar P2′ :

% $ % $ min α (x1 − 1)4 + (x2 − 1)2 + (1 − α) (x1 + 1)2 + (x2 + 1)2 x∈[0.5,1]×[1,2]

para 1000 valores de α distribuidos uniformemente entre 0 y 1.

Para la prueba de ambas versiones del algoritmo se utilizó el mismo punto de inicio factible x0 = [1, 2] con los mismos valores de los parámetros mostrados en la tabla 5.2. En las figuras 5.11 y 5.12 se puede observar un ejemplo del tipo de apróximación que puede generar ACPL1 en el caso con restricciones. En las figuras 5.13 y 5.14 se observan las correspondientes a ACPL2. En total, ACPL1 utilizó 2021 evaluaciones de las funciones objetivo. ACPL2 requirió 103 evaluaciones de las funciones objetivo, 50 evaluaciones de la matriz jacobiana y 50 evaluaciones de las matrices hessianas.

96

1.25 1.2 1.15 1.1 1.05 x

2

1 0.95 0.9 0.85 0.8 0.75 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

x1

Figura 5.9: Verdadero conjunto de Pareto para el problema P2

8.2 8 7.8 7.6

f

2

7.4 7.2 7 6.8 6.6 6.4 0

0.01

0.02

0.03

0.04

0.05

0.06

f1

Figura 5.10: Verdadero frente de Pareto para el problema P2

97

1.03 Verdadero Conjunto de Pareto Aproximación obtenida mediante ACPL1

1.02

x

2

1.01

1

0.99

0.98

0.97 0.5

0.6

0.7

0.8 x

0.9

1

1

Figura 5.11: Aproximación del conjunto de Pareto generada por ACPL1. Problema P 2

8.2 Verdadero frente de Pareto Aproximación obtenida mediante ACPL1

8 7.8 7.6

f2

7.4 7.2 7 6.8 6.6 6.4 0

0.01

0.02

0.03 f1

0.04

0.05

Figura 5.12: Aproximación del frente de Pareto generada por ACPL1. Problema P 2

98

1.01

Verdadero Conjunto de Pareto Aproximación generada popr ACPL2

1.008 1.006 1.004

x

2

1.002 1 0.998 0.996 0.994 0.992 0.99 0.5

0.6

0.7

0.8

0.9

1

x1

Figura 5.13: Aproximación del conjunto de Pareto generada por ACPL2. Problema P 2

8.2 Verdadero Frente de Pareto Aproximación generada por ACPL2

8 7.8 7.6

f

2

7.4 7.2 7 6.8 6.6 6.4 0

0.01

0.02

0.03

0.04

0.05

0.06

f1

Figura 5.14: Aproximación del frente de Pareto generada por ACPL2. Problema P 2

99

5.3.3

Problema no convexo con restricciones de intervalo

En esta sección considere el siguiente problema no convexo propuesto en (Zitzler et al., 1999)   x  1 & (  ' min (5.19) P3 :   x 1 g(x) 1 − g(x) xi ∈ [−5, 5] , i = 2, 3, . . . , 20 x1 ∈ [0, 1]

con g(x) = 191 +

20   i=2

x2i − 10 cos(4πxi )



Este problema posee 2119 frentes locales y el frente global, para el cual g(x) = 1, se obtiene con xi = 0, i = 2, 3, . . . , 20. Se realizaron pruebas con ambas versiones del operador, iniciando la búsqueda a partir del punto   x0 = 0.5 3 3 . . . 3 (5.20)

En la figura 5.15 se muestra un ejemplo de las aproximaciones que pueden obtenerse mediante cada versión del operador.

Aproximación ACPL2 Verdadero frente de Pareto Aproximación ACPL1

250

200

150

100

50

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Figura 5.15: Ejemplo de aproximaciones del frente de Pareto generadas por ACPL1 y ACPL2. Problema P 3

100 Los valores de los parámetros empleados se muestran en la tabla 5.3. Como era de esperar, ambas versiones del operador producen aproximaciones que se encuentran completamente alejadas del verdadero frente, esto motivado por el gran número de frentes locales característicos de este problema. Tabla 5.3: Valores utilizados para los parámetros ACPL. Problema P3 Parámetro Valor ε [0.05, 0.05] 1 εP   0.01 0.1 0.1 . . . 0.1 r Niter 20 nd N [5, 5] T ol 10−4

5.4. Sensibilidad del ACPL El comportamiento del operador ACPL depende de los valores seleccionados para sus parámetros. En esta sección se considera de nuevo el problema P1 formulado en la sección anterior. La influencia del vector r ∈Rn+ se manifiesta en el número de evaluaciones de las funciones objetivo necesarias para alcanzar el Frente de Pareto. En la tabla 5.4 se presenta un ejemplo del número de evaluaciones que emplea ACPL1 antes de que se produzca el primer intento de paso lateral. Note que este número disminuye cuando aumenta r. Mientras menor sea esta norma, menor será el tamaño de los pasos que ejecuta el operador y en consecuencia se requieren más pasos para cubrir la misma distancia. Cuando r aumenta a partir de cierto punto, el operador se hace demasiado impreciso, de manera que se hace muy difícil conseguir direcciones de descenso. En consecuencia, el número de evaluaciones de las funciones objetivo necesarios para alcanzar el frente de Pareto aumenta de manera exponencial.

101 Tabla 5.4: Sensibilidad del operador ACPL1 ante variaciones en el vector r r Evaluaciones de las funciones objetivo [0.001, 0.001] 895 [0.01, 0.01] 135 [0.1, 0.1] 74 [1, 1] 7269

2M −2 regula la distancia entre dos soluComo se ha visto previamente, el vector ε ∈R+ ciones vecinas del conjunto de Pareto, para cada cono de diversidad. En las figuras 5.16, 5.17 y 5.18 se muestran las aproximaciones obtenidas por ACPL2 a partir de tres valores del vector ε. En las mismas es posible apreciar que el espaciamiento entre las soluciones 2M −2. es regulado mediante los valores especificados en el vector ε ∈R+ El valor del parámetro Niter debe ser ajustado en función del número de evaluaciones que es posible invertir durante la ejecución del operador. Mientras más grande sea Niter , mayor cantidad de puntos vecinos serán muestreados en busca de direcciones de descenso o de diversidad. Como regla general, esta cantidad debe ser mayor que la mayor de las cantidades especificadas en el vector Nnd, para permitir que pueda ocurrir al menos un paso lateral en cualquier cono de diversidad. Por su parte, la influencia del vector Nnd se manifiesta en el número de pasos laterales ejecutados por ACPL1 para intentar cubrir el frente de Pareto. En la tabla 5.5 se presenta el número de pasos laterales que son efectivamente ejecutados para diferentes valores del vector Nnd , fijando Niter = 20. 1.5 Verdadero conunto de Pareto Aproximación generada por ACPL2

1

x

2

0.5 0 −0.5 −1 −1.5 −1

−0.8

−0.6

−0.4

−0.2

0 x

0.2

0.4

0.6

0.8

1

Figura 5.16: Aproximación obtenida mediante ACPL2 con ε1=[0.01, 0.01]

102

1

Verdadero conunto de Pareto Aproximación generada por ACPL2

x

2

0.5

0

−0.5

−1 −1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

x1

Figura 5.17: Aproximación obtenida mediante ACPL2 con ε2=[0.1, 0.1]

1 0.8 0.6

Verdadero conjunto de Pareto Aproximación obtenida mediante ACPL2

0.4

x2

0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1

−0.5

0 x1

0.5

1

Figura 5.18: Aproximación obtenida mediante ACPL2 con ε3=[1, 1]

103 Note que mientras menos pasos laterales se ejecutan, menos evaluaciones de las funciones objetivo se invierten en el funcionamiento del operador y la aproximación de la dirección de diversidad es más precisa (se promedian más vectores). La contraparte de esto es que la longitud del frente que se cubre es menor. En efecto, note que en el caso extremo, cuando Nnd = [20, 20], no se ejecuta ningún paso lateral. Tabla 5.5: Sensibilidad del operador ACPL1 ante variaciones en el vector Nnd Nnd Número de pasos laterales N◦ de evaluaciones [2, 2] 202 25024 [5, 5] 45 5591 [10, 10] 30 3919 [20, 20] 0 2000

Las figuras 5.19, 5.20 y 5.21 permiten observar la influencia del parámetro εP en el funcionamiento del operador ACPL2 (la versión que maneja información de gradientes). En general se aprecia que mientras menor sea el valor de εP , más precisa será la aproximación del Conjunto de Pareto, pero por otra parte se invertirá un mayor número de evaluaciones de las funciones objetivo para generarla.

1 0.8 0.6 0.4

Verdadero Conjunto de Pareto Aproximación generada por ACPL2

x2

0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1

−0.5

0 x1

0.5

1

Figura 5.19: Aproximación obtenida mediante ACPL2 con εP = 0.01

104

1 Verdadero Conjunto de Pareto Aproximación generada por ACPL2

0.8 0.6 0.4

x

2

0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1

−0.5

0 x

0.5

1

1

Figura 5.20: Aproximación obtenida mediante ACPL2 con εP = 0.1

1 Verdadero Conjunto de Pareto Aproximación generada por ACPL2

0.8 0.6 0.4

x

2

0.2 0 −0.2 −0.4 −0.6 −0.8 −1

−0.5

0 x

0.5

1

1

Figura 5.21: Aproximación obtenida mediante ACPL2 con εP = 1

105 Para finalizar esta sección, considere el problema de optimización con restricciones P2. La influencia del parámetro T ol se manifiesta en la distancia tolerada para las soluciones de la aproximación con respecto con respecto al borde de la región factible (ver figura 5.22). 1.08 Aproximación obtenida con ACPL1 para Tol = 0.01 Aproximación obtenida con ACPL1 para Tol = 0.1

1.07 1.06 1.05 1.04 1.03 1.02 1.01 1 0.99 0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

Figura 5.22: Aproximaciones obtenidas por ACPL1 para distintos valores del parámetro T ol

5.5. Integración SPEA2-ACPL Tal como se expuso en la sección anterior, los algoritmos meméticos intentan combinar las ventajas de un algoritmo genético con las de un algoritmo de búsqueda local. Sin embargo, no existe en la actualidad un consenso con respecto a la mejor forma de realizar esta integración. Entre los elementos que deben tomarse en cuenta se encuentran: • A cuántos y a cuales individuos debe aplicarse el operador? • En qué momento? Con qué frecuencia? Durante cuántas generaciones? De acuerdo con la investigación que se ha realizado, estas preguntas no han sido respondidas de manera satisfactoria. El compromiso entre búsqueda local vs búsqueda global es similar al que existe entre individualismo y colectivismo. En efecto, cuando es necesario asignar recursos limitados a una población, se debe decidir entre si se quiere lograr mejorar mucho la condición de pocos individuos, o bien mejorar poco la condición de muchos. En este trabajo, se propone utilizar el operador ACPL como un operador de mutación especial, el cual se aplica únicamente a los individuos pertenecientes al archivo externo

106 actualizado por SPEA2. En lo que sigue se denota SPEA2-ACPL1 a la versión híbrida correspondiente al SPEA2 en conjunto con ACPL1 (operador que no utiliza información de gradiente). De forma similar, se denota SPEA2-ACPL2 a la versión híbrida correspondiente al SPEA2 en conjunto con ACPL2 (operador que utiliza información de gradiente). En particular se propone evaluar tres estrategias de integración: • E1: Aplicar el operador en todas las generaciones con probabilidad pACPL a partir de la primera generación. • E2: Aplicar el operador de manera intermitente, una de cada NACPL generaciones, con probabilidad pACP L . • E3: Aplicar el operador con probabilidad pACP L cuando se haya consumido el 50% de la cantidad total de evaluaciones disponibles.

Para todas las pruebas que se presentan a continuación se consideró un presupuesto fijo de 20.000 evaluaciones del vector objetivo. Por una parte esta cantidad de evaluaciones permite que las 30 ejecuciones del algoritmo genético no excedan el tiempo máximo permitido (≈24 horas). Por otra parte esta cantidad de evaluaciones es suficiente para alcanzar una etapa del proceso en la cual las aproximaciones de los frentes evolucionan lentamente con respecto al inicio.

5.5.1

Problema convexo sin restricciones En esta sección considere el siguiente problema convexo   30  4 2  (x1 − 1) + (xi − 1)  i=2  P4 : min30  30    x∈R 2 (xi + 1)

(5.21)

i=1

el cual involucra 30 variables de decisión. El verdadero frente de Pareto puede ser obtenido resolviendo el problema escalar ) * ) 30 * 30   (xi − 1)2 + (1 − α) (xi + 1)2 (5.22) P4′ : min30 α (x1 − 1)4 + x∈R i=2 i=1 con 100 valores de α distribuidos uniformemente entre 0 y 1.

107 La configuración utilizada para la ejecución de las pruebas se muestra en la tabla 5.6. En la tabla 5.7 se presenta el valor medio y la desviación estándar calculada para los indicadores GD, IGD evaluados a partir de 30 ejecuciones de cada algoritmo. El valor resaltado con un asterisco corresponde al algoritmo que consigue estadísticamente los mejores valores, de acuerdo con la prueba de Wilcoxon, con un nivel de significancia α = 0.05. Note que las versiones híbridas permiten, cualquiera sea la estrategia de integración, obtener mejores aproximaciones considerando ambos indicadores. Por otra parte, la estrategia E1, la cual consiste en aplicar el operador a partir de la primera iteración, es la que consigue los mejores resultados. Por lo tanto es posible concluir que para este problema, la mejor opción consiste en utilizar la información de gradiente lo más pronto posible. Tabla 5.6: Configuración utilizada para resolver el problema P4 Parámetros Valor Representación Números reales Operador de cruce SBX Probabilidad de cruce 0.9 Operador de mutación Perturbación gaussiana Probabilidad de mutación 0.1 Probabilidad de aplicación del ACPL 0.5 10 NACPL ε [0.1, 0.1] εP 0.1 r [0.01, 0.01] Niter 20 Nnd [10, 3] T ol 10−4 Condición de parada 200 generaciones Número máximo de evaluaciones 20.000 Tamaño máximo del archivo externo 100 Esquema de selección (100 + 100)

108 Tabla 5.7: Indicadores GD/IGD obtenidos para el problema P4 GD IGD SPEA2 10.5563(5.4469) 16.1034(6.6455) E1 : 6.7728(1.8759)∗ E1 : 3.2782(0.3388)∗ E2 : 11.6404(3.3716) SPEA2-ACPL1 E2 : 5.8222(1.4305) E3 : 4.6148(0.5481) E3 : 11.0787(2.2915) E1 : 0.5350(0.0479)∗ E1 : 0.6797(0.0284)∗ E2 : 2.2375(2.7337) SPEA2-ACPL2 E2 : 0.7679(0.0722) ∗ E3 : 0.6645(0.0306) E3 : 0.9584(0.7102)

5.5.2

Problema no convexo con restricciones de intervalo

En esta sección se retoma el problema P3 anteriormente formulado. La configuración utilizada para estas pruebas se muestra en la tabla 5.8. En la tabla 5.9 se presenta el valor medio y la desviación estándar calculada para los indicadores GD, IGD evaluados a partir de 30 ejecuciones de cada algoritmo. Tabla 5.8: Configuración utilizada para las pruebas del problema P3 Parámetros Valor Representación Números reales Operador de cruce SBX Probabilidad de cruce 0.9 Operador de mutación Perturbación gaussiana Probabilidad de mutación 0.1 Probabilidad de aplicación del ACPL 0.1 NACPL 10 ε [0.1, 0.1] εP 0.1 r [0.01, 0.01] Niter 10 [5, 2] Nnd T ol 10−4 Condición de parada 200 generaciones Número máximo de evaluaciones 20.000 Tamaño máximo del archivo externo 100 Esquema de selección (100 + 100)

109 Tabla 5.9: Indicadores GD/IGD obtenidos para el problema P3 GD IGD SPEA2 3.6011(1.1629) 3.2091(0.9735) E1 : 3.7380(1.2520) E1 : 3.2947(1.0338) E2 : 3.3384(0.9472) SPEA2-ACPL1 E2 : 3.7443(1.1807) E3 : 3.7047(1.1083) E3 : 3.3347(0.8820) E1 : 3.1866(0.9524) E1 : 2.9530(0.8846) E2 : 3.3188(0.9771) SPEA2-ACPL2 E2 : 3.6439(1.0985) E3 : 3.3079(1.0496) E3 : 3.0206(0.9099) Las pruebas estadísticas indican que no existen diferencias significativas entre los resultados obtenidos por los diferentes algoritmos, independientemente de la estrategia de evolución que se adopte. Esto se explica por la gran cantidad de frentes locales, lo cual aumenta la probabilidad de que los operadores de búsqueda local queden atrapados. Ahora se considera una variación del problema P3 tal como se muestra a continuación   x  1 & (  ' (5.23) min P5 :   x g(x) 1 − g(x1 ) xi ∈ [−5, 5] , i = 2, 3, . . . , 20 x1 ∈ [0, 1]

con g(x) = 191 +

20   i=2

x2i − 10 cos(πxi )



Este problema posee 519 frentes locales y el frente global, para el cual g(x) = 1, se sigue obteniendo con xi = 0, i = 2, 3, . . . , 20. Los resultados se muestran en la tabla 5.10. Tabla 5.10: Indicadores GD/IGD obtenidos para el problema P5 GD IGD SPEA2 0.0401(0.0099) 0.0503(0.0062) E1 : 0.0443(0.0103) E1 : 0.0534(0.0063) SPEA2-ACPL1 E2 : 0.0430(0.0076) E2 : 0.0513(0.0052) E3 : 0.0426(0.0096) E3 : 0.0534(0.0054) ∗ E1 : 0.0146(0.0080) E1 : 0.0412(0.0035)∗ E2 : 0.0460(0.0052) SPEA2-ACPL2 E2 : 0.0272(0.0081) E3 : 0.0375(0.0102) E3 : 0.0509(0.0056)

110 Note que para esta variación el algoritmo SPEA2-ACPL2, utilizando la estrategia de integración E1, es nuevamente capaz de obtener mejores valores de los indicadores GD/IGD con respecto al SPEA2, a pesar de la gran cantidad de frentes locales.

5.6. Resumen del capítulo En este capítulo se presentó un nuevo operador local, denominado ACPL, específicamente diseñado para la búsqueda multi-objetivo. Las principales bondades del operador son las siguientes: 1. Si el punto de prueba se encuentra “lejos” de un COLP, entonces se intenta conseguir una dirección de descenso en todos los objetivos. 2. Si el punto de prueba se encuentra “cerca” de un COLP, entonces se intenta conseguir una dirección de diversidad, para intentar mejorar la aproximación del frente. 3. El operador es capaz de utilizar la información de gradiente, en caso de que esta se encuentre disponible. 4. El operador es capaz de tomar en en cuenta las restricciones del problema y de interactuar fácilmente con un algoritmo de búsqueda global. Se consideraron dos versiones del operador ACPL1 (la cual no utiliza información de gradiente) y ACPL2 (que si utiliza esta información). Se presentaron los resultados obtenidos cuando se utiliza el operador, en sus dos versiones, de manera aislada o integrada con un algoritmo de búsqueda global para resolver problemas convexos y no convexos. Se evaluaron tres estrategias de integración que consisten en • E1: Aplicar el operador en todas las generaciones con probabilidad pACPL a partir de la primera generación. • E2: Aplicar el operador de manera intermitente, una de cada NACPL generaciones, con probabilidad pACP L . • E3: Aplicar el operador con probabilidad pACP L cuando se haya consumido el 50% de la cantidad total de evaluaciones disponibles. En el caso del problema convexo P 4, ambas versiones del algoritmo híbrido son capaces de obtener mejores valores para los indicadores GD/IGD con respecto a los obtenidos por

111 el SPEA2 de manera aislada, cualquiera sea la estrategia de integración utilizada. No obstante la estrategia E1 es la que obtiene los mejores resultados. En el caso del problema P 3 con 2119 frentes locales, ninguna versión del algoritmo híbrido es capaz de mejorar los resultados obtenidos por SPEA2. En el caso del problema P 5 caracterizado por 519 frentes locales, únicamente la versión SPEA2-ACPL2, utilizando la estrategia E1, logra mejorar los resultados del SPEA2.

CAPÍTULO VI PROBLEMA H2 /H∞ El proceso de diseño de un sistema de control tiene como objetivo dotar a la interconexión de propiedades que la planta no posee de manera aislada. En el capítulo 2 se definieron algunos índices de desempeño que pueden utilizarse para expresar matemáticamente y verificar en la práctica el cumplimiento de dichas propiedades. Por ejemplo, para un sistema SISO, la norma H2 mide la energía de la señal de respuesta del sistema ante un impulso unitario, mientras que la norma H∞ se relaciona con la amplitud máxima de la misma señal. En este capítulo se propone un nuevo método de diseño llamado Multi-Objective Pole Placement with Evolutionary Algorithms (MOPPEA) y se presentan los resultados obtenidos cuando se aplica dicho método para resolver el problema mixto H2 /H∞ en el caso de cinco modelos extraídos de la biblioteca COMP le ib (Leibfritz, 2004). Como base del proceso de búsqueda se utilizan los algoritmos SPEA2 y SPEA2-ACPL (ver capítulo 4) y se comparan los resultados con respecto a los obtenidos mediante una secuencia de problemas de factibilidad LMI. Parte de los resultados presentados en este capítulo fueron previamente publicados en (Sánchez et al., 2007) y (Sánchez et al., 2008b).

6.1. Formulación Se asume que el modelo de la planta a controlar (ver figura 6.1) es el siguiente:    x˙ = Ax + Iw + Bu   (6.1) z1 = y = Cx     z = Ax

2

el cual fue propuesto originalmente por Herreros (2000) en su tesis de doctorado, para comparar los resultados obtenidos por un algoritmo genético con respecto a la solución calculada mediante LMIs. Una ventaja de este modelo es que únicamente necesita las

113 matrices A, B y C para estar completamente definido, sin perder por ello complejidad en su estructura.

Figura 6.1: Modelo lineal a lazo cerrado Con el fin de analizar el conflicto entre los objetivos de robustez y desempeño del lazo, el problema de diseño mixto PH2/H∞ es formulado en este capítulo como:   G1 (Kc)2   PH2 /H∞ : min (6.2) Kc ∈R(s)nu ×ny G2 (Kc)∞

donde G1 y G2 son las funciones de transferencia a lazo cerrado desde w hacia z1 y z2 respectivamente.

6.2. Antecedentes Uno de los primeros artículos en los que se propone un intento de solución para el problema que nos interesa en este capítulo fue publicado por Khargonekar y Rotea en 1991. Estos autores propusieron un enfoque basado en ecuaciones algebraicas de Riccati para el problema planteado de la forma min

K ∈Rnu ×nx

G1 (K)2

sujeto a

(6.3)

G2(K)∞ ≤ γ donde γ es un escalar positivo. Doyle et al (1994) formularon las condiciones suficientes y necesarias para la existencia de soluciones, proporcionando expresiones matriciales explícitas para la determinación de

114 los controladores, en el caso de un problema formulado de la forma min

K ∈Rnu ×ny

G1(K)2 + γ G2 (K)∞

(6.4)

donde γ es un escalar positivo. Posiblemente uno de los artículos que ha tenido más influencia en el campo de las aplicaciones de LMIs al problema H2 /H∞ fue publicado por Carsten W. Scherer en 1995, el cual de hecho fue previamente mencionado en el capítulo 2. Este artículo ha inspirado y servido de base para muchos otros (Clement, 2001; Hassibi et al, 1999; Hindi et al, 1998; Miyamoto 1997). En cuanto a las soluciones del problema H2 /H∞ basadas en algoritmos genéticos, diversos autores han publicado artículos relacionados con este tema. En su tesis doctoral, Herreros (2000) presenta un estudio sobre el diseño de controladores robustos multi-objetivo por medio de algoritmos genéticos. Propone el algoritmo MRCD, junto con sus respectivos operadores y demuestra de manera gráfica que sus resultados son mejores con respecto a los producidos mediante LMIs. Como ya se ha mencionado, fue precisamente Herreros quien propuso la estructura que se presenta en la figura 6.1. De hecho, propone tres diferentes representaciones del espacio de controladores, con el fin de estudiar la influencia en la convergencia: la forma canónica controlable, la forma canónica diagonal y una forma canónica balanceada. En sus conclusiones, afirma que la forma controlable es la más adecuada para la representación del espacio de controladores ya que tiene más posibilidades de explorar el espacio de búsqueda. Molina-Cristobal (2005) y Popov (2005) utilizan el algoritmo MOGA para reducir el orden de los controladores y mejorar las soluciones con respecto a las obtenidas mediante LMIs. El primero utiliza una representación de los controladores en tiempo discreto, en forma de ceros y polos: K(z) = K

(z + z1) (z + z2) . . . (z + znz ) (z + p1 ) (z + p2) . . . (z + pnp)

(6.5)

y el segundo una representación polinomial en tiempo continuo b0snz + b1 snz−1 + . . . + bnz−1s + bnz K(s) = K np s + a1snp−1 + . . . + anp−1 s + anp

(6.6)

115 Ambos autores señalan que estas representaciones son ineficaces para generar soluciones factibles, lo cual los obliga a aumentar el tamaño de las poblaciones y en consecuencia el número de evaluaciones de las funciones objetivo. Takahashi et al (2004) proponen un algoritmo genético que asegura la generación de aproximaciones consistentes del Conjunto de Pareto y lo aplican al problema de diseño H2 /H∞. Estos autores se ocupan del problema de la realimentación estática de estados y proponen una representación de los controladores de la forma   K = k1 k2 . . . kn (6.7) lo cual les permite utilizar como población inicial la aproximación del Conjunto de Pareto generada mediante LMIs.

6.3. Solución mediante LMIs Una posible solución del problema PH2 /H∞ consiste en utilizar los resultados del trabajo de Scherer (1995) para resolver el problema mixto como una secuencia de problemas de factibilidad H2 y H∞ , de la manera que se describe en el algoritmo 6.1. La idea básica consiste en disminuir linealmente el límite superior de una de las restricciones, aumentando al mismo tiempo el límite de la otra. Note que aunque en este trabajo se plantea una evolución lineal para el ajuste, bien pudiera utilizarse un polinomio de orden superior u otras funciones monótonas decrecientes y crecientes. En el caso que nos ocupa, antes de ejecutar el algoritmo 6.1 es necesario seleccionar los valores de los parámetros NH2 /H∞ ,∆γ 2, ∆γ ∞ , γ min y γ max. El paso más importante consiste en determinar un controlador Kc que satisfaga las condiciones G1 (Kc )2 ≤ γ 2

(6.8)

G2 (Kc)∞ ≤ γ ∞ lo cual puede lograrse resolviendo un conjunto de LMIs, tal como fue mostrado en el capítulo 2. Para resolver el conjunto de LMIs es necesario igualar las matrices de decisión correspondientes a cada norma, razón por la cual el propio Scherer (1995) reconoce que este enfoque genera diseños conservadores. Sin embargo argumenta que “el método posee méritos

116 Entradas : G1 : Función de transferencia entre w y z1 G2 : Función de transferencia entre w y z2 NH2 /H∞ : Número de iteraciones ∆γ 2 : Distancia deseada entre dos soluciones vecinas. Norma H2 ∆γ ∞ : Distancia deseada entre dos soluciones vecinas. Norma H∞ γ min : Mínimo valor deseado para la norma H2 γ max : Máximo valor permitido para la norma H∞ Salidas : P 1 k = 1; 2 Mientras k ≤ NH2/H∞ hacer 3 Calcule Kck solución del sub-problema   G2 (Kck ) ≤ γ min + k∆γ 2 2   G1 (Kck ) ≤ γ max − k∆γ ∞ ∞

4 Pk = ActualizarArchivo(Kck , Pk−1); 5 k = k + 1; 6 Fin

Algoritmo 6.1 : Generación de una aproximación del Frente de Pareto para el problema mixto H2 /H∞ mediante LMIs valiosos con respecto a las alternativas disponibles, puesto que conduce a problemas que pueden ser resueltos numéricamente de manera eficiente”. Dado que las soluciones obtenidas al resolver la secuencia de problemas de factibilidad LMI no son necesariamente consistentes con una aproximación del Frente de Pareto, se propone utilizar un mecanismo para archivar las soluciones no-dominadas (ver capítulo 3) lo cual permite introducir en el proceso de búsqueda cualquier requerimiento que se desee con respecto a la distribución o diversidad deseada en la aproximación del frente. Es importante resaltar que en la investigación bibliográfica que se ha realizado, no ha sido posible encontrar referencias previas del empleo de un algoritmo del tipo que aquí se propone para resolver el problema PH2 /H∞ .

6.4. Método MOPPEA Usualmente la representación genética del espacio de controladores se realiza en formato de función de transferencia s2 + x2 s + x3 Kc (s) = x1 3 s + x4s2 + x5s + x6

(6.9)

117 o de matrices de estado (forma controlable u otra)  x1 x2 x3     1 0 0  Ac Bc = Kc =   0 1 0  Cc Dc  x4 x5 x6

1



  0    0   0

(6.10)

donde x1, x2 , . . . , x6 ∈ R son las variables de decisión (Cesareo, 1997).

Sin embargo, también se han propuesto representaciones más elaboradas, con el objeto de hacer el proceso de búsqueda más eficiente. Por ejemplo, Sánchez y Ferrer (2004) propusieron una representación basada en la parametrización de controladores estabilizantes de Youla, la cual permite eliminar la restricción de estabilidad en la formulación del problema de optimización. Liu et al (2002) incluyen en la representación genética la información correspondiente a los auto-valores y auto-vectores de la matriz de transferencia a lazo cerrado, de la forma p = (λ, R, L)

(6.11)

donde λ es el vector de auto-valores, R es la matriz de auto-vectores derechos y L es la matriz de auto-vectores izquierdos de la matriz de estados a lazo cerrado. Zavala et al (2002) proponen utilizar como variables de decisión la posición de las raices de los polinomios R(s), S(s) que intervienen en la ecuación diofantina característica de un lazo, de la forma F (s) = A(s)S(s) + B(s)R(s) (6.12) Existen aplicaciones que exigen representaciones muy particulares. Por ejemplo, en (Andres et al., 2004) se propuso un algoritmo llamado Variable Interval Multi-Objective (VIM) en el cual se codifica directamente la duración y la amplitud de los intervalos de la señal de control. También en (Villasana & Ochoa, 2003) se utiliza una representación de la señal de control basada en intervalos de duración variable. A continuación se propone un método de diseño, llamado Multi-Objective Pole Placement with Evolutionary Algorithms (MOPPEA) cuya idea principal consiste en asignar

118 una estructura de “Observación + Realimentación de Estados” (ver capítulo 2) a cada individuo, de la forma:  ·  x = (A + BK + LC) x − Ly (6.13)  u = K x

donde x  es la estimación del vector de estado, K es la matriz de realimentación de estados y L la matriz de observación. De manera que el problema de diseño PH2 /H∞ puede replantearse como:   G1 (K, L)2  min n×n  K ∈Rnu ×n ,L∈R y G2 (K, L)∞

(6.14)

y por lo tanto la restricción de estabilidad desaparece del problema de optimización: esto constituye la principal ventaja del método propuesto con respecto a los trabajos de Herreros (2000), Molina-Cristobal (2005) y Takahashi et al (2004). La representación genética propuesta consiste en concatenar las filas de las matrices K y L correspondientes a cada individuo en un vector de valores reales. De esta manera, el número total de variables de decisión se calcula mediante la fórmula nMOP PEA = n × (nu + ny )

(6.15)

donde n es el grado, nu el número de entradas y ny el número de salidas del sistema a lazo abierto (ver figura 6.2).

Figura 6.2: Representación genética propuesta. Método MOPPEA

Esta representación genética permite utilizar los operadores de variación descritos en el capítulo 3 y el operador de búsqueda local ACPL propuesto en el capítulo 4. Se propone en particular utilizar el operador de cruce aritmético. Si los padres seleccionados para el cruce pertenecen a la región de factibilidad, los hijos producidos por este operador son igualmente factibles.

119 Sean (K1, L1 ) y (K2 , L2) las matrices correspondientes a dos individuos estabilizantes, es decir, tales que A + BK1 < 0, A + L1C < 0

(6.16)

A + BK2 < 0, A + L2C < 0     Entonces, se cumple que los hijos K 1 , L1 y K 2 , L2 generados mediante el operador aritmético son ambos igualmente estabilizantes.     En efecto, sean K 1 , L1 y K 2, L2 las matrices que representan a los hijos generados a partir de (K1, L1 ) y (K2, L2 ), a saber K 1 = αK1 + (1 − α)K2

(6.17)

K 2 = (1 − α)K1 + αK2 L1 = αL1 + (1 − α)L2 L2 = (1 − α)L1 + αL2 donde α es un escalar perteneciente al intervalo [0,1]. Posteriormente se multiplican las ecuaciones (6.16) por α y (1 − α) de manera que se obtiene αA + αBK1 < 0, αA + αL1 C < 0

(6.18)

(1 − α)A + (1 − α)BK2 < 0, (1 − α)A + (1 − α)L2 C < 0 (1 − α)A + (1 − α)BK1 < 0, (1 − α)A + (1 − α)L1 C < 0 αA + αBK2 < 0, αA + αL2 C < 0 las cuales se cumplen si α = 0, α = 1. Note que en los casos α = 0 o α = 1 los hijos son idénticos a los padres. Ahora se adicionan las ecuaciones (6.18) para obtener αA + αBK1 + (1 − α)A + (1 − α)BK2 < 0

(6.19)

(1 − α)A + (1 − α)BK1 + αA + αBK2 < 0 αA + αL1 C + (1 − α)A + (1 − α)L2C < 0 (1 − α)A + (1 − α)L1 C + αA + αL2C < 0 lo cual puede reescribirse como A + BK 1 < 0, A + L1C < 0 A + BK 2 < 0, A + L2C < 0

(6.20)

120 de donde se se deduce que los hijos generados de esta forma son ambos estabilizantes. En cuanto al operador de mutación, se propone utilizar un esquema similar al presentado en el capítulo 4. La mutación se realiza fijando α = 1 y a partir de un individuo factible x se genera una mutación mediante x = x + αδ

(6.21)

donde δ es un vector aleatorio de distribución gaussiana. Si el individuo resultante no es factible entonces se disminuye α a la mitad y se genera una nueva mutación hasta que el resultado sea factible. La generación de la población inicial de MOPPEA se basa en el algoritmo recursivo 6.2. Sean pk ∈ Cn , pl ∈ Cn los autovalores de A + BK , A + LC respectivamente. Para asegurar la condición de estabilidad a lazo cerrado, las matrices K y L deben calcularse de tal forma que pk y pl pertenezcan al semi-plano C− , lo cual puede lograrse mediante el método presentado en el capítulo 2, sección 2.4. La idea principal es decidir al azar si se genera un polo real o un par de polos conjugados y a partir de allí formular la recursividad. El objetivo consiste en generar individuos factibles, distribuidos uniformemente en la zona del semi-plano complejo izquierdo delimitada por dos cotas denotadas rmax , imax que deben ser ajustadas en función de la posición de los polos del sistema a lazo abierto. Dado que en general se recomienda perturbar lo menos posible las coordenadas de los polos a lazo abierto, es posible calcular las cotas rmax , imax a partir de dichas coordenadas, dejando sin embargo un cierto márgen para que sea posible la reubicación en caso de que se estime necesario. Una vez generada la población inicial, la fase de exploración puede realizarse mediante cualquier algoritmo genético que permita la búsqueda multi-objetivo. En este trabajo se propone comparar los resultados obtenidos por SPEA2 y SPEA2 - ACPL con respecto al método basado en factibilidad LMI, tal como se presenta en la próxima sección.

6.5. Resultados numéricos Los modelos que se estudian en esta sección fueron extraidos de COMP le ib, librR (Leibfritz, 2004). Una ería que contiene más de 120 sistemas LTI en formato MATLAB justificación para esta selección, es que los modelos estudiados previamente por Herreros,

121

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Entradas : n : grado del polinomio caracerístico del sistema a lazo cerrado rmax : cota superior de los polos en el eje real imax : cota superior de los polos en el eje imaginario Salidas : p Si n = 1 entonces p = −rand ∗ rmax; sino Si n = 2 entonces Si rand < 0.5 entonces p1 = −rand ∗ rmax ; p2 = −rand ∗ rmax ; sino p1 = −rmax ∗ rand + i ∗ imax ∗ rand; p2 = −rmax ∗ rand − i ∗ imax ∗ rand; finsi sino Si rand < 0.5 entonces p = [−rmax ∗ rand, genpol(n − 1, rmax , imax)]; sino p1 = −rmax ∗ rand + i ∗ imax ∗ rand ; p2 = −rmax ∗ rand − i ∗ imax ∗ rand; p = [p1, p2, genpol(n − 2, rmax , imax)] finsi finsi finsi Algoritmo 6.2 : Generación de un vector de polos aleatorio

Molina-Cristobal y Takahashi corresponden a sistemas SISO de bajo orden, más adaptados a la evaluación de controladores PID, para lo cual de hecho fueron propuestos (Astrom et al., 1998). Dado que originalmente COMP leib fue propuesta para la evaluación de soluciones al problema H∞ , para adaptar los modelos al caso H2/H∞ que se estudia en este capítulo, se propone extraer únicamente las matrices A, B y C, conservando la estructura que aparece en la ecuación (6.1). Se propone en concreto evaluar el funcionamiento de tres algoritmos para resolver el problema PH2 /H∞ :

122 • A1: que corresponde al algoritmo SPEA2 utilizando el método MOPPEA para la generación y el mantenimiento de individuos estabilizantes. • A2: que corresponde al algoritmo memético SPEA2-ACPL-E1 (ver capítulo 5) utilizando el método MOPPEA para la generación y el mantenimiento de individuos estabilizantes. • A3: que corresponde al algoritmo basado en factibilidad de LMIs (algoritmo 6.1). En la tabla 6.1 se presentan los datos correspondientes a los modelos seleccionados para las pruebas. Cada modelo COMP le ib se caracteriza por una nomenclatura específica. Para este trabajo se consideran cinco en particular, a saber: AC1, AC6, NN10, WEC1 y AC9. Note que los mismos fueron seleccionados de forma que el número de variables de decisión oscile entre 30 y 90, con el fin de analizar la influencia de esta cantidad sobre los resultados obtenidos por cada algoritmo. También se muestran en la tabla 6.1 los valores de los parámetros NH2 /H∞ ,∆γ 2, ∆γ ∞ , γ min y γ max correspondientes a la configuración del algoritmo basado en LMIs. Los mismos fueron ajustados mediante ensayo y error tomando en cuenta las soluciones extremas que se obtienen considerando cada objetivo de manera aislada, utilizando los comandos h2syn R. y hinfsyn de MATLAB Los parámetros rmax , imax que también aparecen en la tabla 6.1 corresponden a las cotas utilizadas para la generación de la población inicial de MOPPEA, calculadas para incluir los polos de cada modelo a lazo abierto. Tabla 6.1: Datos correspondientes a COMP le ib n nu ny nMOPP EA NH2 /H∞ AC1 5 3 3 30 100 AC6 7 2 4 42 100 NN10 8 3 3 48 100 WEC1 10 3 4 70 100 AC9 10 4 5 90 100

los problemas COMP leib ∆γ 2 ∆γ ∞ γ min γ max rmax 0.01 0.01 1.5 5.5 10 0.1 1 1 110 20 0.01 0.02 3 15 10 0.05 1 2 200 100 0.1 0.1 50 100 100

imax 10 20 10 100 100

123 Es importante resaltar que adicionalmente a los algoritmos A1, A2 y A3 antes mencionados se realizaron pruebas utilizando la representación tradicional en forma polinomial y en espacio de estados. Sin embargo, no fue posible conseguir individuos estables utilizando estos formatos, por lo cual no se incluyen en los resultados que se presentan a continuación. Se realizaron 30 ejecuciones de cada algoritmo. En las figuras 6.3, 6.4, 6.5, 6.6 y 6.7 se puede apreciar un ejemplo de las aproximaciones obtenidas mediante cada uno.

2.5 Ejemplo de aproximación obtenida con A1 Ejemplo de aproximación obtenida con A2 Ejemplo de aproximación obtenida con A3

2.4 2.3 2.2

H

∞

2.1 2 1.9 1.8 1.7 1.6

0.8

1

1.2

1.4

1.6

1.8

2

H2

Figura 6.3: Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema AC1

124

160 Ejemplo de aproximación obtenida con A1 Ejemplo de aproximación obtenida con A2 Ejemplo de aproximación obtenida con A3

140 120

H∞

100 80 60 40 20

2

3

4

5 H

6

7

8

2

Figura 6.4: Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema AC6

20

18

H

∞

16

14

12 Ejemplo de aproximación obtenida por A1 Ejemplo de aproximación obtenida por A2 Ejemplo de aproximación obtenida por A3

10

8 2.5

3

3.5

4

4.5

5

H2

Figura 6.5: Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema NN10

125

Ejemplo de aproximación obtenida mediante A1 Ejemplo de aproximación obtenida mediante A2 Ejemplo de aproximación obtenida mediante A3

180

160

H∞

140

120

100

80

3

4

5 H

6

7

8

2

Figura 6.6: Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema WEC1

800 700 600

Ejemplo de aproximación obtenida por A1 Ejemplo de aproximación obtenida por A2 Ejemplo de aproximación obtenida por A3

H∞

500 400 300 200 100 0

0

50

100

150

200

250

300

350

400

450

H

2

Figura 6.7: Ejemplo de las soluciones obtenidas por A1, A2 y A3. Problema AC9

126 En estas figuras es posible apreciar que entodos los casos A2 produce individuos que en su mayoría dominan a los de A1. Para el problema AC1 (30 variables), se aprecia que A2 mejora los resultados de A3. En el caso de AC6, las soluciones de A2 y A3 son mutuamente no-dominadas. En los casos NN10, WEC1 y AC9 los resultados del algoritmo basado en LMIs mejoran claramente los resultados del genético, considerando un presupuesto fijo de 20.000 evaluaciones del vector objetivo. Esta cantidad de evaluaciones permite que las 30 ejecuciones del algoritmo genético no excedan el tiempo máximo permitido para completar el proceso de diseño (≈24 horas). La configuración utilizada para la ejecución de las pruebas se muestra en la tabla 6.2. Los resultados obtenidos para los indicadores C, DMAX y ESS se muestran en las tablas 6.3 y 6.4. Independientemente del problema, es posible observar que los resultados obtenidos por A2 son mejores que los de A1 con respecto al operador de cobertura, lo cual ilustra la utilidad del operador ACPL para la búsqueda local en estos casos. Por otra parte, el algoritmo A2 logra mejorar las soluciones de A3 en el caso del problema AC1 e igualar en el caso AC6, con respecto al operador de cobertura. En los otros dos casos, el algoritmo basado en LMIs es claramente superior, tal como fue posible constatar de manera gráfica. Tabla 6.2: Parámetros utilizados por SPEA2 y SPEA2-ACPL Parámetros Valor Representación K+L Operador de cruce Aritmético Probabilidad de cruce 0.9 Operador de mutación Gaussiana Probabilidad de mutación 0.1 Probabilidad de aplicación del ACPL 0.25 ε [0.1, 0.1] r [0.1, 0.1] Niter 10 Nnd [5, 5] T ol 10−4 Condición de parada 200 generaciones Número máximo de evaluaciones 20.000 Tamaño máximo del archivo externo 20 Esquema de selección (100 + 100)

127

C(Ai , Aj ) A1

A2

A3

Tabla 6.3: Resultados del indicador de cobertura A1 A2 A3 AC1:0.0200(0.0535) AC1:0.8619(0.1326) AC6:0.0217(0.0773) AC6:0(0) NN10:0.3650(0.3996) NN10:0(0) − WEC1:0.2500(0.2603) WEC1:0(0) AC9:0.2978(0.3548) AC9:0(0) AC1:0.8517(0.1850) AC1:0.9810(0.0494) AC6:0.6750(0.2589) AC6:0(0) NN10:0.4600(0.4438) NN10:0(0) − WEC1:0.3867(0.2846) WEC1:0(0) AC9:0.4850(0.4329) AC9:0(0) AC1:0(0) AC1:0(0) AC6:0(0) AC6:0(0) NN10:1(0) NN10:1(0) − WEC1:1(0) WEC1:1(0) AC9:1(0) AC9:1(0)

Tabla 6.4: Resultados de los indicadores ESS y DMAX DMAX ESS AC1:0.7589(0.0535) AC1:0.0190(0.0069) AC6:120.4263(67.0600) AC6:3.5698(4.0164) NN10:0.1026(0.0407) A1 NN10:3.5940(0.7046) WEC1:45.3320(7.7541) WEC1:1.0784(0.4881) AC9:564.2026(190.4456) AC9:21.5037(17.8440) AC1:0.9638(0.3524) AC1:0.0253(0.0121) AC6:73.0991(15.1735) AC6:1.7792(0.5986) NN10:0.0998(0.0298) A2 NN10:3.7460(0.6381) WEC1:44.4208(6.2985) WEC1:1.0812(0.3809) AC9:261.0278(169.0322) AC9:14.7817(19.8238) AC1:0.1618(0) AC1:0.0122(0) AC6:2.8065(0) AC6:0.0152(0) NN10:1.0200(0) NN10:0.0188(0) A3 WEC1:0.1540 WEC1:36.0976(0) AC9:59.2017(0) AC9:3.3800(0)

128 De igual forma, con respecto al indicador ESS, el algoritmo A3 logra mejorar los resultados obtenidos por los otros dos, lo cual se explica por el caracter determinista del algoritmo basado en LMIs, el cual intenta producir una distribución uniforme en el frente, caracterizada por los parámetros ∆γ 2 , ∆γ ∞. En relación con el índice DMAX no se aprecian claras diferencias entre A1 y A2, aunque ambos algoritmos producen frentes mucho más amplios que A3, para todos los problemas. Tampoco con respecto a ESS se aprecian notables diferencias entre A1 y A2.

6.6. Resumen del capítulo En este capítulo se propuso un nuevo método de diseño de controladores multiobjetivo llamado MOPPEA, que consiste fundamentalmente en: 1. Asignar a cada individuo una representación de la forma (K, L), donde K es la matriz de realimentación y L es la matriz de observación de estados. 2. Utilizar un algoritmo recursivo para generar la población inicial de individuos factibles. 3. Utilizar un operador de cruce aritmético que permite mantener a la población dentro de la región factible y un algoritmo genético multi-objetivo para llevar a cabo el proceso de búsqueda. En particular se presentó una comparación entre los resultados obtenidos por SPEA2 y SPEA2 - ACPL con respecto al método basado en factibilidad LMI, cuando se intenta resolver el problema mixto H2 /H∞ en el caso de cinco modelos de COMP le ib. Las pruebas estadísticas realizadas permiten afirmar que, considerando un presupuesto fijo de 20.000 evaluaciones, el enfoque genético es eficiente únicamente para los problemas AC1 y AC6 caracterizados por 30 y 42 variables de decisión respectivamente. Por el contrario, en el caso de los pro-blemas NN10, WEC1 y AC9 (con 48, 70 y 90 variables) las soluciones obtenidas mediante LMIs dominan a las producidas por SPEA2 y SPEA2 ACPL. Por otra parte, en todas las pruebas realizadas, los valores del indicador de cobertura obtenidos por el algoritmo memético SPEA2-ACPL superan a los obtenidos por el SPEA2, lo cual confirma la conveniencia de utilizar el operador de búsqueda local en esta clase de problemas.

CAPÍTULO VII CONCLUSIONES

La principal contribución de este trabajo consiste en la propuesta y el análisis de nuevos operadores, especialmente concebidos para intentar mejorar el proceso de diseño de controladores automáticos lineales mediante algoritmos genéticos. Sin embargo, algunas estrategias pueden ser extrapoladas para resolver problemas de ingeniería en otras áreas tales como identificación de sistemas y análisis de señales. Por ejemplo, el operador ACPL puede ser adaptado para intentar resolver otros problemas de optimización multi-objetivo, previa adecuación de sus parámetros y escalamiento de variables. En este sentido, se presentaron evidencias que permiten concluir que para problemas que no poseen muchos frentes locales, la utilización del ACPL conduce a mejoras estadísticamente significativas del indicador de cobertura, con respecto al algoritmo genético puro. Este resultado confirma las conclusiones obtenidas por otros autores en trabajos previos similares. También con respecto al operador ACPL, es importante destacar que la codificación binaria propuesta para la identificación de los conos de diversidad y el cálculo de la dirección correspondiente a cada cono son aspectos sobre los cuales no se han conseguido referencias previas en la revisión bibliográfica realizada. Por otra parte, aunque los mecanismos de integración entre el algoritmo genético y el operador de búsqueda local explorados en el transcurso de esta investigación fueron muy sencillos, los resultados obtenidos son prometedores y alientan la continuación en esta línea de investigación. El método propuesto para el diseño PID basado en el cálculo exacto de la región estabilizante produjo resultados satisfactorios en un caso de estudio basado en un modelo a lazo abierto de orden seis y de fase no-mínima, de una complejidad superior a los analizados por otros autores. Es posible que su aplicación se justifique en un caso práctico, puesto que permitiría mejorar considerablemente la calidad de las soluciones y disminuir el tiempo del proceso de diseño. Una extensión posible consiste en evaluar el método en aplicaciones en

130 línea (auto-tuning y control adaptivo), puesto que es posible asegurar que los controladores obtenidos son al menos estabilizantes y por ende no pondrían en riesgo la integridad de la planta controlada. El tema de la inestabilidad es un punto sobre el cual se hizo mucho énfasis en este trabajo. Es curioso que ninguno de los investigadores que adoptaron previamente la metodología genética para resolver problemas de control haya tomado en cuenta este aspecto tan fundamental con mayor interés, optando en la mayoría de los casos por un enfoque basado en una función de penalidad. En este sentido, el método MOPPEA permite generar y mantener una población compuesta exclusivamente por individuos estabilizantes: algo que, de acuerdo con la revisión de la literatura realizada, ningún otro algoritmo genético había garantizado plenamente. Y esto gracias a las ventajas de la representación matricial (K, L) para la exploración del espacio de controladores estabilizantes, con respecto a la representación polinomial o en espacio de estados utilizada por otros autores. Las pruebas realizadas confirmaron que en efecto estas últimas no son nada eficientes en el caso de sistemas MIMO de alto orden. La representación (K, L) propuesta en este trabajo es preferible ya que toma en cuenta de manera intrínseca la restricción de estabilidad. Es posible afirmar, a partir de los resultados obtenidos, que la representación genética de los controladores en forma polinomial o en forma controlable únicamente se justifica si el sistema a lazo abierto es SISO y de bajo orden. En nuestro conocimiento, ninguno de los investigadores dedicados previamente al problema mixto H2/H∞ ha publicado información relacionada con los límites de la eficiencia de los algoritmos genéticos con respecto a las soluciones LMI. En este trabajo se abordó este tema a partir de casos de estudio en los cuales el número de variables de decisión oscila entre 30 y 90, considerando un presupuesto fijo de 20.000 evaluaciones. Se pudo concluir que MOPPEA es capaz de mejorar las soluciones LMI en el caso del problema de 30 variables y de igualarlas en el caso de 42 variables. Para un mayor número mayor de variables las soluciones LMI son claramente superiores, tomando en cuenta el presupuesto limitado de evaluaciones del vector objetivo.

131 Estos resultados confirman que, como era de esperarse, mientras más grande es el número de variables de decisión, menos competitivo se hace el enfoque genético con respecto a las LMIs. Queda como punto pendiente analizar lo que sucede si se aumenta el presupuesto de evaluaciones, para lo cual se requiere la disponibilidad de equipos con mayor capacidad de cómputo, de manera que se pueda obtener resultados en un tiempo adecuado. Esto se debe a que es necesario repetir las pruebas tantas veces como sea posible, con el fin de obtener medidas estadísticamente significativas. Otro punto de investigación que queda abierto consiste en lograr una formulación del problema LMI en términos de las matrices K y L de la estructura “Observación + Realimentación de Estados”. Esta formulación facilitaría una interacción que podría conducir a mejores soluciones con respecto a cada método de diseño considerado de manera aislada. Una innovación interesante de MOPPEA consiste en el algoritmo recursivo para la generación de la población inicial. En su versión actual este algoritmo permite generar lazos cerrados cuyos polos pertenecen a una región rectangular delimitada por dos cotas. Sin embargo, en el futuro es posible adaptar el algoritmo para que los polos se distribuyan en otras regiones del plano complejo, tales como discos o conos. También es importante señalar que aunque los indicadores de desempeño utilizados a lo largo del trabajo (Cobertura, Distancia Máxima y Espaciado) son ampliamente conocidos por los investigadores en el área de algoritmos genéticos multi-objetivo, tampoco ha sido posible conseguir referencias previas relacionadas con el cálculo de indicadores de desempeño y pruebas estadísticas similares a las presentadas en este trabajo, en el caso de problemas de diseño de controladores. Es recomendable que los procedimientos de evaluación de las soluciones multi-objetivo se conviertan en un estándar para este tipo de aplicaciones, ya que permite la comparación entre los resultados obtenidos por diversos autores. A lo largo del trabajo se puso mucho énfasis en la generación de la población inicial. De hecho, en este sentido se propusieron tres novedosos métodos con esta finalidad: el método exacto que calcula los vértices de la región estabilizante en el caso PID, el método que calcula una aproximación mediante muestreo aleatorio y el método recursivo para la generación de individuos factibles (MOPPEA). Es interesante notar que estos métodos se inspiran de diversas fuentes tales como la teoría clásica de control, la programación

132 matemática, la geometría, la topología y los algoritmos genéticos. Cada disciplina fue abordada intentando obtener lo mejor de cada una para alcanzar el objetivo final, el cual no es otro que conseguir las mejores soluciones posibles para los problemas de diseño. De esta forma se ha confirmado la necesidad de incorporar conocimiento específico en las aplicaciones de los algoritmos genéticos, como consecuencia del famoso teorema “No free lunch” (Wolpert & Macready, 1997). Muchos puntos quedan sin embargo por abordar, entre los cuales es posible mencionar los siguientes: • Implementar algoritmos genéticos adaptivos, para disminuir la necesidad de ajustes de parámetros por parte del diseñador. • Continuar el estudio sobre la interacción de los operadores de búsqueda local con los operadores de búsqueda global (explotación vs exploración). • Considerar conceptos novedosos de dominancia, para permitir la inclusión de las preferencias del diseñador, antes o después del proceso de búsqueda. • Mejorar las interfaces con el usuario y del proceso de selección de soluciones (toma de decisiones). • Desarrollo de métodos multi-objetivo de optimización estructural y mixto-entera. • Continuar el estudio sobre representaciones del espacio de controladores, principalmente en el caso no-lineal de múltiples entradas y salidas. • Consideración de problemas de modelaje e identificación de la planta en el contexto del diseño multi-objetivo. • Resolver problemas de control que necesiten optimización multi-objetivo en línea, considerando las dificultades tecnológicas que ello implica.

133

APÉNDICE Se denomina señal a cualquier función f : R → Cn que represente una magnitud física o lógica variable en el tiempo. A continuación se presentan algunos resultados teóricos relacionados con señales, operadores y normas en diversos espacios de interés para el análisis de sistemas dinámicos. Se denota Ln2 (−∞, ∞) al espacio lineal compuesto por las señales f : R → Cn continuas y deterministas, tales que:
+∞ f ∗(t)f (t)dt < ∞ (7.1)

−∞

La función ·2 : Ln2 (−∞, ∞) → R+ 0 tal que

 
+∞   f 2 =  f ∗(t)f (t)dt, −∞

∀f ∈ Ln2

(7.2)

define una norma sobre Ln2 (−∞, ∞). Esta cantidad puede ser interpretada como la energía total contenida en f . Se denota C+ al semi-plano complejo abierto derecho y C cerrado derecho, esto es

+

al semi-plano complejo

C+ ≡ {c ∈ C; Re(c) > 0} C

+

≡ {c ∈ C; Re(c) ≥ 0}

ˆ n2 (jR) al espacio de Hilbert compuesto por las funciones matriciales Se denota L fˆ : jR → Cn

(7.3)

dotado del producto interno +

,
+∞ ˆ ˆ fˆ1∗(jω)fˆ2 (jω)dω f1 , f2 = −∞

(7.4)

p×m

El operador lineal acotado Γ : L2

ˆ n2 (jC) definido por (−∞, ∞) → L


+∞ f (t)e−jωtdt Γ(f ) = fˆ(jω) = −∞

134

(7.5)

se denomina Transformada de Fourier. Γ es un operador unitario e invertible, esto es, Γ ˆ n2 (jC). define un isomorfismo entre Ln2 (−∞, ∞) y L Se denota H2p×m al espacio de funciones matriciales complejas

+ fˆ : C → Cp×m tales que:

+ • fˆ es analítica en C . • Para cualquier ω ∈ R , excepto en un conjunto finito {ω 1 , ω 2, ...ω N } , se cumple: fˆ(jω) = lim+ fˆ(σ + jω) σ→0

(7.6)

• La siguiente cantidad es finita
+∞ sup traza[fˆ∗ (σ + jω)fˆ(σ + jω)]dω σ≥0 −∞ La función ·2 : H2p×m → R+ 0 tal que  
+∞    1   ˆ traza[fˆ∗ (σ + jω)fˆ(σ + jω)]dω f  = sup 2 σ≥0 2π −∞

(7.7)

(7.8)

∀f ∈ H2p×m , define una norma sobre H2p×m , denotada H2 .

+ Una función matricial fˆ : C → Cp×m se denomina real racional, si cada uno de sus elementos fˆij (s) es de la forma: bnu snu + bnu−1snu−1 + .. + b1 s + b0 Nij (s) fˆij (s) = nd = s + and−1snu−1 + .. + a1s + a0 Dij (s) donde todos los coeficientes de los polinomios Nij (s) y Dij (s) son reales.

(7.9)

135 El conjunto de funciones matriciales fˆ : C → Cp×m reales racionales se denota R(s)p×m . Si ∀fˆij (s), nd ≥ nu , fˆ es llamada propia. Si ∀fˆij (s), nd > nu , fˆ es lla+ mada estrictamente propia. Para cualquier fˆij (s), se denomina “polo” a cualquier p ∈ C     tal que fˆij (s) → ∞, cuando s → p. Una función real racional fˆ pertenece a H2p×m si y + solo si es estrictamente propia y no posee polos en C . Se denota RH2p×m al sub-espacio de funciones matriciales reales racionales pertenecientes a H p×m .

+

2

El operador lineal acotado Λ : Ln2 [0, ∞) → H2n definido por:
+∞ f (t)e−stdt Λ(f ) = fˆ(s) =

(7.10)

0

se denomina Transformada Unilateral de Laplace. Λ define un isomorfismo entre Ln2 [0, ∞) y H2n .

Se denota Ln∞ (−∞, ∞) al espacio lineal de señales continuas y determinísticas, comn puesto por las funciones f : R+ 0 → C tales que: ess sup t∈(−∞,∞)



f ∗(t)f (t) < ∞

(7.11)

La función ·∞ : Ln∞ (−∞, ∞) → R+ 0 tal que f ∞ = ess sup t∈(−∞,∞)



f ∗ (t)f (t)

(7.12)

∀f ∈ Ln∞ define una norma sobre Ln∞(−∞, ∞) que puede ser interpretada como el módulo máximo de la señal f para t ∈ (−∞, ∞). p×m al espacio de funciones matriciales G ˆ : C+ → Cp×m tales que: Se denota H∞

ˆ es analítica en C+ . • G • Para cualquier ω ∈ R , excepto en un conjunto finito {ω 1 , ω 2, ...ω N }, se cumple: ˆ + jω) ˆ G(jω) = lim+ G(σ σ→0

(7.13)

• La siguiente cantidad es finita: ˆ ess sup σ[G(s)] + s∈C

(7.14)

136 ˆ donde σ representa el máximo valor singular de la matriz G(s), es decir, ' ˆ ˆ ∗G(s)] ˆ σ[G(s)] =max λi [G(s) i

con,

(7.15)

ˆ ∗G(s)] ˆ ˆ ∗G(s) ˆ ≡ i-ésimo autovalor de G(s) λi [G(s)

p×m → R+ tal que La función ·∞ : H∞ 0    ˆ ˆ G =sup σ[G(jω)], ∞ ω∈R

p×m ˆ ∈ H∞ ∀G

(7.16)

p×m denotada H . De hecho suponga u ∈ H nu . Entonces, G ˆ define una norma sobre H∞ ∞ 2 ny ×nu n y ˆ ∈ H . El se comporta como un operador lineal multiplicativo, esto es, y = Gu ∈ H∞ 2 ˆ operador G es comunmente llamado “matriz de transferencia” desde u hacia y. Además, ˆ ∈ H p×m define un operador mediante el operador Transformada de Laplace, la matriz G ∞ p m lineal acotado y causal G de L2 [0, ∞) → L2 [0, ∞), de la forma:

ˆ G = Λ−1GΛ

(7.17)

p×m si y solo si G ˆ es propia y no posee ˆ ∈ R(s)p×m , pertenece a H∞ Una función racional G + p×m se polos en C . El conjunto de funciones matriciales racionales pertenecientes a H∞ p×m . denota RH∞

REFERENCIAS Ackermann, J. & Kaesbauer, D. (2003). Stable polyhedra in parameter space. Automatica, (39):937—943. Adra, S. ; Amody, A. ; Griffin, I. & Fleming, P. (2005a). A hybrid multi-objective evolutionary algorithm using an inverse neural network for aircraft control system design. In Proceedings of the IEEE Congress on Evolutionary Computation, pp 1—8. Adra, S. ; Griffin, I. & Fleming, P. (2005b). Hybrid multi-objective genetic algorithm with a new adaptive local search process. In Hans-Georg-Beyer et al. editor, Genetic and Evolutionary Computation Conference GECCO’2005, volume 1, pp 1009—1010. Andres, B. ; Jiron, J. ; Fernandez, P. ; Lopez, J. & Besada, E. (2004). Multiobjective optimization and multivariable control of the beer fermentation process with the use of evolutionary algorithms. Journal of Zhejiang University - Science, pp 378—389. Apkarian, P. ; Noll, D. & Rondepierre, A. (2008). Mixed H2 /H∞ control via nonsmooth optimization. SIAM Journal on Control and Optimization, 47(3):1516—1546. Arakawa, M. ; Nakayama, H. ; Hagiwara, I. & Yamakawa, H. (1998). Multi-objective optimization using adaptive range genetic algorithms with data envelopment analysis. In A Collection of Technical Papers of the 7th Symposium on Multidisciplinary Analysis and Optimization, pp 2074—2082. Astrom, K. & Hagglund, T. (1995). PID controllers: theory, design and tuning. Instruments Society of America. Astrom, K. ; Panagopoulos, H. & Hagglund, T. (1998). Design of PI controllers based on non-convex optimization. Automatica, 34(5):585—601. Bertalanffy, L. (1950). An outline of general systems theory. British Journal for the Philosophy of Science, 1(2). Boyd, S. ; El-Ghaoui, L. ; Feron, E. & Balakrishnan, V. (1994). Linear matrix inequalities in system and control theory. SIAM Studies in Applied Mathematics. Brown, M. & Smith, R. (2005). Directed multi-objective optimization. International Journal on Computers, Systems and Signals, 6(1):3—17. Broyden, C. (1970). The Convergence of a Class of Double-rank Minimization Algorithms. Journal of the Institute of Mathematics and Its Applications, 6:76—90. Calafiore, G. & Dabbene, F. (2006). Probabilistic and randomized methods for design under uncertainty. Springer-Verlag, London.

138 Caponio, A. ; Cascella, G. ; Neri, F. ; Salvatore, N. & Sumner, M. (2007). A fast adaptive memetic algorithm for online and offline control design of PMSM drives. IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, 37(1):28—41. Cesareo, J. (1997). Estrategias evolutivas y su aplicación en la síntesis de controladores. PhD thesis, Universidad de Vigo. Dpto. de Ingeniería de Sistemas. Chipperfield, A. ; Bica, B. & Fleming, P. (1998). Fuzzy scheduling control of a gas turbine aero-engine: a multi-objective approach. IEEE Transactions on Industrial Electronics, 49(3):536—548. Chipperfield, A. & Fleming, P. (1996). Multi-objective gas turbine engine controller desing using genetic algorithms. IEEE Transactions on Industrial Electronics, 43(5):583—587. Clement, B. (2001). Synthese multiobjectifs et sequencement de gains: Application au pilotage d’un lanceur spatial. PhD thesis, Service Automatique de SUPELEC - Gif sur Yvette, France. Coello, C. ; Lamont, G. & Van-Veldhuizen, D. (2007). Evolutionary algorithms for solving multi-objective problems. Springer. Genetic and Evolutionary Computation Series. Second Edition. Dawkins, R. (1976). The selfish gene. Oxford University Press. Deb, K. (2001). Multiobjective optimization using evolutionary algorithms. Jhon Wiley and Sons Ltd. Deb, K. & Agrawal, R. (1995). Simulated binary crossover for continuous search space. Complex Systems, 9(2):115—148. Deb, K. ; Agrawal, S. ; Pratab, A. & Meyarivan, T. (2000). A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II. In Proceedings of the Parallel Problem Solving from Nature VI Conference, pp 849—858. Dellino, G. ; Lino, P. ; Meloni, C. & Rizzo, A. (2007). Enhanced evolutionary algorithms for multidisciplinary design optimization: a control engineering perspective. In Hybrid Evolutionary Algorithms - Studies in Computational Intelligence, pp 173— 216. Springer Berlin / Heidelberg. Dennis, J. & Schnabel, R. (1983). Numerical methods for unconstrained optimization and nonlinear equations. Prentice-Hall. New York. Doyle, J. ; Glover, K. ; Khargonekar, P. & Francis, B. (1989). State-space solutions to standard H2 and H∞ control problems. IEEE Transactions on Automatic Control, 34(8):831—847. Duarte, N. ; Ruano, A. ; Fonseca, C. & Fleming, P. J. (2000). Accelerating multiobjective control system design using a neuro-genetic approach. Proceedings of the 8th IFAC Symposium on Computer Aided Control Systems Design. University of Salford, U.K. Dupuis, A. ; Ghribi, M. & Kaddouri, A. (2004). Multi-loop PI tuning using a multiobjective genetic algorithm. In Proceedings of the IEEE International Conference on

139 Industrial Technology, pp 1505—1510. Eiben, A. & Smith, J. (2003). Introduction to evolutionary computing. Springer-Verlag. Fleming, P. (2001). Genetic Algorithms in Control System Engineering. Technical report, No. 789. University of Sheffield. UK. Fleming, P. J. & Purshouse, R. (2002). Evolutionary algorithms in control systems engineering: a survey. Control Engineering Practice (10), pp 1223—1241. Fonseca, C. (1995). Multi-objective genetic algorithms with application to control engineering problems. PhD thesis, Departement of Automatic Control and Systems Engineering. University of Sheffield. Fonseca, C. & Fleming, P. (1994). Multiobjective optimal controller design with genetic algorithms. In Proc. IEE Control International Conference, volume 1, Warwick, U.K, pp 745—749. Gambier, A. ; Wellenreuthera, A. & Badreddin, E. (2006). A new approach to design multi-loop control systems with multiple controllers. In Proceedings of the 45th IEEE Conference on Decision and Control. San Diego, CA, USA, pp 416—423. Goel, T. & Deb, K. (2002). Hybrid methods for multi-objective evolutionary algorithms. In Proceedings of the 4th Asia-Pacific conference on Simulated Evolution and Learning, pp 188—192. Granado, E. ; Colmenares, W. & Perez, O. (2007). Diseño de PID multivariable usando LMI iterativas. Revista de la Fac. Ing. UCV, 22(3):5—12. Griffin, I. ; Schroder, P. ; Chipperfield, A. & Fleming, P. (2000). Multi-objective Optimization Approach to the Alstom Gasifier Problem. In Proceedings of the I MECH E Part I Journal of Systems and Control Engineering, pp 453—468. Hart, W. (1994). Adaptive Global Optimization with Local Search. PhD thesis, University of California, San Diego. USA. Hassibi, A. ; Howy, J. & Boyd, S. (1999). Low-Authority Controller Design via Convex Optimization. AIAA Journal of Guidance, Control and Dynamics, 22(6):862—872. Heo, J. ; Lee, K. & Garduno-Ramirez, R. (2006). Multiobjective control of power plants using particle swarm optimization techniques. IEEE Transactions on Energy Conversion, 21(2):552—561. Herreros, A. (2000). Diseño de controladores robustos multiobjetivo por medio de algoritmos genéticos. Tesis Doctoral. Departamento de Ingenieria de Sistemas y Automatica. Universidad de Valladolid, España. Hillermeier, C. (2001). Nonlinear multiobjective optimization - A generalized homotopy approach. Birkhauser. Hindi, H. A. ; Hassibi, B. & Boyd, S. P. (1998). Multi-objective H2/H∞ optimal control via finite dimensional Q-parametrization and linear matrix inequalities. In Proceedings of the American Control Conference, pp 3244—3249.

140 Hohenbichler, N. & Abel, D. (2008). Calculating all KP admitting stability of a PID control loop. In Proceedings of the 17th World IFAC Congress. Seoul, Korea, July 6-11, pp 5802—5807. Hwang, C. & Hsiao, C.-Y. (2002). Solution of a non-convexoptimization arising in PIPID control design. Automatica 38, pp 1895—1904. Ishibuchi, H. & Murata, T. (1996). Multi-objective genetic local search algorithm. In Proc. of the IEEE Int. Conf. on Evol. Computation, Nagoya, Japan, pp 119—124. Joos, H. (2002). A multi-objective optimisation-based software environment for control systems design. In Proceedings of the IEEE International Symposium on Computer Aided Control System Design, pp 7—14. Kalman, R. (1960). Contributions to the Theory of Optimal Control. Boletín de la Sociedad Matemática Mexicana, (5):102—119. Kalman, R. (1963). Mathematical Description of Linear Dynamical Systems. SIAM Journal on Control, 1(2):152—192. Khor, E. ; Tan, K. & Lee, T. (2002). Learning the search range for evolutionary optimization in dynamic environments. Knowledge and Information Systems, 4(2):228—255. Knowles, J. & Corne, D. (2000). M-PAES: a memetic algorithm for multiobjective optimization. In Proceedings of the IEEE Congress on Evolutionary Computation, pp 325—332. Kumar, S. & Subramanian, S. (2002). Mutation rates in mammalian genomes. In Proceedings of the National Academy of Sciences of the United States of America, pp 803—808. Lagunas, J. (2004). Sintonización de controladores PID mediante un algoritmo genético multi-objetivo (NSGA-II). PhD thesis, Departamento de Control Automático. Centro de Investigación y de Estudios Avanzados. México. Lara, A. ; Sánchez, G. ; Coello-Coello, C. & Schüetze, O. (2010). HCS: A New Local Search Strategy for Memetic Multiobjective Evolutionary Algorithms. IEEE Transactions on Evolutionary Computation, 14(1):112—132. Leibfritz, F. (2004). COMPlib - Constrained Matrix-optimization Problem Library. Technical report, University of Trier. Department of Mathematics, Germany. Liu, G. ; Yang, J. & Whidborne, J. (2002). Multiobjective optimisation and control. Research Studies Press, Baldock, UK. Liu, T.-K. ; Ishihara, T. & Inooka, H. (1995). Multiobjective Control Systems Design by Genetic Algorithms. In Proceedings of the 34th SICE annual conference. Hokkaido, Japan, pp 1521—1526. Lyapunov, A. (1892). The general problem of stability motion. PhD thesis, University of Kharkov, Russia. Malpica, J. M. ; Fraile, A. ; Moreno, I. ; Obies, C. I. ; Drakec, J. W. & García-Arenal, F. (2002). The Rate and Character of Spontaneous Mutation in an RNA Virus.

141 Genetics, 162:1505—1511. Martinez, M. ; Sanchis, J. & Blasco, X. (2006). Algoritmos genéticos aplicados al diseño de controladores robustos. Revista Iberoamericana de Automática e Informática Industrial, 3(1):39—51. Miyamoto, S. (1997). Robust control design - a coprime factorization and LMI approach. PhD thesis, University of Cambridge. Molina-Cristobal, A. (2005). Multi-objective control: linear matrix inequality techniques and genetic algorithms approach. PhD thesis, University of Sheffield, UK. Moscato, P. (1989). On evolution, search, optimization, genetic algorithms and martial arts: toward memetic algorithms. Technical report, Caltech Concurrent Computation Program. Murata, T. ; S.Kaige & Ishibuchi, H. (2003). Generalization of Dominance Relation Based Replacement Rules for Memetic EMO Algorithms. In Proceedings of the Genetic and Evolutionay Computation Conference, pp 1234—1245. O’Dwyer, A. (2006). PI and PID controller tuning rules: an overview and personal perspective. In Proceedings of the IET Irish Signals and Systems Conference, pp 161—166. Osyczka, A. (1985). Multicriteria optimization for engineering design. Academic Press, Cambridge, MA. Pareto, V. (1906). Manuale di economia politica. Societa Editrice Libraria, Milan, Italy. Pedersen, G. (2005). Towards Automatic Controller Design using Multi-Objective Evolutionary Algorithms. PhD thesis, Department of Control Engineering. Aalborg University, Denmark. Popov, A. ; Farag, A. & Werner, H. (2005). Less Conservative Mixed H2/Hinf ControllerDesign Using Multi-Objective Optimization. Technical report, Technische Universitat Hamburg-Harburg. Reyes, O. ; Sánchez, G. & Strefezza, M. (2009). Multiobjective GA fuzzy logic controller. In Proceedings of the ICINCO 2009 Conference. Milan, Italy, pp 384—389. Rotea, M. & Khargonekar, P. (1991). H2-optimal Control with an Hinf-constraint: The State Feedback Case. Automatica, 27(2):307—316. Roy, B. (1971). Problems and methods with multiple objective functions. Mathematical Programming, 1(2):239—266. Sánchez, G. & Ferrer, J. (2004). A new approach for approximating the free transfer function in Q-parametrization. In Proceedings of the IEEE Conference Computer Aided Control Systems Design.Taipei,Taiwan, pp 339—343. Sánchez, G. ; Reyes, O. & Gonzalez, G. (2008a). Diseño de controladores PID multiobjetivo mediante algoritmos genéticos incluyendo el cálculo de los límites de la región de controladores estabilizantes. In XIII Congreso Latinoamericano de Control Automático. 25-28 de Noviembre. Mérida. Venezuela.

142 Sánchez, G. ; Reyes, O. & Strefezza, M. (2009). A Muti-objective Approach to Approximate the Stabilizing Region for Linear Control Systems. In Proceedings of the ICINCO 2009 Conference. Milan, Italy, pp 153—158. Sánchez, G. ; Villasana, M. & Strefezza, M. (2007). Multi-Objective Pole Placement with Evolutionary Algorithms. Springer-Verlag. Berlin. Series: Lecture Notes in Computer Science, 4403:417—427. Sánchez, G. ; Villasana, M. & Strefezza, M. (2008b). Solving Multi-Objective Linear Control Design Problems Using Genetic Algorithms. In Proceedings of the 17th IFAC World Conference, pp 12324—12329. Sasaki, D. ; Obayashi, S. & Nakahashi, K. (2002). Navier-Stokes optimization of supersonics wings with four objectives using evolutionary algorithms. Journal of Aircraft, 39(4):621—629. Scherer, C. ; Gahinet, P. & Chilali, M. (1997). Multi-objective output-feedback control via LMI optimization. IEEE Transactions on Automatic Control, 42(7):896—910. Schäffler, S. ; Schultz, R. & Weinzierl, K. (2002). A stochastic method for the solution of unconstrained vector optimization problems. J. Optimization Theory Applicat., 114(1):209—222. Schott, J. (1995). Fault tolerant design using single and multicriteria genetic algorithm optimization. Master’s thesis, Department of Aeronautics and Astronautics, MIT. Cambridge, MA. USA. Schütze, O. ; Sánchez, G. & Coello, C. A. C. (2008). A New Memetic Strategy for the Numerical Treatment of Multi-Objective Optimization Problems. Genetic and Evolutionary Computation Conference. July 12-16 ,Atlanta, Georgia, USA. Sharma, D. ; Kumar, A. ; Deb, K. & Sindhya, K. (2007). Hybridization of SBX based NSGA-II and sequential quadratic programming for solving multi-objective optimization problems. In Proceedings of the IEEE Congress on Evolutionary Computation, Singapore, pp 3003—3010. Shyr, W.-J. ; Wang, B.-W. ; Yeh, Y.-Y. & Su, T.-J. (2002). Design of optimal PID controllers using memetic algorithm. In Proceedings of the American Control Conference, pp 2130—2131. Silva, V. ; Fleming, P. ; Sugimoto, J. & R.Yokoyama (2008). Multiobjective optimization using variable complexity modelling for control system design. Applied Soft Computing 8, pp 392—401. Takahashi, R. ; Palhares, R. ; Dutra, D. & Goncalves, L. (2004). Estimation of Pareto sets in the mixed H2/H∞ control problem. International Journal of System Science, 35(1):55—67. Tan, K. ; Khor, E. & Lee, T. (2005). Multiobjective evolutionary algorithms and applications. Springer-Verlag. Villasana, M. & Ochoa, G. (2003). Heuristic Design of Cancer Chemotherapies. IEEE Transactions on Evolutionary Computation, Vol. 8, No. 6.

143 Wang, J. ; Brackett, B. T. & Harley, R. G. (2009). Particle Swarm-Assisted State Feedback Control: From Pole Selection to State Estimation. In Proceedings of the American Control Conference.St. Louis, MO, USA. June 10-12, pp 1493—1498. Wang, J. ; Wang, S. ; Wang, Z.-P. ; Tao, L. & Chen, J.-J. (2008). Robust controller design for a main steam pressure based on SPEA2. Journal of Zhengzhou University. Wanner, E. ; Guimaraes, F. ; Takahashi, H. & Fleming, P. (2006). A Quadratic Approximation-Based Local Search Procedure for Multiobjective Genetic Algorithms. In IEEE Congress on Evolutionary Computation. Vancouver, BC, Canada, pp 3361—3368. Wanner, E. ; Guimaraes, F. ; Takahashi, R. & Fleming, P. (2008). Local search with quadratic approximations into memetic algorithms for optimization with multiple criteria. Evolutionary Computation, 16(2):185—224. Wolpert, D. & Macready, W. (1997). No free lunch theorems for optimization. IEEE Transactions on Evolutinoary Computation, 1(1):67—82. Zakian, V. & Al-Naib, V. (1973). Design of dynamical and control Systems by the method of inequalities. Proceedings of the Institution of Electrical Engineers, pp 1421—1427. Zavala, J. ; Stewart, P. & Fleming, P. (2002). Multiobjective automotive drive by wire controller design. In Proceedings of the International Conference on Computer aided Control System Design, Glasgow, Scotland, pp 69—73. Ziegler, J. & Nichols, N. (1942). Optimum settings for automatic controllers. Transactions of ASME, pp 759—768. Zitzler, E. (1999). Evolutionary algorithms for multiobjective optimization: methods and applications. PhD thesis, Swiss Federal Institute of Technology. Zurich, Switzerland. Zitzler, E. ; Deb, K. & Thiele, L. (1999). Comparison of Multiobjective Evolutionary Algorithms on Test Functions of Different Difficulty. In Wu, A. S., editor, Proceedings of the 1999 Genetic and Evolutionary Computation Conference, pp 121—122, Orlando, Florida. Zitzler, E. ; Laumanns, M. & Thiele, L. (2001). SPEA2: improving the Strength Pareto Evolutionary Algorithm. Technical report, Computer Engineering and Networks Laboratory - TIK,ETH Zurich.