S.E.P.

S.E.S.

D.G.E.S.T.

CENTRO NACIONAL DE INVESTIGACIÓN Y DESARROLLO TECNOLÓGICO

cenidet

“EVALUACIÓN COMPARATIVA DE OTRAS ALTERNATIVAS AL ALGORITMO CLÁSICO DE MÍNIMOS CUADRADOS PARA LA IDENTIFICACIÓN DE SISTEMAS”

T E S I S QUE PARA OBTENER EL GRADO DE: MAESTRO EN CIENCIAS EN INGENIERÍA ELECTRÓNICA P R E S E N T A: ING. EDSON LÓPEZ MARTÍNEZ

DIRECTOR DE TESIS DR. VÍCTOR MANUEL ALVARADO MARTÍNEZ CUERNAVACA, MOR.

AGOSTO 2006

Dedicatoria Al creador, creador motor inmóvil y fuente generadora de energía que ilumina mi sendero. Por darme siempre una nueva oportunidad. A mis padres, padres Leodegario y Aida Delia: Por su inquebrantable amor y continuo aopoyo en la cimentación de mi estructura moral, espiritual y profesional; sin ustedes nada de esto sería posible, los amo. A mis Hermanas, Hermanas Alondra, Coral y Aideé: Ejemplos invaluables que me motivan e inspiran a ser un mejor ser humano. Siempre están ahí cuando las necesito, las amo. A mis Abuelos, Abuelos Román, Alicia, María y Néstor: Viejos, gracias por dejar en mí su legado difícil de superar y enseñarme a que lo importante no es vivir la vida su no saber cómo vivirla. A mis Tíos: Tíos: Por cada uno de sus atinados consejos en las diferentes etapas de mi vida. A mis Primos: Primos Por sus ánimos para alcanzar mis metas. A mis amigos, Carlos, Orlando, Hernando, Gustavo, Roberto II, Max, Ernesto, Jose Luis, Elyzabet y Mayith: Por su lealtad incondicional. A mi hija, hija Marysol López Alcántara: Alfa de mi proyecto de vida… A mi alma mater, Universidad Americana de Acapulco.

…Porque si no tengo amor, de nada me sirve…aunque tuviese profecía, y entendiese todos los misterios y toda ciencia…” 1 Corintios 13; 1:13.

Agradecimientos A Dios por darme la fortaleza para superar los obstáculos que se me presentaron durante esta etapa formativa de mi vida, ahora sé que es el único que nunca me fallará.

Al Dr. Víctor Manuel Alvarado Martínez, por su apoyo, confianza y paciencia brindada durante el desarrollo de este trabajo.

A los miembros del comité revisor de tesis: Dr. Marco Antonio Oliver Salazar, Dr. Alejandro Rodríguez Palacios y al Dr. Carlos Astorga Zaragoza por sus valiosos comentarios que enriquecieron este trabajo.

A mis profesores: Dra. Patricia Caratozzolo Martelliti, Dr. Hugo Calleja Gjumlich, Dr. Gerardo Vicente Guerrero Ramírez, Dr. Enrique Quintero-Mármol Márquez, Dr. Luis Gerardo Vela Valdés y al Dr. Jaime Arau Roffiel.

A mis compañeros de generación: Rafael Máxim Méndez Ocaña, José Luis Rullán, Ernesto Elías Vidal Rosas, Luis Adrian Sorcia Navarro, Israel Uribe Hernández, Gerardo Vázquez Guzmán, Abraham Cortés Dorantes y Javier Alejandro Molina Coronel, por su amistad y apoyo. Amigos la meta se alcanzó.

A mis demás compañeros con los cuales pasé buenos momentos: Chaca, chupis, chino, backs, cesarín, gracia, chikion, jojojó, rouse, ansioso, matis, paloma, pacheco, pita y el peque. Gracias por su amistad, apoyo y sobre todo gracias por su confianza.

Sin duda, he dejado de mencionar a muchas otras personas que hicieron más placentera mi estancia en Cuernavaca. No obstante, a todos ustedes les reitero mi más sincero agradecimiento.

De igual manera le agradezco al Consejo Nacional de Ciencia y Tecnología (CONACYT) y a la Secretaría de Educación Pública (SEP) por el apoyo económico que me brindó para la realización de este trabajo.

Finalmente, le agradezco al Centro Nacional de Investigación y Desarrollo Tecnológico CENIDET) por la oportunidad que me brindó para enriquecer mi formación profesional y personal.

Evaluación Comparativa de Otras Alternativas al Algoritmo Clásico de Mínimos Cuadrados para La Identificación de Sistemas

Edson López Martínez

Resumen La demanda de calidad en los procesos industriales ha motivado al surgimiento de nuevas técnicas de estimación de parámetros. El algoritmo tradicional para estimación de parámetros ha sido durante décadas el algoritmo de mínimos cuadrados. Sin embargo, el problema principal de este algoritmo radica en que presenta desempeño numérico muy pobre cuando la matriz de coeficientes es sobredimensionada, este trabajo de tesis aborda cuatros algoritmos alternativos: factorización LU, factorización QR, factorización UD Aumentada (AUDI) y la norma L1. Estos algoritmos alternativos se evalúan numéricamente con la finalidad de verificar que las factorizaciones: LU y QR se obtienen correctamente. Para evaluar los algoritmos: AUDI y la norma L1 se plantea un modelo sencillo tipo ARX para estimar sus parámetros. Las estructuras lineales que se emplean para la estimación paramétrica son ARX y ARMAX. Los cuatro casos de estudio para la validación de los algoritmos son: Sistema de tercer orden, el banco de nivel RCN100 que se encuentra en CENIDET, secadora de pelo y un sistema de calentamiento. El primer proceso es un sistema del cual se obtienen los datos de entrada/salida por medio de simulación, Por otro lado los datos de estos últimos tres procesos se obtuvieron de forma experimental. Las pruebas realizadas en la etapa de evaluación y validación muestran que la mayoría de los algoritmos son capaces de estimar los parámetros y caracterizar lo suficientemente bien a la dinámica del proceso real.

COMPARATIVE EVALUATION OF OTHER ALTERNATIVES TO THE CLASSICAL ALGORITHM OF LEAST-SQUARES FOR SYSTEM IDENTIFICATION

Edson López Martínez

Abstract

Demanding of quality in the industrial process has motivated to the rising of new parameter estimation techniques. The classical algorithm to parameter estimation during several decades has been the Least-Squares. However, the principal problem of this algorithm, is its poor numerical performance, when the coefficient matrix is overdeterminated; this thesis deal four alternative algorithms: LU decomposition, QR decomposition, AUDI identification, and L1 norm. These algorithms are numerically evaluated with a single goal, to verify that the factorizations are well done. To evaluate the algorithm based on L1 norm a single ARX model is established to estimate its parameters. The linear structures that are used to the parametric estimation are ARX and ARMAX. The process selected to validate the algorithms are four: 3th order system, the RCN100 apparatus level control that is in the automatic control laboratory of CENIDET, a dryer, and heating system. The first one is a system in which the input/output data are obtained by simulation, the other three handle data are obtained experimentally. Evaluation and validation results demonstrate that the algorithms are capable of estimating the parameters and characterize sufficiently the dynamics of real processes.

Contenido Lista de figuras

v

Lista de tablas

vii

Nomenclatura

ix

Introducción

1

Capítulo 1

5

Preliminares

5

1.1 Introducción 5 1.2 Antecedentes y problemática asociada 5 1.3 Algoritmos para la estimación de parámetros 6 1.3.1 Mínimos cuadrados 6 1.3.2 Algoritmos alternativos para la estimación de parámetros. 8 1.3.2.1 Factorización LU 8 1.3.2.2 Factorización QR 9 1.3.2.3 Factorización AUDI o UD aumentada 10 1.3.2.4 Norma L1 11 1.4 Hipótesis 12 1.5 Objetivos 12 1.5.1 Objetivo particular 12 1.5.2 Objetivos particulares 12 1.6 Metodología 13 1.7 Alcances y limitaciones del trabajo 13 1.8 Señal de excitación persistente 14 1.9 Proceso de identificación 14 1.9.1 Etapas a seguir para la identificación de un modelo 14 1.9.2 Selección y validación de la estructura del modelo. 17 1.9.2.1 Selección de la estructura del modelo 18 1.9.2.2 Análisis preeliminares de los datos para la selección de la estructura del modelo 18 1.9.2.3 Criterios para la comparación de la estructura de los modelos 19 1.9.2.4 Validación del modelo 21 1.10 Conclusiones 22 Capítulo 2

23

Factorización LU

23

2.1 Introducción 2.2 Definiciones y conceptos Generales 2.3 Eliminación de Gauss con intercambio (con pivote) i

23 23 26

Contenido

2.3.1 Eliminación de Gauss-Jordan 2.3.2 Inversa mediante la eliminación de Gauss 2.4 Existencia de soluciones a sistemas de ecuaciones 2.5 LU y eliminación de Gauss 2.5.1 LU y eliminación de Gauss sin intercambios 2.5.2 LU y eliminación de Gauss con intercambios 2.6 Métodos para encontrar la factorización LU 2.7 LU para resolver el problema de mínimos cuadrados 2.8 Función LU en MATLAB® 2.9 Evaluación del algoritmo LU 2.10 Conclusiones

28 28 29 29 29 30 31 33 34 35 37

Capítulo 3

39

Factorización QR

39

3.1 Introducción 3.2 Definiciones y conceptos generales 3.3 Bases ortogonales y ortonormales 3.4 Creación de bases ortogonales y ortonormales: Gram-Schmidt. 3.5 Proyecciones ortogonales 3.6 Descomposición QR 3.6.1 Gram-Schmidt 3.6.2 Rotaciones de Givens 3.6.3 Reflexiones de Householder 3.7 Aplicación de la Factorización QR a mínimos cuadrados 3.8 Función QR en MATLAB®. 3.9 Evaluación del algoritmo QR 3.10 Conclusiones

39 40 43 44 45 47 51 51 53 54 56 57 63

Capítulo 4

65

Factorización AUDI

65

4.1 4.2 4.3 4.4 4.5 4.6 4.7

Introducción Definiciones y conceptos generales Factorización UD Aumentada (AUDI) Algoritmo AUDI-LS Algoritmo AUDI-ELS Evaluación del algoritmo AUDI Conclusiones

65 65 72 76 78 80 82

Capítulo 5

83

Norma L1

83

5.1 Introducción 5.2 Definiciones y conceptos generales 5.3 Aplicación de la norma L1 a identificación

ii

83 83 84

Contenido

5.4 Error de predicción final (fpe) para L1 5.5 Evaluación de la norma 1 para estimación de parámetros 5.6 Conclusiones

85 86 87

Capítulo 6

89

Validación de algoritmos

89

6.1 Introducción 6.2 Estructura del modelo 6.2.1 Estructura ARX. 6.2.2 Estructura ARMAX. 6.3 Identificación de los procesos 6.3.1 Sistema de tercer orden 6.3.1.1 Estructura ARX 6.3.1.2 Estructura ARMAX 6.3.2 Banco de nivel (RCN 100) 6.3.2.1 Estructura ARX 6.3.2.2 Estructura ARMAX 6.3.3 Secadora de pelo (Dryer) 6.3.3.1 Estructura ARX 6.3.3.2 Estructura ARMAX 6.3.4 Sistema de calentamiento (Heating System) 6.3.4.1 Estructura ARX 6.3.4.2 Estructura ARMAX 6.4 Conclusiones

89 90 91 92 93 93 93 99 101 103 107 109 109 113 115 115 119 121

Capítulo 7

127

Conclusiones y trabajos futuros

127

7.1 Conclusiones 7.2 Trabajos futuros

127 135

Referencias

137

Apéndice A

143

Descripción de Procesos

143

A.1 A.2 A.3 A.4

Sistema de tercer orden. Banco de nivel RCN 100. Secadora de pelo (Dryer) Sistema de calentamiento (Heating System)

143 143 144 145

Apéndice B

147

Selección del periodo de muestreo

147

B.1 Obtención del periodo de muestreo para el sistema de tercer orden. iii

147

Contenido

iv

Lista de figuras Figura i.1 Diagrama a bloques de uno de los métodos para la identificación de parámetros

2

Figura 1.1 Etapas para la identificación de un proceso [53].

15

Figura 1.2 GUI del toolbox de identificación

16

Figura 3.1 Rotación de Givens.

52

Figura 3.2 Matriz Householder H , veces el vector x .

53

Figura 4.1 Mínimos cuadrados vs Modelos múltiples de mínimos cuadrados [48].

68

Figura 4.2 Árbol de aplicación del estimador MMLS/AUDI [48].

82

Figura 6.1 Esquema de bloques del modelo general

90

Figura 6.2 Datos de entrada/salida para el sistema de tercer orden.

94

Figura 6.3 Respuesta al escalón y al impulso del sistema de tercer orden.

95

Figura 6.4 Respuesta en frecuencia del sistema de tercer orden.

95

Figura 6.5 Estimaciones paramétricas para el sistema de tercer orden para un modelo tipo ARX. 98 Figura 6.6 Estimaciones paramétricas para el sistema de tercer orden para un modelo tipo ARMAX.

101

Figura 6.7 Datos de entrada/salida del proceso banco de nivel RCN100.

103

Figura 6.8 Datos de entrada/salida después de eliminar tendencias del proceso banco de nivel RCN100. 103 Figura 6.9 Respuesta al escalón y al impulso del banco de nivel RCN100.

104

Figura 6.10 Respuesta en frecuencia del banco de nivel RCN100.

104

Figura 6.11 Estimaciones paramétricas del proceso banco de nivel RCN100 para un modelo tipo ARX. 107 Figura 6.12 Estimaciones paramétricas del proceso banco de nivel RCN100 para un modelo tipo ARMAX. 109 Figura 6.13 Datos de entrada/salida del proceso secadora de pelo.

110

Figura 6.14 Datos de entrada/salida después de eliminar tendencias del proceso secadora de pelo. 110 Figura 6.15 Respuesta al escalón y al impulso de la secadora de pelo.

111

Figura 6.16 Respuesta en frecuencia de la secadora de pelo.

111

Figura 6.17 Estimaciones paramétricas del proceso secadora de pelo para un modelo tipo ARX.113 Figura 6.18 Estimaciones paramétricas del proceso secadora de pelo para un modelo tipo ARMAX. 115 Figura 6.19 Datos de entrada-salida del sistema de calentamiento.

115

Figura 6.20 Datos de entrada/salida después de eliminar tendencias del proceso sistema de calentamiento.

116

v

Lista de figuras

Figura 6.21 Respuesta al escalón y al impulso de la secadora de pelo.

116

Figura 6.22 Respuesta en frecuencia de la secadora de pelo.

116

Figura 6.23 Estimaciones paramétricas del proceso sistema de calentamiento para un modelo tipo ARX. 119 Figura 6.24 Estimaciones paramétricas del proceso sistema de calentamiento para un modelo tipo ARMAX 121 Figura 7.1 Tiempos de ejecución de los algoritmos para una estructura ARX.

130

Figura 7.2 Tiempos de ejecución de los algoritmos para una estructura ARMAX

131

Figura A. 1 Sistema de tercer orden

143

Figura A. 2 Banco de nivel RCN100.

144

Figura A. 3 Secadora de pelo (dryer).

145

Figura A. 4 Sistema de calentamiento (heating system).

145

vi

Lista de tablas Tabla 3.1 Algoritmos para obtener la factorización QR mediante rotaciones de Givens.

61

Tabla 3.2 Algoritmos para obtener la factorización QR mediante transformaciones householder. 61 71

Tabla 4.1 Métodos de factorización para MMLS [48]. Tabla 4.2 Métodos de descomposición para la solucionar mínimos cuadrados[48].

72

Tabla 4.3 Relación entre la factorización QR y las factorizaciones LDL , LU y Cholesky[48].

73

Tabla 4.4 Algoritmo de mínimos cuadrados recursivos (RLS)[48].

74

Tabla 4.5 Algoritmo para obtener la factorización UD de Bierman[48].

75

Tabla 4.6 Algoritmo AUDI-LS recursivo[48].

77

Tabla 4.7 Algoritmo recursivo AUDI-ELS [48].

79

Tabla 6.1 Casos especiales más comunes de (6.1) [53].

90

T

Tabla 6.2 Porcentaje de desviación de los parámetros estimados para el sistema de tercer orden con estructura ARX. 95 Tabla 6.3 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el sistema de tercer orden con estructura ARX. 96 Tabla 6.4 Tiempo de ejecución de QR vs ARX, LS, AUDI y L1 para el sistema de tercer orden con estructura ARX. 96 Tabla 6.5 Porcentaje de desviación de los parámetros estimados para el sistema de tercer orden con estructura ARMAX. 99 Tabla 6.6 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el sistema de tercer orden con estructura ARXMAX. 99 Tabla 6.7 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para el sistema de tercer orden con estructura ARXMAX. 99 Tabla 6.8 ECM y EAM de los modelos estimados para el proceso banco de nivel RCN100 con estructura ARX. 105 Tabla 6.9 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el banco de nivel RCN100 con estructura ARX. 105 Tabla 6.10 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para el banco de nivel RCN100 con estructura ARX. 105 Tabla 6.11 ECM y EAM de los modelos estimados para el proceso banco de nivel con estructura ARMAX. 107 Tabla 6.12 Tiempo de ejecución de ARMAX, LS, LU, AUDI y L1 para el banco de nivel RCN100 con estructura ARMAX. 107

vii

Lista de tablas

Tabla 6.13 Tiempo de ejecución de ARMAX, LS, QR, AUDI y L1 para el banco de nivel RCN100 con estructura ARMAX. 107 Tabla 6.14 ECM y EAM de modelos estimados para el proceso secadora de pelo con estructura ARX. 111 Tabla 6.15 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para la secadora de pelo con estructura ARX.

111

Tabla 6.16 Tiempo de ejecución de ARX, QR, LU, AUDI y L1 para la secadora de pelo con estructura ARX.

111

Tabla 6.17 Desviación del error en los modelos estimados para el proceso secadora de pelo con estructura ARMAX. 113 Tabla 6.18 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para la secadora de pelo con estructura ARX.

113

Tabla 6.19 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para la secadora de pelo con estructura ARX.

113

Tabla 6.20 ECM y EAM de los modelos estimados para el proceso sistema de calentamiento con estructura ARX. 117 Tabla 6.21 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el sistema de calentamiento con estructura ARX. 117 Tabla 6.22 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para el sistema de calentamiento con estructura ARX. 117 Tabla 6.23 Desviación del error en los modelos estimados para el proceso sistema de calentamiento con estructura ARMAX.

119

Tabla 6.24 Tiempo de ejecución de ARMAX, LS, LU, AUDI y L1 para el sistema de calentamiento con estructura ARMAX. 119 Tabla 6.25 Tiempo de ejecución de ARMAX, LS, QR, AUDI y L1 para el sistema de calentamiento con estructura ARMAX. 119 Tabla 6.26 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos LU para una estructura tipo ARX. 122 Tabla 6.27 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos LU para una estructura tipo ARMAX. 123 Tabla 6.28 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos QR para una estructura tipo ARX. 123 Tabla 6.29 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos QR para una estructura tipo ARMAX. 123 Tabla 7.1 Criterios para la selección y validación de modelos.

128

Tabla 7.2 Orden máximo de cada algoritmo para cada proceso con una estructura ARX.

129

Tabla 7.3 Orden máximo de cada algoritmo para cada proceso con una estructura ARMAX.

130

viii

Nomenclatura Lista de símbolos y notación utilizada en el texto AT

Transpuesta de la matriz A .

A −1

Inversa de la matriz A .

A

Matriz de datos aumentada.

arg min f ( x)

Valor de x que minimiza f ( x)

D

Matriz diagonal.

C

Matriz de covarianzas.

C

Matriz de covarianzas aumentada.

.

()

⊥

Complemento ortogonal.

diag ( A)

Elementos de la diagonal de una matriz A .

E

Matriz de residuos.

Ω

Matriz de Givens.

e(t )

Perturbación en un tiempo t .

ε (t , θˆ )

Error de predicción.

H

Matriz de reflexión o transformación Householder.

j

Vector de datos, formado con los vectores entrada/salida.

J,V

Función de costo.

L

Matriz triangular inferior.

ix

Nomenclatura

. .

Norma de un vector. Norma de una matriz.

Po

Matriz de proyección ortogonal.

Q

Matriz ortogonal.

R

Matriz triangular superior.

rango( A)

Rango de una matriz A .

θˆ

Vector de parámetros estimados.

θ

Vector aumentado de parámetros estimados.

Θ

Vector de parámetros para el sistema aumentado.

yˆ (t | θ)

Salida estimada en un tiempo t .

Acrónimos ACM:

Augmented Covariance Matriz (Matriz de covarianzas aumentada).

AIM:

Augmented Information Matrix (Matriz de Información Aumentada).

AIC

Akaike Information Criteria (Criterio de información de Akaike).

ARARMAX :

AutoRegresive AutoRegresive Moving Average with eXternal input (Auto regresivo Auto regresivo con Media Móvil con entrada exógena).

ARARX:

AutoRegresive AutoRegresive with eXternal input (Auto regresivo Auto regresivo con entrada exógena).

ARMA:

AutoRegresive Moving Average (Auto regresivo con Media Móvil).

ARMAX:

AutoRegresive Moving Average with eXternal input (Auto regresivo de Media Móvil con entrada exógena).

x

Nomenclatura

ARX:

AutroRegresive with eXternal input (Auto regresivo con entrada exógena).

AUDI:

Augmented UD Identification (Identificación UD Aumentada).

AUDI-LS:

Augmented UD Identification-Least-Squares Aumentada con Mínimos Cuadrados).

AUDI-ELS:

Augmented UD Identification-Extended Least-Squares (Identificación UD Aumentada con Mínimos Cuadrados Extendidos).

BJ:

Box-Jenkins Model Structure (Estructura de modelo Box-Jenkins).

(Identificación

UD

DEGQRF: ELS:

Extended Least-Squares (Mínimos Cuadrados Extendidos).

IV:

Instrumental Variables (Variables Instrumentales).

FIR:

Finite Response Impulse (Respuesta al Impulso Finita).

FPE:

Final Prediction Error (Error de Predicción Final).

LAPACK:

Linear Algebra Package (Paquete de Álgebra Linear).

LS:

Least-Squares (Mínimos Cuadrados).

MATLAB:

Matrix Laboratory (Laboratorio de Matrices).

MIMO:

Multiple Input Multiple Output (Múltiples Entradas Múltiples Salidas).

MMLS:

Multiple Model Least-Squares (Mínimos Cuadrados de Modelos Múltiples).

OE:

Output Error Model Estructure (Estructura de modelo Output Error).

RLS:

Recursive Least-Squares (Mínimos Cuadrados Recursivos).

SISO:

Single Input Single Output (Una Entrada Una Salida).

SVD:

Singular Value Decomposition (Descomposición en Valores Singulares)

xi

Nomenclatura

TE:

Tiempo de Ejecución.

xii

Introducción El funcionamiento de los procesos industriales ha cambiado drásticamente en las últimas décadas. Esto se debe principalmente a la evolución de la tecnología de las computadoras. Como consecuencia la automatización de los procesos ha requerido un aumento en la productividad de algunos sectores industriales, obligando a la industria a adaptarse a las demandas de mercado y aumentar su nivel de competencia. Para lograr la eficiencia que los procesos industriales exigen actualmente ha sido necesario desarrollar nuevas técnicas que permitan maximizar la eficiencia de los procesos, desarrollando controladores de gran calidad, y maximizar la flexibilidad de los procesos con el menor ajuste de la planta [12]. Por lo tanto es primordial conocer el comportamiento dinámico del proceso, ya que para mejores controles se necesitan mejores modelos. La identificación de sistemas es el conjunto de teorías, métodos y algoritmos que permiten a partir de datos experimentales de entradas y salidas obtener un modelo matemático que represente la dinámica del sistema en cuestión [04]. Los modelos se utilizan en diferentes áreas como: bioingeniería, construcción, economía, meteorología, procesos químicos, por mencionar algunos. El campo en el cual dichos modelos se pueden utilizar es muy amplio, destacando aplicaciones como: control, supervisión, predicción, simulación, optimización, etc. La identificación de sistemas tiene sus raíces en técnicas estadísticas estándar y algunas de sus rutinas básicas igualmente conocidas como métodos estadísticos son mínimos cuadrados y máxima verosimilitud. La comunidad de control tomó una parte activa en el desarrollo y aplicación de estas técnicas básicas para sistemas dinámicos justo después del nacimiento de la “teoría de control moderna” a principios de los años 60’s. Åström y Bohlin en 1965 aplicaron la estimación de máxima verosimilitud a ecuaciones en diferencias lo cual dio pauta al surgimiento de un amplio rango de técnicas de estimación y parametrizaciones de modelos [51]. El objetivo principal de la identificación de sistemas es determinar un modelo matemático que represente la dinámica de cualquier proceso a partir de la recolección de datos entradasalida.

1

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Para determinar el modelo matemático de un proceso hay distintas formas de hacerlo, según se consideren modelos no paramétricos o modelos paramétricos [12]. La identificación puede ser aplicada tanto a sistemas SISO (single input-single output) como a sistemas MIMO (multiple input-multiple output), en este trabajo nos enfocaremos solo a modelos SISO, lineales e invariantes en el tiempo. Con la finalidad de comprender lo que es la identificación de sistemas en la Figura i.1 se muestra el diagrama a bloques de uno de los métodos existentes para realizar la identificación de parámetros para un proceso cualquiera. La señal de excitación u(t ) es inyectada a la entrada del proceso y del modelo, la señal de salida del proceso y (t ) será comparada con la señal de salida del modelo yˆ (t ) . La diferencia que existe entre las dos salidas, ε(t ) = y (t ) − yˆ (t ) , llamada error de estimación de predicción, es utilizada como señal de ajuste para corregir el modelo. Mediante un algoritmo de adaptación paramétrica se minimiza este error bajo una función pérdida con la finalidad de obtener un modelo mas aproximado.

Figura i.1 Diagrama a bloques de uno de los métodos para la identificación de parámetros

Existen varios criterios a partir de los cuales se pueden clasificar los métodos de identificación. En función del modelo obtenido pueden ser técnicas de identificación no paramétricas, de las que obtenemos modelos no paramétricos, y técnicas de identificación paramétricas, obteniéndose modelos paramétricos. Dentro de las técnicas de identificación no paramétricas podemos mencionar como las más importantes: •

Análisis de la respuesta transitoria, la cual se basa en la obtención de la respuesta del sistema a un impulso o escalón.

•

Análisis de correlación, genera la función de correlación entre las variables de interés.

2

Intorducción

•

Técnicas frecuenciales, las cuales estiman la respuesta frecuencial del sistema.

Estas técnicas de identificación son aplicables a sistemas lineales o linealizables, en la cual no se debe suponer ningún tipo de estructura para el modelo y los resultados son de tipo gráfico. Por otro lado en las técnicas de identificación paramétricas se debe escoger qué tipo de estructura determinada se utilizará para estimar el modelo. Los parámetros estimados se calculan minimizando el error ε(t ) , existente entre el modelo estimado y el proceso. En general podemos distinguir dos técnicas para el análisis de modelos paramétricos: •

Técnicas frecuenciales, las cuales minimizan el error entre la respuesta frecuencial real del proceso y la respuesta frecuencial del modelo.

•

Técnicas temporales, las cuales minimizan el error temporal, error de predicción o error de salida, entre el modelo y el proceso.

Ambas técnicas se pueden utilizar en la estimación de parámetros de modelos continuos o discretos. En este trabajo se analizan sistemas en tiempo discreto. Este trabajo de tesis se organiza de la siguiente manera. En el capítulo 1 se presentan los objetivos generales y particulares de trabajo y la revisión del estado del arte acerca de algunos de los algoritmos existentes para la estimación de parámetros. En el capítulo 2 se aborda el tema sobra la factorización LU, así como algunos de los métodos existentes para obtener tal factorización; también se muestra su aplicación a la solución del problema que presenta mínimos cuadrados y la etapa de evaluación. En el capítulo 3 se analiza la factorización QR, los métodos más confiables para lograr dicha factorización, además de la aplicación de dicha factorización a la solución del problema que presenta mínimos cuadrados y la etapa de evaluación. En el capítulo 4 se presenta la familia de algoritmos AUDI ó factorización UD Aumentada, la cual se basa en el modelo general de mínimos cuadrados, así como la utilización de esta factorización para solucionar el problema de mínimos cuadrados y también se presenta la evaluación del algoritmo. El capítulo 5 se enfoca en la utilización de la norma L1 como criterio de minimización del error y la evaluación de esta rutina de optimización para determinar parámetros de un sistema dado.

3

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

El capítulo 6 presenta los resultados obtenidos de la aplicación de los algoritmos mencionados anteriormente a la identificación de sistemas utilizando dos estructuras paramétricas, ARX y ARMAX. Para finalizar el documento se presenta un apartado que contiene las conclusiones así como los posibles trabajos futuros que se derivan de este trabajo de tesis.

4

Capítulo 1 Preliminares 1.1 Introducción En este capítulo se presenta la problemática que aborda este trabajo de tesis y la revisión del estado del arte sobre los algoritmos alternativos (LU, QR, AUDI, norma L1) para sustituir al algoritmo clásico de mínimos cuadrados aplicados a la identificación de sistemas. En la sección 1.2 se incluyen los antecedentes y la problemática asociada. En la sección 1.3, se muestra el marco teórico del algoritmo clásico utilizado en estimación de parámetros así como las alternativas planteadas anteriormente. En la sección 1.4 se presenta la hipótesis del trabajo de investigación. En la sección 1.5 se presentan el objetivo general y objetivos particulares de este trabajo. En la sección 1.6 se presenta la metodología que se siguió para desarrollar el trabajo de tesis. En la sección 1.7 se presentan los alcances y las limitaciones del trabajo. La sección 1.8 describe la condición que deben de cumplir las señales de excitación que se utilizan en tareas de identificación. La sección 1.9 trata del proceso de identificación que se seguirá en la etapa de evaluación de los algoritmos. Para finalizar en la sección 1.10 se presentan las conclusiones del capítulo.

1.2 Antecedentes y problemática asociada El algoritmo de mínimos cuadrados (LS) se ha logrado imponer como el estándar para la búsqueda de parámetros en el proceso de identificación de sistemas. A pesar de que los datos no estén bien filtrados o bien tomados este algoritmo es eficiente en la obtención de los parámetros debido a su robustez, estabilidad y convergencia a pesar del ruido [17].

5

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Por otro lado algunas de sus desventajas son que es un algoritmo lento cuando los datos a estimar son varios, además que tiende a quedarse en mínimos locales; es decir, no logra extraer de manera eficiente la información que proporcionan los datos [50]. Además de que LS involucra la inversión de una matriz y de esta manera está sujeto a problemas numéricos [44]. Por lo tanto se pueden definir dos problemas serios que existen con el estimador de mínimos cuadrados: •

Desempeño numérico muy pobre, sobre todo sobreparametrizado, que es el caso de la identificación.

cuando

el

modelo

es

•

No utiliza de manera eficiente la información contenida en la matriz de covarianzas [50].

Debido a esta problemática que presenta el algoritmo de mínimos cuadrados no siempre proporciona un modelo lo suficientemente representativo del sistema que se esté analizando, esto se debe a que en muchas ocasiones existen dinámicas que no están bien modeladas como consecuencia de que no se tiene un buen aprovechamiento de los datos. Por consiguiente es necesario analizar otras alternativas para realizar el proceso de identificación, como factorizaciones LU, QR, UD Aumentada ó AUDI por mencionar algunas, además del cambio de norma de L2 a L1 como criterio para minimizar el error.

1.3 Algoritmos para la estimación de parámetros 1.3.1 Mínimos cuadrados El problema de mínimos cuadrados se reduce a la resolución de un conjunto ecuaciones simultáneas lineales sobredeterminados de la forma  y1   a11 a12  y  a  2   21 a22  y3  =  a31 a32    ⋮ ⋮   ⋮  y  a a  m   m1 m 2

a13 a23 a33 ⋮ am 3

⋯ a1n   θ1  ⋯ a2 n  θ 2  ⋯ a3n  θ 3    ⋱ ⋮  ⋮  ⋯ amn  θ m 

(1.1)

o en forma compacta y = AT θ

donde: AT ∈ ℜ mxn θ ∈ ℜ nx1 y ∈ ℜ mx1

Es la matriz de coeficientes (pares de datos entrada/salida). Es el vector de parámetros desconocidos. Es el vector de observación (salida del sistema).

6

(1.2)

Capítulo 1:Preliminares

El algoritmo de mínimos cuadrados permite encontrar la solución

θˆ = θˆ1 ,θˆ2 ,...,θˆn 

que

minimice la suma de los cuadrados de la diferencia entre los datos observados y sus valores estimados, la cual es θˆ = arg min y − Aθ θ

2

(1.3)

donde • 2 denota la norma Euclidiana de un vector. De acuerdo al teorema de mínimos cuadrados de Gauss [48] la solución a la ecuación (1.2) se puede encontrar mediante la ecuación normal −1

 AT A  θ = AT y

(1.4)

θˆ = ( AT A) −1 AT y

(1.5)

Conceptualmente, la solución a

y la función de costo que minimiza el error esta dada por J=

1 2 ∑ [ y (t ) − yˆ (t )] 2

(1.6)

donde: yˆ (t ) = AT (t )θ(t )

(1.7)

es la salida estimada. Si observamos la ecuación (1.5) vemos que el algoritmo requiere una inversión de matriz como se mencionó anteriormente, lo cual resulta muy complejo numéricamente y vuelve muy lento el cálculo computacional cuando la cantidad de datos es demasiado grande. Existen alternativas al algoritmo clásico de mínimos cuadrados para la optimización mediante algoritmos de identificación paramétrica que permiten un mejor aprovechamiento de la información contenida en el par de datos entrada-salida. Estas alternativas aun no han reportado su aporte en la literatura en el caso de la identificación, por lo que en este trabajo se estudian y se implementan tales alternativas como herramientas de optimización para la identificación de sistemas. Los algoritmos alternativos estudiados son: • • • •

Factorización LU Factorización QR Factorización AUDI o UD aumentada Norma L1

7

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

1.3.2 Algoritmos alternativos para la estimación de parámetros. 1.3.2.1 Factorización LU La utilización de la factorización o descomposición LU, que también es conocida como eliminación Gaussiana, data de varias décadas, y su aplicación ha sido muy variada durante todo este tiempo. Generalmente, esta factorización se utiliza en la solución de sistemas de ecuaciones lineales “sparse”, las cuales se presentan principalmente en sistemas eléctricos y sistemas de potencia [07], [16], [26], [27], [57]. Teniendo como objetivo principal la reducción del tiempo de procesamiento en la ejecución de operaciones. Otra de su aplicación ha sido en la descomposición de una matriz de coeficientes en matrices triangulares inferior y superior, para la realización de filtros digitales 2-D [05]. Ahora bien, la factorización LU en términos generales consiste en lo siguiente. Dado un sistema no homogéneo (1.8)

Ax = b

si la matriz A es de rango completo se puede factorizar de tal manera que A = LU

(1.9)

LUx = b

(1.10)

podemos representar al sistema como

donde: L y U son matrices triangulares inferiores y superiores respectivamente.  l11 0   l21 l22 L =  l31 l32  ⋮ ⋮ l  n1 ln 2

0 0 l33 ⋮ ln3

⋯ 0  u11 u12   ⋯ 0  0 u22 ⋯ 0 , U =  0 0   ⋱ ⋮  ⋮  ⋮   ⋯ lnn  0  0

u13 u23 u33 ⋮ 0

⋯ ⋯ ⋯ ⋱ ⋯

u1n   u2 n  u3n   ⋮  unn 

(1.11)

de manera que si Ux = z , el sistema puede representarse como Lz = b

por lo tanto, el sistema a resolver mediante esta factorización está dado por

8

(1.12)

Capítulo 1:Preliminares

l11 z1

= b1

l21 z1 + l22 z2 l31 z1 + l32 z2 + l33 z3

= b2 = b3

l41 z1 + l42 z2 + l43 z3 + l44 z4

= b4

⋮

⋮

(1.13)

⋮

así, el sistema se resuelve de forma progresiva y para determinar los valores de x se resuelve el sistema de manera regresiva. Existen diferentes formas de obtener la factorización LU y algunas de ellas se analizarán en el Capítulo 2.

1.3.2.2 Factorización QR Otra alternativa para solucionar el problema que presenta mínimos cuadrados en (1.3) es la descomposición ortogonal QR, la cual permite obtener mucho más información que el algoritmo tradicional de mínimos cuadrados, y utilizando una descomposición adecuada se puede obtener una estructura simple y fácil de implementar [48]. En [11] se desarrolla una librería, con base en el uso de la factorización QR en forma recursiva, que permite usar la memoria jerárquica eficientemente en sistemas de multiprocesamiento. La factorización LU se utiliza ampliamente en la solución de sistemas lineales “sparse” debido a que el tiempo de cálculo y el uso de memoria es menor, a pesar de que la factorización QR presenta una mejor estabilidad numérica en [41] se plantea una solución a este tipo de sistemas utilizando la factorización QR. En [06] se desarrolla un algoritmo basado en la factorización QR utilizando la estrategia de pivote local, teniendo los mismos resultados en menor tiempo de ejecución al igual que utilizar la estrategia de pivote global. En [39] se propone un algoritmo en base a la factorización QR para descomponer un conjunto de matrices las cuales tienen columnas en común. Aunque existen diversas aplicaciones de la factorización QR y LU no se ha reportado el uso de las mismas en tareas de identificación. Análogamente a la factorización LU, dado el sistema (1.8) si la matriz A es de rango completo puede proponerse la factorización A = QR

donde:

9

(1.14)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

R es una matriz triangular superior definida positiva, es decir, rij > 0 R

Q es una matriz ortogonal, es decir QT Q = I

(1.15)

Utilizando la expresión anterior, el sistema puede representarse como

Rx = QT b

(1.16)

Si QT b = d , el sistema queda representado por Rx = d

(1.17)

de la misma forma que en la factorización LU el sistema se resuelve de manera progresiva y se obtienen los valores de x de manera regresiva. Existen varios métodos para obtener tal factorización que se explican a detalle en el Capítulo 3.

1.3.2.3 Factorización AUDI o UD aumentada

También existe otro tipo de factorización para la identificación de sistemas, el cual consiste en obtener múltiples modelos con base en el algoritmo de mínimos cuadrados. Este método tiene una estructura recursiva y produce múltiples modelos de orden 1 a un valor máximo “n” especificado por el usuario [44]. Una de las características que distingue a este algoritmo o estructuras es que son mucho menos vulnerables con respecto al cálculo numérico, además que evita la sobre-parametrización [50]. Su forma recursiva aplicada a la identificación de sistemas se le conoce como Identificación UD (AUDI). Una de las aplicaciones que se encuentra en la literatura se muestra en [47], la cual consiste en el desarrollo de un método para detectar problemas de identificabilidad de parámetros debido a excitación no-persistente, sobre-parametrización y/o retroalimentación en la salida del sistema identificado [20]. La factorización AUDI de manera resumida consiste en dado un sistema y (t ) + a1 y (t − 1) + ... + an y (t − n) = b1u (t − 1) + ... + bn u (t − n) + u (t ) + v(t )

(1.18)

en el formato del algoritmo de mínimos cuadrados tenemos que y (t ) = XT (t )θ(k ) + v (t )

yX ()T= + kv onde: dtΘ 10

(1.19)

Capítulo 1:Preliminares

n

es el orden del modelo.

u(t ) y y (t )

Son las entradas y salidas del sistema respectivamente.

v (t )

es un ruido blanco o cualquier tipo de perturbación.

θˆ (t )

es el vector de parámetros desconocidos.

este método se basa en los principios del algoritmo de mínimos cuadrados, pero con una implementación más eficiente [50]. El vector de datos aumentado está definido por Ψ = [−y (t − n), u(t − n), ⋯ , − y (t − 1), u (t − 1), − y (t )]T

(1.20)

y el vector de parámetros desconocidos es θ = [an , bn , ⋯ a1 , b1 , 1]T

(1.21)

la matriz de covarianzas es t

c(t ) = [∑ ψ ( j )ψ T ( j )] = U(t )D(t )UT (t )

(1.22)

j =1

donde: U = U(t ) es la matriz de parámetros de forma triangular superior, la cual contiene

los parámetros de los modelos estimados de orden 1 hasta “n”. D = D−1 (t ) es la matriz que contiene la suma de los cuadrados de los residuos (“loss function”). En el Capítulo 4 se presentan a detalle la familia de algoritmos AUDI, así como su evaluación correspondiente. 1.3.2.4 Norma L1 La norma L1 generalmente tiene su aplicación en el diseño de controladores robustos. En [35] se presenta un método para sintonizar controladores que logren un desempeño robusto globalmente óptimo mediante el uso de la norma L1. En [40] la norma L1 se usa como técnica alternativa para el cálculo de subespacio robusto, donde la descomposición en valores singulares (SVD) es el algoritmo estándar. Una de las aplicaciones de esta norma en un ámbito diferente es en el rubro de economía [03], teniendo como objetivo minimizar la suma de las desviaciones absolutas y realizar

11

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

una comparación con la norma que se utiliza comúnmente, la norma L2, como el método de ajuste de la regresión en problemas económicos. Por otro lado en [21] se aplica la norma L1 para la identificación y validación de un modelo. El modelo obtenido utilizando este criterio se compara con otro modelo que se obtuvo utilizando el criterio clásico que se basa en la norma L2; el caso de estudio es un experimento de laboratorio, el cual consiste en propagar un sonido a través de un ducto acústico [17]. La utilización de la norma L1 consiste en encontrar una solución, para el caso de identificación, que minimice la suma del valor absoluto de la diferencia de los datos observados, es decir θˆ = arg min y − Xθ θ

(1.23)

En el Capítulo 5 se aborda la aplicación de este criterio en la estimación de parámetros mediante el uso de una rutina de optimización.

1.4 Hipótesis El desarrollo de algoritmos de optimización como alternativas al algoritmo clásico de mínimos cuadrados aplicados a la identificación de sistemas, permitirá un aprovechamiento más eficiente de los datos entrada-salida para determinar el modelo matemático de un sistema dinámico.

1.5 Objetivos 1.5.1 Objetivo particular Contribuir en los algoritmos de optimización aplicados a la identificación de modelos para mejorar el entendimiento de la dinámica de los procesos, además de realizar un esfuerzo de integración de los diferentes desarrollos en algoritmos de optimización puntuales existentes que no han sido aplicados a la identificación. 1.5.2 Objetivos particulares Los objetivos particulares son los siguientes: 1. Revisar el estado del arte mediante el cual se conozcan y analicen algunas de las alternativas existentes al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. 2. Programar y validar cada una de las alternativas seleccionadas.

12

Capítulo 1:Preliminares

3. Realizar una serie de identificaciones de procesos conocidos, de los cuales se conozca su modelado para comprobar los programas elaborados. 4. Comparar los resultados obtenidos con los algoritmos de procesos conocidos aplicados a la identificación de sistemas como ARX, ARMAX.

1.6 Metodología Para alcanzar los objetivos planteados, el trabajo se divide es las siguientes etapas: 1. Selección de algoritmos. Apoyado en la revisión de literatura dentro del contexto de algoritmos de estimación de parámetros, se seleccionaron 4 algoritmos alternativos: Factorización LU, factorización QR, factorización UD Aumentada (AUDI) y la norma L1. 2. Programación Programación en MATLAB® de cada uno de los algoritmos seleccionados. 3. Evaluación Realización de pruebas numéricas para verificar que los algoritmos seleccionados obtengan los datos correctos y así comprobar la programación de cada uno. 4. Validación Se utilizaron datos obtenidos experimentalmente o mediante simulación proporcionados por la base de datos Daysi1 para validar la aplicación de cada uno de los algoritmos.

1.7 Alcances y limitaciones del trabajo 1. Probar como mínimo cuatro métodos diferentes para la optimización. 2. Programar en MATLAB® los métodos de optimización seleccionados. 3. Comparar los métodos programados con datos experimentales de procesos conocidos con los que cuenta el CENIDET en el laboratorio de control o de procesos reportados en la literatura.

1

Daysi: Base de datos que proporciona datos experimentales y/o de simulación de algunos procesos. 13

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

1.8 Señal de excitación persistente Dentro del contexto de identificación de sistemas, un aspecto muy importante es la selección de la señal de excitación de entrada al sistema a identificar. Para ello es importante cumplir con la condición de unicidad, la cual indica que el mínimo obtenido es único y esta dado por (1.5), siempre y cuando la matriz AT A sea no singular, es decir, de rango completo. Esta condición se conoce como condición de excitación [23]. Definición 1.1 Excitación persistente

Una señal cuadrada u es de excitación persistente (PE, Persistently Exciting) de orden n , si la matriz AT A es positiva definida.

1.9 Proceso de identificación 1.9.1 Etapas a seguir para la identificación de un modelo Las etapas que se tienen que seguir para la identificación de un modelo son las siguientes: 1. Diseño del experimento y adquisición de datos. En esta etapa principalmente se tiene que decidir que tipo de señales de excitación se van a utilizar, el periodo de muestreo adecuado para la adquisición de datos y la cantidad de datos necesarios. 2. Observación y mejora de la calidad de los datos. Antes de utilizar los métodos de estimación se necesita observar y reparar los datos erróneos, eliminar “offsets” y tendencias, etc. 3. Determinación de la estructura del modelo. En esta etapa es necesario definir los tipos de modelos a utilizar, que pueden ser: continuos o discretos, tipos de ruido, lineales o no lineales, regresiones, redes neuronales, etc. Y dependiendo de la estructura del modelo, utilizar un procedimiento adecuado para determinar el orden del modelo. 4. Estimación de los parámetros. Esta etapa tiene generalmente tiene una vinculación con la etapa anterior; en ella se presenta el problema de decidir qué método o métodos de estimación de parámetros se va utilizar para calcular el valor de los mismos. Normalmente se puede elegir entre técnicas en el dominio temporal o técnicas en el dominio frecuencial. 5. Validación del modelo. Es la etapa en la que se cuestiona si el modelo obtenido representa la dinámica del proceso estudiado. En él se debe definir un criterio para

14

Capítulo 1:Preliminares

evaluar la calidad. Comúnmente se tiene un conjunto de modelos candidatos y se debe elegir el mejor mediante un criterio específico. En la Figura 1.1 se muestran las etapas para la identificación que se describieron anteriormente, en la cual los rectángulos son dependientes directamente del computador y los óvalos son en general responsabilidad del usuario.

Figura 1.1 Etapas para la identificación de un proceso [53].

Actualmente existen programas comerciales de ayuda al ingeniero en las etapas de identificación, durante las etapas de estimación de parámetros y evaluación de las propiedades del modelo estimado. Dentro de esos programas la herramienta de trabajo para nuestro caso es MATLAB® 7.0 que contiene el Toolbox de Identificación [33]. El Toolbox de Identificación permite construir modelos matemáticos de un sistema dinámico a partir de datos medidos. La validación de los modelos estimados se basa en la cercanía que tenga la salida del modelo estimado con respecto a la salida de los datos medidos. El Toolbox de Identificación de MATLAB® 7.0 contiene las técnicas comunes para la estimación de parámetros para todo tipo de modelos lineales. También permite examinar las propiedades de los modelos y verificar si son buenos modelos. Si se asume que el modelo a estimar es lineal, con el Toolbox de identificación se puede obtener de manera directa su respuesta al impulso o al escalón utilizando análisis de correlación o su respuesta infrecuencia usando análisis espectral. Este Toolbox contiene una interfaz gráfica que permite al usuario una utilización más completa de las diversas funciones que contiene para la estimación de modelos lineales. En la Figura 1.2 se muestra una imagen de la GUI de identificación.

15

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Figura 1.2 GUI del toolbox de identificación

En términos generales la principal información que se presenta en la ventana principal se puede dividir como sigue: Tablas de trabajo: •

Una tabla acerca de los conjuntos de datos disponibles, cada uno se representa con un icono.

•

Una tabla referente a los modelos que se han creado, cada uno se representa con un icono.

Estas tablas se conocen como “Data Board” y “Model Board”. El usuario puede introducir un conjunto de datos de tres formas diferentes: •

Abriendo sesiones que se guardaron con anterioridad.

•

Importando los datos desde el espacio de trabajo de MATLAB®.

•

Creando los modelos mediante la eliminación de tendencias, filtrado, etc., de otro conjunto de datos que se encuentren disponibles en el “Data Board”.

Las importaciones se manipulan mediante el menú Data, este menú proporciona las diferentes formas de importar un conjunto de datos: •

Datos en el dominio del tiempo.

16

Capítulo 1:Preliminares

•

Datos en el dominio de la frecuencia.

•

Objeto de datos creado en la estación de trabajo “workspace”.

•

Ejemplo que contiene el toolbox para fines explicativos.

Mientras que la creación de nuevos conjuntos de datos se manipulan mediante el menú Preprocess. Ese menú proporciona rutinas para el pretratamiento del conjunto de datos como •

Eliminación de desplazamientos (offsets).

•

Eliminación de tendencias.

•

Filtrado

•

Remuestreo

•

Etc.

Las importaciones se manipulan mediante el menú Models, este menú proporciona la forma de importar modelos que no se generaron con la GUI de identificación. Mientras que los diferentes esquemas de estimación se manipulan mediante el menú Estimate. El cual proporciona las siguientes diferentes técnicas para estimar modelos: •

Estimación de modelos paramétricos, utilizando las estructuras de identificación paramétrica: ARX, ARMAX, OE, BJ.

•

Estimación de modelos de procesos representados en modelos de función transferencia.

•

Estimación de modelos utilizando el análisis espectral.

•

Estimación de modelos utilizando el análisis de correlación.

Básicamente son los rubros principales con los que cuenta la GUI de identificación para la estimación de modelos. 1.9.2 Selección y validación de la estructura del modelo. Otro de los aspectos muy importantes además de los métodos de estimación de parámetros es el de encontrar la estructura apropiada para el modelo (ver Tabla 6.1). Para ello es necesario entender los métodos de estimación y tener conocimiento del proceso que se va a identificar. Una vez que se determina la estructura del modelo, los métodos de estimación que se utilicen proporcionan un modelo particular de dicha estructura. Por consiguiente cuando ya

17

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

se tiene identificado el modelo la siguiente etapa es cuestionarse si el modelo obtenido es bueno, por lo que se hecha mano de las técnicas de validación. 1.9.2.1 Selección de la estructura del modelo Para poder abordar el problema de identificación es necesario tener en mente tres aspectos importantes: 1. Elegir el tipo de conjunto de modelos. 2. Selección del orden de los modelos. 3. Métodos de identificación. La selección de tipos de modelos se refiere a buscar entre diferentes clases de modelos lineales o no lineales, o entre modelos de entrada-salida, caja negra o parametrizados físicamente, entre otros. La búsqueda del tipo de modelo es un tanto subjetiva y habitualmente es un compromiso entre diferentes factores, como la calidad del modelo y el esfuerzo de modelado. Para obtener resultados aceptables con pocos parámetros es necesario conocimiento previo del proceso y un poco de intuición. La forma de evaluar al error y su minimización influye mucho en el esfuerzo de programación y el tiempo de cálculo. Las regresiones lineales son los modelos más simples y en los que la minimización del criterio es robusta. La selección del orden del modelo implica la selección del número de parámetros de una determinada clase de modelos. Por lo general este paso requiere de alguna evaluación del conjunto de datos. La selección de método de identificación depende de la clase de modelo elegido, por ejemplo un sistema lineal puede representarse de diferentes formas, por su respuesta transitoria o frecuencial o modelo paramétrico. Es claro que se puede transformar unas en las otras, pero estas transformaciones están muy mal condicionadas en el contexto de que pequeños errores en alguna representación pueden producir un gran error en la otra. Es decir la búsqueda del modelo está relacionada con el objetivo de la identificación.

1.9.2.2 Análisis preeliminares de los datos para la selección de la estructura del modelo Ciertas evaluaciones de los datos permiten producir una estructura del modelo. En el caso de los métodos paramétricos, el problema consiste en la selección del orden del modelo. El

18

Capítulo 1:Preliminares

orden de un sistema lineal de puede estimar de diferentes formas. Basado en el análisis de datos preliminares se pueden mencionar tres categorías. 1. La evaluación del análisis espectral estimado proporciona información sobre el rango de frecuencias de interés del proceso y el nivel de ruido, por lo tanto no ayuda a diseñar correctamente la señal de entrada y los filtros. En cambio la estimación no paramétrica de la función de transferencia proporciona más información sobre los picos de resonancia y el retardo de fase. Todo esto proporciona una orientación sobre el orden requerido para una adecuada descripción de las dinámicas del sistema. 2. Evaluar el rango de la matriz de covarianza es un buen método par estimar el orden del sistema. 3. El método de correlación también permite obtener información acerca de qué variables incluir en la estructura del modelo. 1.9.2.3 Criterios para la comparación de la estructura de los modelos Una aproximación natural para determinar la estructura del modelo es simplemente probar diferentes estructuras y comparar el resultado entre ellas. Algunos de los métodos de prueba más comunes son: •

Prueba de la función error o pérdida

•

Prueba de cancelación de polos y ceros

•

Prueba del residuo

1.9.2.3.1 Prueba de la función error o pérdida. Este método es el más simple ya que consiste en mirar la función pérdida V (θˆ ) ya que ésta se encuentra directamente vinculada con el orden del modelo. Cuando incrementa el orden, la función pérdida decrece hasta que se mantiene constante o cambia lentamente. Otros métodos se basan en pruebas estadísticas de la función pérdida o en la evaluación de diferentes criterios que tiene en cuenta la complejidad del modelo. a) Un método estadístico para evaluar si la función pérdida disminuye significativamente cuando incrementa el orden del modelo es el F-test. Si la reducción de la función es significativa cuando el número de parámetros aumenta de p1 a p2 se utiliza:

19

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

t=

V −1 (θ ) − V 2 (θ ) N − p2 V 2 (θ ) p1 − p2

(1.24)

Esta cantidad tiene distribución F [ p1 − p2 , N − p2 ] , entonces para un suficiente número de muestras N , se cumple que: t > Fα [ p1 − p2 , N − p2 ] entonces se cumple la hipótesis que la reducción de la función pérdida es significativa con el aumento de parámetros y por lo tanto el nuevo parámetro es aceptado con un nivel de confianza de 1 − α . b) El criterio de análisis de complejidad puede llevarse acabo simplemente agregando a la función pérdida un término extra que penalice la complejidad del modelo. En este caso se selecciona el mejor modelo que minimice al criterio. El criterio de información de Akaike (AIC) disminuye con la función pérdida y aumenta con el número de parámetros, se define como AIC ( p) = N log V (θ ) + 2 p

(1.25)

donde: V (θ ) =

1 N

N

∑ε

2

(t , θ )

(1.26)

t =1

c) Otro criterio propuesto por Akaike es el error de predicción final (FPE (Final Predictor Error)), que se define como FPE ( p ) =

N+p V (θ ) N−p

(1.27)

Una propiedad atractiva tanto del criterio AIC y FPE es que el orden del modelo se determina como resultado del valor mínimo del criterio y no es necesario evaluarlo en función de unos niveles de confianza como es el caso del método estadístico.

1.9.2.3.2 Prueba de cancelación de polos-ceros Esta prueba se basa en que la suposición de que el orden del modelo estimado es superior al orden del proceso real, y éstos originan pares de polos-ceros muy próximos que pueden ser cancelados.

1.9.2.3.3 Prueba de residuo Es el resultado de comparar la salida con la predicción del modelo, es decir ε(t , θˆ ) = y (t ) − yˆ (t , θ)

20

(1.28)

Capítulo 1:Preliminares

La cual representa una forma de representar las diferencias entre las variables observadas y el comportamiento del modelo estimado. En general los residuos deben de ser blancos o aproximadamente blancos, y no correlados con la entrada, si el modelo es una buena representación del sistema.

1.9.2.4 Validación del modelo Una vez identificado del modelo la pregunta crucial es si éste modelo es suficientemente adecuado para los objetivos considerados. La validación permite comprobar si el modelo identificado representa la dinámica del sistema real, tomando en cuenta las limitaciones de los métodos de identificación y el objetivo final. El problema se puede plantear de diferentes maneras 1. El modelo, ¿Satisface suficientemente bien los datos observados? 2. ¿Es el modelo suficientemente bueno para los requerimientos propuestos? 3. ¿Describe el modelo al sistema? Para responder a estas preguntas es necesario confrontar el modelo con más información, como conocimientos a priori, datos experimentales y experiencia, sobre el sistema real. Algunos de los criterios típicos a la hora de descartar o elegir unos modelos respecto a otros son: a) Validación en base a la aplicación del modelo Debido a que en la práctica es imposible determinar si un modelo responde exactamente a un sistema real, solo es necesario comprobar que el modelo es capaz de resolver el problema para el cual ha sido hallado (simulación, predicción, diseño de un controlador, etc.). b) Comprobación de parámetros físicos Para una determinada estructura parametrizada en función de magnitudes físicas, un método importante de validación es comparar el valor estimado de dichos parámetros y el que sería de esperar mediante el conocimiento previo que se tiene de la planta. c) Coherencia en el comportamiento de entrada-salida Para determinar si el comportamiento de entrada-salida esta lo suficientemente caracterizado, puede ser necesario recurrir a diferentes métodos de investigación y

21

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

comparar los resultados obtenidos. Por ejemplo: diagramas de Bode, el método de las variables instrumentales y/o por análisis espectral. d) Reducción de modelo Si una reducción en el orden del modelo no produce alteraciones apreciables en el comportamiento de entrada-salida del mismo, entonces el modelo original era innecesariamente complejo. e) Intervalos de fiabilidad de parámetros Otro método para determinar si el modelo bajo estudio tiene demasiados parámetros consiste en comparar los parámetros estimados con su desviación estándar. Si el intervalo de confianza de un parámetro contiene valor cero, se debe considerar eliminar dicho parámetro. f) Simulación Otro procedimiento muy habitual que puede ser considerado como otra técnica de validación consiste en simular el modelo con un conjunto de entradas distintas a las utilizadas para identificación, y comparar la respuesta del modelo con la obtenida del sistema real.

1.10 Conclusiones Este trabajo está enfocado hacia la evaluación comparativa de algunas de las alternativas existentes al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. La problemática puede abordarse mediante el estudio de las alternativas anteriormente planteadas las cuales solo se han estudiado en un nivel comparativo con respecto a las propiedades que ofrecen cada una de ellas respecto al algoritmo de mínimos cuadrados. El reto que presenta este trabajo es utilizar las alternativas de optimización al algoritmo de mínimos cuadrados para la identificación de sistemas. Las alternativas que se describieron anteriormente se validarán principalmente con procesos reales y/o procesos sintéticos de los cuales se conozca su estructura, con la finalidad de realizar una comparación entre las alternativas propuestas y el algoritmo tradicional de mínimos cuadrados. Para la identificación del sistema en cuestión, se seguirán los pasos presentados en la Sección 1.9. Los modelos resultantes de cada algoritmo serán elegidos tomando en cuenta los criterios de selección y validación que se presentaron en la Sección 1.9.2. En los capítulos posteriores se analizan a detalle cada uno de los algoritmos.

22

Capítulo 2 Factorización LU 2.1 Introducción La factorización o descomposición LU es el primero de los cuatro algoritmos alternativos que se estudian en este trabajo para sustituir al algoritmo clásico de mínimos cuadrados aplicados a la identificación de sistemas. En este capítulo se divide de la siguiente manera: La sección 2.2 contiene algunos conceptos básicos necesarios para la comprensión de operaciones con matrices. La sección 2.3 presenta el proceso de eliminación de Gauss con intercambio, la eliminación Gauss-Jordan y la obtención de la inversa de una matriz mediante eliminación de Gauss. La sección 2.4 describe la existencia de soluciones en sistemas de ecuaciones lineales. En la sección 2.5 se presenta como obtener la factorización LU utilizando la eliminación de Gauss. La sección 2.6 aborda los dos métodos más comunes para obtener la factorización LU, Eliminación Gaussiana y métodos directos: Doolittle, Crout y Cholesky. En la sección 2.7 se presenta la aplicación de la factorización LU para resolver el problema de mínimos cuadrados. En la sección 2.8 se describe la función lu que contiene MATLAB®. La sección 2.9 muestra la etapa de evaluación del algoritmo LU, donde se comprueba que los diferentes métodos realizar una factorización correcta. Por último la Sección 2.10 presenta las conclusiones del capítulo.

2.2 Definiciones y conceptos Generales En este apartado se establecen los conceptos y definiciones necesarias para el análisis de la factorización LU. Estas definiciones y conceptos generales se tomaron de [02] y [43].

23

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Definición 2.1 Matriz cuadrada

Una matriz cuadrada A de (m, n ) con m = n ; los elementos (i , i ) para 1 ≤ i ≤ m forman la diagonal principal de la matriz A . Definición 2.2 Matriz aumentada.

La matriz separada

[A

b ] es la matriz aumentada del sistema de ecuaciones que se

describe como Ax = b . Definición 2.3 Operaciones de renglón

Una operación elemental de renglón (a veces llamada solo operación de renglón) en una matriz A es cualquiera de los tres tipos de operaciones siguientes: •

Intercambio de dos renglones de A .

•

Reemplazo de un renglón r de A por cr , donde c es una constante, para algún número de c ≠ 0 .

•

Reemplazo de un renglón r1 de A por la suma r1 + cr2 de ese renglón y el múltiplo de otro renglón r2 de A .

Definición 2.4 Matriz elemental.

Una matriz elemental, mxn es una matriz que se produce aplicando exactamente una operación elemental de renglón a I m . E ij es la matriz elemental que se obtiene mediante el intercambio de los renglones i -ésimo y j -ésimo de I m . Ei (c ) es la matriz elemental obtenida multiplicando el i -ésimo renglón de I m por c ≠ 0 . Ei (c ) es la matriz elemental que se obtiene sumando c veces el j -ésimo renglón de I m , donde i ≠ j . Basándonos en las definiciones Definición 2.3 y Definición 2.4 tenemos el siguiente Teorema 2.1 Operaciones de renglón y matrices elementales

Suponga que E es una matriz elemental (m, n ) que se produce aplicando una operación elemental de renglón a I m y que A es una matriz (m, n ) . Entonces EA es la matriz que resulta de la aplicación de esa misma operación elemental de renglón en A .

24

Capítulo 2:Factoriozación LU

Teorema 2.2 Matrices elementales y no singularidad.

Cada matriz elemental es no singular, y su inversa, por sí misma es una matriz elemental. De modo más preciso: E ij −1 = E ji (= E ij ) E i (c )−1 = E i (1/ c )

1 con c ≠ 0 2

(2.1)

E ij (c )−1 = E ij ( −c ) con i ≠ j

Definición 2.5 Matriz de permutación.

Una matriz de permutación Pmxm es cualquier matriz que resulta de permutar, es decir, reordenar el orden de los renglones de I m . De modo más preciso: Cada renglón de P contiene exactamente un elemento diferente de cero, es decir, 1; y cada columna de P contiene exactamente un elemento diferente de cero, es decir, 1. Teorema 2.3 Matriz de permutación

Sea P una matriz de permutación. Entonces: •

Para cualquier A , PA puede obtenerse a partir de A permutando los renglones de A exactamente como se permutaron los renglones de I para obtener P .

•

P es no singular, y P −1 = PT : PP T = P T P = I .

Definición 2.6 Matriz tridiagonal.

Una matriz tridiagonal es una matriz cuadrada cuyos elementos que no están en la diagonal principal son iguales a cero. Por diag(d1 ,..., dm ) se entiende que la matriz diagonal mxm cuyos elementos (i , i ) equivalen a di para 1 ≤ i ≤ m . Definición 2.7 Matriz triangular inferior.

Una matriz triangular inferior L de mxn satisface L

ij

= 0 si i < j , para 1 ≤ i ≤ m y

1≤ j ≤n.

Definición 2.8 Matriz triangular superior.

Una matriz triangular superior U de mxn satisface U 1≤ j ≤n.

25

ij

= 0 si i > j , para 1 ≤ i ≤ m y

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Definición 2.9 Matriz triangular inferior o superior unidad T.

Es una matriz triangular inferior (o superior) que satisface T

ii

= 1 para 1 ≤ i ≤ min(m, n ) .

Eliminación de Gauss para matrices aumentadas en general Sea el sistema: a11 x1 a21 x1 ⋮

am1 x1

+ a12 x2 + a22 x2 ⋮ ⋮ + am 2 x2

+ + ⋮ +

⋯ +

a1 n xn

⋯ +

a2 n xn

⋱ ⋮

⋮

⋯ + amn x n

= b1 = b2 ⋮ ⋮ = bm

(2.2)

que consiste de m ecuaciones lineales con n incógnitas será equivalente a Ax = b para A de (m, n ) con A

ij

= aij , x nx 1 con x

j

= x j , y bmx 1 con b i = b i . La eliminación de

Gauss para matrices es posible realizarla en la matriz aumentada [ A

b] .

La eliminación Gaussiana se puede llevar a cabo mediante dos formas distintas:

•

Eliminación de Gauss con intercambios

•

Eliminación de Gauss sin intercambios

2.3 Eliminación de Gauss con intercambio (con pivote) Para poder llevar a cabo la eliminación Gaussiana con intercambios es necesario seguir el siguiente procedimiento: 1. Sea j = 1 y r1 = 1 ; utilice el renglón rj para eliminar en la columna j de la matriz aumentada actual como se indica en los pasos del 2 al 6. 2. Seleccione un renglón de entre los renglones numerados rj , rj + 1,..., m para su uso en la eliminación de la columna j ; llame a este renglón i , de modo que el elemento (i , j ) -llamado pivote- en la matriz aumentada actual sea diferente de cero. Si no

hay elementos diferentes de cero en esta parte inferior de la columna j , entonces no se requiere de eliminación; ponga rj +1 = rj para utilizar el mismo renglón y salte directamente al paso 6.

26

Capítulo 2:Factoriozación LU

3. Intercambie los renglones i-ésimo y rj -ésimo. 4. Reemplace este nuevo renglón rj -ésimo por sí mismo dividido entre el pivote (su elemento diferente de cero en su columna j-ésima). 5. Utilice este nuevo renglón rj -ésimo para eliminar los elementos en la columna jésima en los renglones rj , rj + 1,..., m . Ponga rj +1 = rj + 1 para utilizar el siguiente renglón. 6. Aquí se presentan dos casos: a. Si j ≤ n y rj +1 ≤ m , entonces aún es posible seguir eliminando; incremente j en 1 y vuelva al paso 2.

b. En caso contrario, se ha completado la eliminación Gauss. 7. La matriz reducida final se interpreta como la matriz aumentada para un sistema de ecuaciones y proceda a encontrar las soluciones, si existen, mediante la sustitución en reversa. Cabe señalar que el algoritmo descrito anteriormente es para sistemas lineales de m ecuaciones y n incógnitas donde m no necesariamente es igual a n . Una forma de seleccionar el pivote puede ser utilizando la estrategia de escalar el sistema de ecuaciones nuevamente por una constante diferente de cero y reemplazar cada variable por un múltiplo diferente de cero de una nueva variable. Por lo que en algunos casos es necesario escalar el sistema de ecuaciones original antes de utilizar el proceso de eliminación Gaussiana. Con base a esto existen dos maneras de evitar los pivotes pequeños.

•

Pivote parcial: Las incógnitas se eliminan del modo común x1 , x2 ,..., xm en el cual el pivote para la eliminación con xi es el coeficiente de xi en las ecuaciones enumeradas i , i + 1,..., m que tiene el mayor valor absoluto.

•

Pivote completo: Las incógnitas no se eliminan necesariamente en el orden acostumbrado. La primera variable eliminada es aquella variable xi con el coeficiente mayor (en valor absoluto), y ese coeficiente se utiliza como pivote. La variable r-ésima eliminada es aquella variable con el coeficiente mayor (en valor absoluto) que el de todas aquellas m − r + 1 variables aun sin eliminar, en las ecuaciones m − r + 1 restantes y ese coeficiente se utiliza como pivote.

27

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

En [37] se presentan una comparación de algunas estrategias de pivote que se pueden aplicar a las estrategias tradicionales utilizadas en el proceso de eliminación Gaussiana. Pero en la práctica es suficiente utilizar el pivote parcial [02]. Por otro lado la eliminación Gaussiana presenta el mismo procedimiento explicado en la sección 2.3, excepto que no se realizan intercambios de renglones lo cual genera la existencia de pivotes pequeños en el transcurso de las operaciones elementales con filas y columnas.

2.3.1 Eliminación de Gauss-Jordan El proceso de eliminar renglones (ecuaciones) tanto arriba como debajo de la diagonal principal de una matriz A se le conoce como eliminación de Gauss-Jordan.

2.3.2 Inversa mediante la eliminación de Gauss El cálculo de la inversa de una matriz mediante el uso de la técnica de eliminación Gaussiana se basa en el siguiente teorema. Teorema 2.4 Matrices Inversas

Sea A una matriz de (m, n ) . 1. La matriz R nxm , separada en columnas como

R = [r1 , r2 ,...rm ] , es una inversa

derecha de A si y solo si las ri , resuelven las ecuaciones Ari = ci para 1 ≤ i ≤ m . 2. La matriz L nxm , con transpuesta LT mxn

particionada en columnas como

LT = [ I 1 , I 2 ,..., I n ] , es una inversa izquierda de A si y solo si las I i resuelven las

ecuaciones A T I i = ei para 1 ≤ i ≤ n . 3. Si m = n de modo que A sea una matriz cuadrada, entonces la matriz X mxm es una inversa (bilateral) de A si y solo si X satisface las condiciones para R en el punto 1 y para L en el punto 1. A partir del Teorema 2.4 se define el siguiente teorema: Teorema 2.5 Inversa bilateral.

Sea A una matriz, entonces 1. Si existen para A una inversa izquierda L y una inversa derecha R , entonces son iguales y son una inversa (bilateral). 2. Cualesquiera dos inversas (bilaterales) de A son idénticas.

28

Capítulo 2:Factoriozación LU

2.4 Existencia de soluciones a sistemas de ecuaciones El análisis más simple para un sistema de ecuaciones de la forma Ax = b , es tomando en cuenta la existencia de una incógnita ax = b y una variable x ; por lo que la solución esta dada por x = b / a . Para la cual existen tres posibilidades: 1. Si a ≠ 0 , entonces x = b / a tiene sentido, y ésta es la única solución de esta ecuación. 2. Si a = 0 , existen dos posibilidades: a. Si b ≠ 0 , entonces en la ecuación, se puede encontrar x = 0 | x = b ≠ 0 , y no existe una solución x . Se dice que “no existe solución” o que “la ecuación es inconsistente” ya que implica que 0 = b ≠ 0 , lo cual es contradicción. b. Si b = 0 , entonces hay infinitamente muchas soluciones: Cada número x es una solución, ya que 0 x = b = 0 sin importar el valor de x .

2.5 LU y eliminación de Gauss 2.5.1 LU y eliminación de Gauss sin intercambios Teorema 2.6 Suponga que A es una matriz de (m, n ) . Entonces:

a) La eliminación de Gauss sin intercambios puede llevarse a término con pivotes no cero a modo de producir una forma reducida unitaria triangular superior U , si y solo si las submatrices principales A k son no singulares para 1 ≤ k ≤ m , donde A k es la esquina superior izquierda, (k , k ) de A : A k

ij

= A

ij

para 1 ≤ i , j ≤ k

(en particular A es no singular). b) La eliminación de Gauss sin intercambios puede completarse como en el inciso a), si y solo si A = LU , donde U es una matriz reducida unitaria triangular superior de a) y L es una matriz triangular inferior con los pivotes α ij de la eliminación de Gauss sobre la diagonal principal y los negativos de los multiplicadores mij de la eliminación de Gauss en las subdiagonales; L A = Lo Uo donde:

29

ii

= α ii , L

ji

= −m ji . Equivalente a

(2.3)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

L 0 = LD0 −1 : Es una matriz unitaria triangular inferior y U 0 = D0 U es una

matriz triangular superior con los pivotes α ii sobre la diagonal principal.

D0 = diag(α 11 ,..,α mm ) y α ii = U 0

ii

= L

ii

(2.4)

si la eliminación de Gauss en la matriz A de (m, m) puede completarse satisfactoriamente pero sin intercambios, esto implica que el proceso es equivalente a escribir A como un producto LU . Cuando A es de (m, m) . Con base al Teorema 2.6 L 0 se encuentra a partir de L dividiendo cada columna de L entre el pivote en esa columna, mientras que U 0 se encuentra a partir de U multiplicando cada renglón de U por el pivote de ese renglón.

En programas computacionales, por lo general no se divide entre este pivote, ya que éste se comienza a eliminar. El efecto es que la forma reducida que se obtiene es la U 0 en lugar de U . Los elementos en L 0 por debajo de la diagonal principal son tan solo los negativos de

los multiplicadores utilizados durante esta forma de eliminación. De manera inversa al tener L 0 y U 0 , se podría encontrar L = L o D0 y U = D0 −1 U 0 , donde D0 es la matriz diagonal que

se forma de los elementos de la diagonal principal de U .

2.5.2 LU y eliminación de Gauss con intercambios Con intercambios, toda matriz A no singular puede ser reducida mediante eliminación de Gauss a una matriz triangular superior. En caso de que la matriz A sea singular, el proceso estándar de eliminación de Gauss puede llevarse a cabo pero no producirá una matriz triangular superior. En este caso se hecha mano del Teorema 2.6 para descomponer A = L 0 U 0 , ya que no es posible la factorización A = LU , con base a este teorema con L 0

unitaria triangular inferior, si se permiten elementos ceros sobre la diagonal principal de U . Esto se puede compactar en el siguiente teorema. Teorema 2.7 Suponga que A es de (m, m) . Entonces:

a) A se puede definir mediante la eliminación de Gauss con intercambios de la manera usual, a una matriz unitaria triangular superior U si y solo si A es no 30

Capítulo 2:Factoriozación LU

singular. Cuando se obtiene dicha U , existe una matriz de permutación P y una matriz triangular inferior tal que (2.5)

PA = LU y A = PT LU

Existe una matriz de permutación P , una matriz unitaria triangular inferior L 0 y una matriz triangular superior U 0 tales que

(2.6)

PA = L 0 U 0 y A = P T LU

b) A es no singular si y solo si U 0 es no singular, para A no singular L 0 y U 0 se relacionan con las L y U de a) por (2.7)

U = D0 − 1 U 0 y L = L 0 D 0

donde: D0 = diag(α ii ,...,α mm ) y α ii = U 0

ii

= L

ii

(2.8)

En la siguiente sección se presentan algunos de los métodos más comunes que se utilizan para encontrar la factorización LU.

2.6 Métodos para encontrar la factorización LU Existen dos métodos comunes de realizr una factorización LU: 1. Eliminación Gaussiana. 2. Cálculo directo. a. Forma Doolittle b. Forma Crout c. Forma Cholesky.

El primer método descrito anteriormente crea un matriz triangular inferior durante el proceso de eliminación, la matriz triangular inferior, con 1’s en la diagonal principal, es generada con los multiplicadores utilizados durante la eliminación. Suponiendo el caso en el cual la matriz A sea una matriz tridiagonal, la factorización LU se puede encontrar muy eficientemente a partir de los vectores que contienen los elementos de la diagonal principal, los elementos arriba de la diagonal y los elementos debajo de la diagonal [28]. Donde la matriz tridiagonal tiene la forma 31

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

 a11  a A =  21  ⋮   am 1

⋯ a1 n   ⋯ a2 n  ⋱ ⋮   ⋯ amn 

a12 a22 ⋮ am 2

(2.9)

Por lo tanto la factorización LU de A asumiendo que A ∈ℜ4 x 4  1  bb2 L =  0   0

0

0

1 bb3

0 1

0

bb4

0  0 0  1 

y

 dd1  0 U =  0   0

a1

0

dd2 0

a2 dd3

0

o

0   0  a3   dd4 

(2.10)

El segundo método es otra forma de obtener una factorización LU, el cual se pueden clasificar en tres formas más comunes, con base a suposiciones hechas acerca de los elementos contenidos en la diagonal de L y de U . Estos métodos son [19]: 1. La forma Doolittle, ésta forma elementos en la diagonal principal de la matriz triangular inferior iguales a la unidad. Lo cual proporciona los mismos factores producidos por la eliminación Gaussiana. Para llevar a cabo esta forma de factorizar una matriz A , partimos de que LU = A , donde

1 0   l21 1  l31 l32   l41 l42

0 0 1 l43

0   u11 u12  0   0 u22 0 0 0  1 0 0

Y el procedimiento para obtener

u13 u23 u33 0

u14   a11   u24   a21 = u34   a31   u44   a41

a12

a13

a22

a23

a32 a42

a33 a43

a14   a24  a34   a44 

es el siguiente: En el primer

(2.11)

paso obtenemos

inmediatamente u11 = a11 y entonces se resuelven los elementos restantes mediante el procedimiento de sustitución hacia delante.

2. La forma Crout, este método factoriza una matriz A , en dos matrices LU triangulares inferior y superior respectivamente. Pero la característica principal de este método de factorización es que la matriz triangular superior contiene elementos igual a la unidad en la diagonal principal [13]. Siguiendo con la suposición de que A ∈ℜ4 x 4 para fines explicativos, tenemos que LU = A puede ser expresado por

32

Capítulo 2:Factoriozación LU

 l11 0   l21 l22  l31 l32   l41 l42

0 0 l33 l43

0   1 u12  0 0 1 0 0 0  l44   0 0

u13 u23 1 0

u14   a11   u24   a21 = u34   a31   1   a41

a12

a13

a22

a23

a32 a42

a33 a43

a14   a24  a34   a44 

(2.12)

El procedimiento para obtener los elementos de las matrices L y U es la multiplicación sucesiva de los renglones de L con cada una de las columnas de la matriz U , de tal manera que de forma progresiva (sustitución hacia delante) obtener cada uno de los factores.

3. La forma Cholesky, este método de factorizar una matriz A requiere que los elementos correspondientes a la diagonal de L y U sean iguales, lo cual es muy útil cuando el algoritmo se implementa para matrices definidas positivas simétricas, ya que en tal caso la matriz preserva la simetría y produce una factorización L = U T . 0 0   α11 0   α l 0 0 22 , L =  21  l31 l32 α 33 0     l41 l42 l43 α 44 

 α11 u12 u13 u14    0 α 22 u23 u24  U=  0 0 α 33 u34    0 0 α 44   0

(2.13)

Para una matriz simétrica definida positiva, una modificación a la forma Cholesky proporciona una factorización simétrica de la forma [01] : LDLT

(2.14)

2.7 LU para resolver el problema de mínimos cuadrados En el modelo (1.5) de mínimos cuadrados, el cálculo de A T A y la solución de este sistema de ecuaciones a menudo produce errores cuando los cálculos se hacen en computadoras con aritmética de punto flotante. La factorización LU es una alternativa para resolver este problema, considérese las ecuaciones A T Ax = A T y

(2.15)

A T A = LU

(2.16)

cuando

33

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

tenemos que (LU )x = A T y

(2.17)

por lo tanto la solución para A T A en (2.15) mediante LU es ˆ = U −1 L−1 A T y x

(2.18)

2.8 Función LU en MATLAB ® MATLAB® cuenta con la función lu que se incluye en el paquete LAPACK. Esta función expresa una matriz A como el producto de dos matrices triangulares, una de ellas es una permutación de una matriz triangular inferior y la otra una matriz triangulas superior. En algunas ocasiones a la factorización LU se le llama factorización LR; la matriz A puede ser cuadrada o rectangular. La sintaxis para obtener la factorización LU, la cual se obtuvo de [34], es la siguiente:

[L,U] = lu(A), regresa una matriz triangular superior en U y una matriz triangular inferior L , la cual es un producto de una matriz triangular inferior y una matriz de permutación P tal que A = LU .

[L,U,P] = lu(A), regresa una matriz triangular superior en U , una matriz triangular inferior con diagonal unitaria en L y una matriz de permutación en P , tal que LU = PA .

Y = lu(A), si A es “sparse”, la function lu regresa una matriz triangular inferior L estricta, por ejemplo, sin que la diagonal de L sea unitaria; y la matriz triangular superior U incluida en la misma matriz Y , tal que YU + L − speye( size( A)) . La matriz de permutación no existe.

[L,U,P,Q] = lu(A), para una matriz “sparse” no vacía A , la función lu regresa una matriz triangular inferior L , una matriz triangular superior U , una matriz de permutación de filas P y una matriz Q de reordenamiento de columnas; tal que PAQ = LU . Esta sintaxis

ocupa menor tiempo y mejor uso de la memoria en el proceso de cálculo con respecto a las anteriores. Si A es una matriz vacía y no “sparse”, lu despliega un mensaje de error.

[L,U,P] = lu(A,thresh), con el uso de esta sintaxis se controla el pivote en matrices “sparse”, donde tres es un pivote que se mantiene en el intervalo [0.0,1.0]. EL pivote ocurre cuando la diagonal completa en una columna tiene magnitud menor que trheshveces la magnitud de cualquier subdiagonal completa en esa columna. Si thresh = 0.0 fuerza a que el pivote sea el elemento que se encuentra en la diagonal. EL valor por defecto es thresh = 1.0 .

34

Capítulo 2:Factoriozación LU

En la siguiente sección se presenta la etapa de evaluación del algoritmo, así como los diferentes métodos de obtener la factorización LU. El valor del pivote es el que tiene por defecto.

2.9 Evaluación del algoritmo LU Para evaluar la programación del algoritmo LU se analizan dos casos, el primero de ellos consta en realizar una factorización y comprobar que dicha factorización proporciona la matriz original, y el segundo caso es aplicar este algoritmo para identificar los parámetros de una sistema propuesto del cual se conocen sus parámetros. Sea la matriz  1 2 3   A =  4 5 6  7 8 0  

(2.19)

Obteniendo su factorización LU con la función lu de MATLAB®, tenemos que 0 0  1   L =  0.1429 1 0  0.5714 0.5000 1   

Además de una matriz

P

y

8 0  7   U =  0 0.8571 3  0 0 4.5  

(2.20)

de permutación 0 0 1   P = 1 0 0  0 1 0  

Por lo tanto para comprobar que la factorización obtenida es correcta probamos 7 8 0 7 8 0     P * A =  1 2 3,L * U =  1 2 3  4 5 6      4 5 6

(2.21) P*A = L*U

(2.22)

Con lo cual concluimos que la factorización está bien hecha. Esta factorización se realiza siguiendo el procedimiento de la eliminación Gaussiana, que es el método convencional, para obtener la matriz L , y la matriz U se forma con los multiplicadores utilizados como pivote. Ahora evaluamos las diferentes formas (Doolittle, Crout y cholesky) de obtener la factorización LU para la misma matriz (2.19). Para el primer método Doolittle, que se programó en código C, las matrices L y U son

35

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.  1 0 0 1 2 3      L =  4 1 0  , U =  0 −3 −6  7 2 1  0 0 −9     

(2.23)

 1 2 3   L * U =  4 5 6  7 8 0  

(2.24)

Comprobando que A = LU ,

Por lo tanto se puede decir que la factorización es correcta ya que la forma tanto de la matriz L como de la matriz U cumple con la propiedad de tal factorización.

Para obtener la forma Crout de A también se codificó en lenguaje C el procedimiento que se tiene que realizar. Las matrices L y U resultantes son 1 0 0    L =  4 −3 0   7 −6 −9   

y

1 2 3   U = 0 1 2 0 0 1  

(2.25)

Ahora comprobando que A = LU  1 2 3   A = LU =  4 5 6   7 8 0  

(2.26)

De acuerdo con la propiedad de L y U para factorizar una matriz A cualquiera se cumplen, la factorización es correcta. Es importante señalar que MATLAB® no cuenta con funciones para llevar acabo las dos factorizaciones anteriores, ya que solo permite obtener la factorización LU mediante eliminación Gaussiana y de la forma cholesky, que es el último método que se analizó para poder obtener una factorización triangular. La forma cholesky de la matriz (2.19) no se puede obtener debido a que no cumple con la propiedad de ser simétrica definida positiva. Por o que se propone una matriz A *  4 8 12    A* =  8 20 20  12 20 41   

(2.27)

utilizando la función chol de MATLAB® se obtiene sus respectivas matrices triangular inferior y triangular superior  2 0 0   L =  4 2 0  6 −2 1   

y

2 4 6    U = L ' =  0 2 −2  0 0 1   

36

(2.28)

Capítulo 2:Factoriozación LU

Comprobando que A* = LL '  4 8 12    LL ' =  8 20 20  12 20 41   

(2.29)

También es posible obtener la matriz A = U H DU , donde D = diag (L) . Con esto se puede concluir que la factorización es correcta.

2.10 Conclusiones En este capítulo se analizaron los conceptos mínimos necesarios para poder entender las diferentes formas de factorizar una matriz A en una matriz triangular inferior L y una matriz triangular superior U . Dentro de las formas más comunes de encontrar dicha factorización existen dos: •

La eliminación Gaussiana

•

Los métodos directos que básicamente son tres: a) La forma Doolittle b) La forma Crout c) La forma Cholesky

Las técnicas que se encuentran dentro del contexto de métodos directos se aplican con base a suposiciones hechas acerca de los elementos que se encuentran en la diagonal tanto de la matriz L y la matriz U . Se realizó una etapa de evaluación numérica a cada una de formas de obtener la factorización LU con la finalidad de verificar que la programación de los algoritmos es correcta, teniendo éxito todos los casos. Estas formas de factorización aplicadas a la identificación de sistemas permitirán tener una reducción en el tiempo de cálculo numérico, tal como se ha reportado en la literatura.

37

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

38

Capítulo 3 Factorización QR 3.1 Introducción En este capítulo se analiza la factorización QR, así como los métodos más comunes (GramSchmidt, Rotaciones de Givens y Reflexiones Househoder) que se utilizan para encontrar dicha factorización. El capítulo se organiza de la siguiente manera: En la Sección 3.2 se presentan los conceptos y definiciones generales de la factorización QR. En la Sección 3.3 se presenta la teoría básica respecto a las bases ortogonales y ortonormales, debido a que la factorización QR tiene como característica principal que la matriz Q es ortogonal. En la Sección 3.4 se presenta el método más común para generar bases ortogonales partiendo de de un conjunto generador. En la Sección 3.5 se aborda el tópico de proyecciones ortogonales, éste concepto es útil cuando la factorización QR se obtiene mediante Rotaciones de Givens o Reflexiones de Householder. En la Sección 3.6 se describe el algoritmo para obtener la factorización QR, además de los métodos más comunes para encontrar esta factorización. En la Sección 3.7 se presenta la aplicación de la factorización QR al problema de mínimos cuadrados. En la Sección 3.8 se muestra la función QR que contiene MATLAB® así como la sintaxis para obtener tal factorización. En la Sección 3.9 se realiza la etapa evaluativo del algoritmo con sus diferentes métodos de obtención. Y para finalizar en la Sección 3.10 se exponen las conclusiones del capítulo.

39

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

3.2 Definiciones y conceptos generales Esta sección describe las definiciones y conceptos generales utilizados durante el análisis de la factorización QR y los diferentes métodos para encontrar dicha factorización, extraídos de [02]. Definición 3.1 Espacio columna y espacio renglón.

El espacio columna real (o complejo) de A es el subespacio de ℝ m (o de ℂm ) que está generado por el conjunto de columnas de A . El espacio renglón real (o complejo) de A es el subespacio vectorial real (o complejo) de todas las matrices (1, n ) que está generado por el conjunto de renglones de A .

Definición 3.2 Rango de una matriz.

El número de renglones diferentes de cero en una matriz A obtenida mediante operaciones elementales de renglón en A se llama rango de A .

Definición 3.3 Matriz Transpuesta.

Sea A una matriz de (n, m) . La transpuesta A T de A que se obtiene intercambiando los renglones y columnas de A , el primer renglón se vuelve la primera columna, y así sucesivamente. Esto es, AT

ij

= A

ji

para 1 ≤ i ≤ n y 1 ≤ j ≤ m

(3.1)

Definición 3.4 Matriz transpuesta Hermitiana.

La transpuesta hermitiana A H de A es la matriz (n, m) que se obtiene tomando el complejo conjugado de sus elementos en A T , la primera columna de A H consiste en los complejos del primer reglón de A , y así sucesivamente. Esto es, AH

ij

= A

ji

para 1 ≤ i ≤ n y 1 ≤ j ≤ m

(3.2)

Teorema 3.1 Espacios renglón/columna y rango.

Si A tiene rango k , entonces a) El espacio columna de A tiene dimensión k , y una base de él consiste en el conjunto de columnas de A correspondientes a las columnas delanteras de cualquier forma reducida de Gauss de A .

40

Capítulo 3:Factoriozación QR

b) El espacio renglón de A tiene dimensión k , y una base de él consiste en los renglones no cero en cualquier forma reducida de Gauss de A . c) La matriz A de (m, m) , es no singular si y solo si sus columnas forman un conjunto linealmente independiente, en cuyo caso el espacio columna tiene dimensión p . d) La matriz A (m, m) , es no singular si y solo si sus renglones forman un conjunto linealmente independiente, en cuyo caso el espacio renglón tiene dimensión p .

Corolario 3.1 Rango de renglón = rango de columna.

Los rangos de A , A T y A H son iguales.

Teorema 3.2 Leyes para las transpuestas. ( A T )T = A y ( A H )H = A ( A ± B)T = A T ± BT y ( A ± B)H = A H ± B H

(3.3)

( AB)T = BT A T y ( AB)H = B H A H

Definición 3.5 Matriz simétrica.

Una matriz A (cuadrada) para la cual A T = A , de modo que A

ij

= A

ij

i y j , se dice

que es simétrica. Una matriz B (cuadrada) para la cual BT = B , de modo que B

ij

= B

ij

para toda i y j se dice que es hermitiana.

Definición 3.6 Producto interior.

Para el caso real Sea V un espacio vectorial real. Un producto interno en V es una función que asigna a cada par ordenado de vectores u y v en V , un número real que se denota (u, v ) y que satisface las siguientes propiedades: 1) (u, v ) = ( v , u) para toda u y v en V . 2) (α u + β v + w ) = α ( u, w ) + β ( v , w ) para toda u, v y w en V y todos los números reales α y β . 3) (u, u) > 0 si u ≠ 0 , y (u, u) si y solo si u = 0 .

41

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

El ángulo θ entro dos vectores diferentes de cero u y v se define por su coseno: cos θ =

( u, v ) ( u, u )1 /2 ( v , v )1 /2

(3.4)

Para el caso complejo Sea V un espacio vectorial complejo. Un producto interno en V es una función que asigna a cada par ordenado de vectores u y v en V , un número posiblemente complejo que se denota (u, v ) y que satisface las siguientes condiciones: 1. (u, v ) = ( v , u) para toda u y v en V . 2. (α u + β v + w ) = α (u, w ) + β ( v , w ) y (α u + β v + w ) = α ( u, w ) + β ( v , w ) para toda u , v y w en V y todos los números complejos α y β .

3. (u, v ) > 0 si u ≠ 0 , y (u, v ) = 0 si y solo si u = 0 .

Teorema 3.3 Normas de productos interiores.

Sea (u, v ) un producto interior en V , y defínase a v = ( u, v )1 /2 . Entonces

.

es una

norma en V (se dice que es la norma inducida por el producto interior).

Teorema 3.4 Desigualdad de Cauchy-Schwarz.

Sea (u, v ) un producto interior en el espacio vectorial real o complejo V . Entonces (u, v ) ≤ (u, u )1/ 2 ( v , v )1 /2 para toda u , v en V .

(3.5)

En el caso matricial, Sean x y y matrices columna mx1 . Entonces xH y ≤ x

Definición 3.7 Ortogonalidad.

..

Sea ( , ) un producto interior en V y sea

2

y

2

(3.6)

. su norma inducida según el teorema 2.

a) Se dice que dos vectores u , v son ortogonales si y solo si (u, v ) = 0 . b) Se dice que un conjunto de vectores es ortogonal si y solo si cada par de vectores del conjunto es ortogonal (u, v ) = 0 para toda u ≠ v de ese conjunto.

42

Capítulo 3:Factoriozación QR

c) Si se usa un vector no cero u para producir v = u / u tal que v = 1 , entonces se dice que u ha sido normalizada para producir el vector normalizado v . d) Se dice que un conjunto de vectores es ortonormal si y solo si el conjunto es ortogonal y v = 1 para todo v en el conjunto. Teorema 3.5 Proyecciones ortogonales.

Sea V un espacio vectorial con un producto interior. Sea Vo el subespacio de V generado por un conjunto ortogonal (3.7)

S = {v1 , v2 , ⋯ vn }

de vectores diferentes de cero. Defina la proyección ortogonal Po sobre Vo como sigue: para cualquier v en V , sea Po v = α1 v1 + ... + α n v n donde α1 =

( vi , v) vi , vi

(3.8)

Entonces: a) v − Po v es ortogonal a todo vector v o en Vo . b) Po ( u + v ) = Po u + Po v para todo u y v en V . c) Po (α v ) = α Po v para todo escalar α y todo vector v en V .

Teorema 3.6 Mejor aproximación.

Sea V un espacio vectorial con un producto interior y con una norma inducida

. , y sea V

o

el subespacio generado por el conjunto ortogonal de vectores diferentes de cero v 1 ,..., v n . Entonces, para cualquier v , Po v es el único punto más cercano en Vo a v , y v − Po v es la distancia de v a Vo , en el sentido de que Po v en Vo y v − Po v < v − v o para todo v o ≠ Po v en Vo

(3.9)

3.3 Bases ortogonales y ortonormales Todo conjunto ortogonal de vectores diferentes de cero es linealmente independiente. Suponga que

43

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

S = { v1 , v 2 , ⋯ v n }

(3.10)

c1 v1 + ... + cn v n = 0

(3.11)

Es ese conjunto y que

Tomando le producto interior con cada v i resulta

0 = ( v j , 0) = ( v j , c1 v1 + ... + cn v n ) = c j ( v j , v j )

(3.12)

Lo que significa que c j = 0 porque ( v j , v j ) > 0 . Por lo tanto, cualquier conjunto generador ortogonal de vectores diferentes de cero es un conjunto generador ortogonal linealmente independiente, esto es, una base ortogonal para Vo . Se puede escribir, como

α1v1 + ... + α n v n

(3.13)

Si sucediera que v mismo está en Vo , entonces, verdaderamente v es la mejor aproximación a sí misma en Vo , y se obtiene

v = α1 v1 + ... + α n v n

(3.14)

Es posible expresar un vector v como una combinación lineal de vectores en una base ortogonal sin tener que resolver ecuaciones para determinar los coeficientes; simplemente, se evalúan algunos productos interiores para obtener directamente los coeficientes.

Teorema 3.7 Bases ortogonales. Sea B = {v 1 , v 2 ,..., v n } una base ortogonal (u ortonormal). Entonces, la representación de cualquier vector v con respecto a la base ortogonal B puede escribirse inmediatamente:

v = α1 v1 + ... + α n v n donde α i =

( vi , v) ( vi , vi )

(3.15)

3.4 Creación de bases ortogonales y ortonormales: Gram-Schmidt. El proceso de Gram-Schmidt produce una base ortogonal partiendo de un conjunto generador y detecta si el conjunto generador es linealmente dependiente.

Teorema 3.8 Gram-Schmidt por proyecciones.

Sea S = {v 1 , v 2 ,..., v n } un conjunto generador del espacio vectorial V que tiene un producto interior. El proceso de Gram-Schmidt es como sigue: Sea Vi el subespacio de V generador

44

Capítulo 3:Factoriozación QR

por

{v 1 ,..., v i } y sea

Pi la proyección ortogonal sobre Vi (calculada

como indica el

Teorema 3.5 mediante los conjuntos ortogonales siguientes, Bi ; si Vi = 0 sea Pi v = 0 para toda v . 1. Defina u1 = v 1 . 2. Para 2 ≤ i ≤ n , defina u i = v i − Pi −1 v i , Lo que sigue es válido para los vectores producidos mediante este procedimiento: 1. B = {u1 , u2 ,..., u n } es un conjunto ortogonal que genera a V . 2. Para 1 ≤ i ≤ n , Bi = {u1 ,..., u i } . 3. u i = 0 si y solo si v i es linealmente dependiente de los vectores {v 1 ,..., v i } . 4. Puede obtener una base ortogonal para V a partir del conjunto generador ortogonal B omitiendo las columnas iguales a cero, si las hubiere. 5. Si S es una base para V , entonces B es una base ortogonal para V .

Procedimiento clásico (tradicional) de Gram-Schmidt. Dado un conjunto generador

{v 1 ,..., v n } : a. Defina u1 = v 1 . b. Para 2 ≤ i ≤ n defina ui = v i − α1u1 − ... − α i −1u i −1

(3.16)

en donde α ji = (u j , v i )/( u j , u j ) si u j ≠ 0 y α ji = 0 si u j = 0 .

3.5 Proyecciones ortogonales A continuación se reformularán las proyecciones ortogonales en ℝ n y ℂ n mediante terminología matricial. Suponga que Vo es un subespacio de ℝ n generado por el conjunto ortogonal de vectores S = { v 1 , v 2 ,..., v n } y que v es otra matriz mx1 ; se desea calcular la proyección ortogonal Po v . Se supondrá que S es un conjunto ortonormal; es decir que v i

condiciones de ortogonalidad 0 = ( v i , v j ) = v i T v j para i ≠ j .

45

2

= 1 , además de las

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Defínase la matriz Q de (m, n ) , por Q = [ v1 , v 2 ,..., v n ]

(3.17)

QT Q = [ v1 , v 2 ,..., v n ]T [ v1 , v 2 ,..., v n ]

(3.18)

si calcula QT Q , encontramos que

tiene como elemento (i , j ) , justamente a v i T v j ; que es igual a uno si y solo si uno si i = j , y es igual a 0 si i ≠ j . Esto es, QT Q = I n ; ésta es la reformulación matricial del hecho de que S es ortonormal. De acuerdo al Teorema 3.5. Po v = α1 v1 + ... + α n v n

(3.19)

Para α i adecuados en notación matricial, esto significa que Po v = Qα , donde Q está definida en la expresión (3.17) y

α = [α1 , α 2 , ⋯ α n ]T

(3.20)

Nuevamente de acuerdo al Teorema 3.5 α i se calcula como

αi =

( vi , v) = v iT v ( vi , vi )

(3.21)

En notación matricial, esto dice que α = QT v . Por lo tanto se obtiene Po v = QQT v

(3.22)

En general QQT ≠ I m ; QT es una inversa izquierda de Q porque QQT = I m , pero generalmente no es una inversa derecha porque no está suponiendo que Q sea cuadrada. Esta fórmula para Po demuestra que la proyección ortogonal no es más que da multiplicación por una matriz especial QQT , llamada matriz de proyección.

Teorema 3.9 Matrices de proyección. Sea Q una matriz de (m, n) que tiene columnas ortonormales en ℝ m (o en ℂ m ), esto es QT Q = I n o (Q H Q = Tn )

(3.23)

y sea Vo es subespacio generado por la base ortonormal para Vo formada por la columna de Q . Entonces:

46

Capítulo 3:Factoriozación QR

a) La proyección ortogonal Po en Vo , como se describió en los

Teorema 3.5 y

Teorema 3.6 se calcula con la matriz de proyección QQT ó (QQT ) como Po v = Po v

(3.24)

en donde Po es la matriz de (m, m ) , Po = QQT o bien (Po = QQT ) . b) La matriz de proyección Po satisface lo siguiente: a. Po es simétrica (o hermitiana). b. Po 2 = Po c. Po (I m − Po ) = (I m − Po )Po = 0 d. (I m − Po )Q = 0

3.6 Descomposición QR La factorización QR es una alternativa más de solución para resolver sistemas lineales de la forma Ax = b . Primero mencionaremos el teorema de base.

Teorema 3.10 Descomposición QR. Sea A una matriz de (m, n ) y rango k . Entonces: 1. A puede escribirse en su descomposición QR no normalizada como A = Qo R o , en donde: a. Qo es de (m, n) y tiene columnas ortogonales (de las cuales k son diferentes de cero y n − k son cero) que generan el espacio columna de A. b. R o es (m, n ) , triangular superior unitaria y no singular. c. La norma 2 de la i-ésima columna de Q o es igual a la distancia de la i-ésima columna de A al espacio generado por las primeras i − 1 columnas de A . 2. A puede escribirse en su descomposición QR normalizada como A = QR , en donde: a. Q es de (m, k ) y tiene columnas ortonormales que generan el espacio columna de A . b. R es triangular superior de (k , n) y tiene rango k .

47

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

c. Si k = n , entonces R

ii

es igual a la distancia desde la i-ésima columna de

A hasta el espacio generado por las primeras i − 1 columnas de A .

Con base al Teorema 3.8, se puede calcular un conjunto generador ortogonal partiendo de un conjunto generador general. Suponga que v 1 , v 2 ,..., v n son matrices columna (m, n ) , y considere la implementación tradicional del procedimiento de Gram-Schmidt. A partir de la definición de matrices columna u j en (3.16), se tiene v j = v1α1 j + v 2α 2 j + ... + v j −1α j −i , j + u

(3.25)

Representando en notación de matrices separadas, es tan solo

[ v1

⋯ v n ] = [ v1

1 α12 0 1 ⋯ vn ]  ⋮ ⋮  0 0

⋯ α1q  ⋯ α 2 q  ⋱ ⋮   ⋯ 1 

(3.26)

Sea Q o la matriz (m, n) , [ v 1 ,..., v n ] de la representación matricial en (3.25) y sea R o la matriz unitaria triangular superior (n, n ) allí mismo; como A = [ v 1 ,..., v n ] , es posible escribir de la misma forma como A = Qo R o

(3.27)

Las ui de Qo se construye en (Gram-Schmidt por proyecciones), y por lo tanto son mutuamente ortogonales; es decir, que Q o tiene columnas ortogonales, algunas de las cuales podrán ser iguales a cero. Sean Q y R o las matrices obtenidas para la eliminación de columnas de ceros de Q o y los renglones correspondientes de R o , y mediante la división de cada columna no cero de Q o entre su norma 2, y multiplicando cada renglón correspondiente de R o por la misma norma 2. Por lo tanto A = QR

(3.28)

con R como triangular superior y Q con columnas ortonormales. La forma (3.26) se llama descomposición no normalizada QR de A , mientras que (3.27) es la descomposición QR normalizada.

48

Capítulo 3:Factoriozación QR

Si se supone una matriz A de (m, n) , donde m = 4 y n = 2 la factorización QR correspondiente tiene la forma A • •  •  •

= •  • •  • = •  •   •  •

Q • • •  • • • •  0 • • • 0  • • • 0

R • •  0  0

(3.29)

donde • indica un elemento con valor diferente de cero. Como Q es ortogonal, por lo tanto QT Q = I y, de aquí que, Q −1 = QT . La matriz R es triangular superior; y por lo general los últimos elementos m − n de R son elementos igual a cero [43]. Llamemos q j , r j , y a j a las columnas de Q , R , A respectivamente, y el producto de Q y R lo escribimos

[q1 | q 2 | ⋯ | q m ] = ...r j ... = ...a j ...

(3.30)

desde el punto de vista matricial, la tercera columna de A esta definida como una combinación lineal de las primeras tres columnas de Q . Es decir, r1,3q1 + r2,3q 2 + r3,3q 3 = a3

(3.31)

Ya que ri , j = 0 para i > j , solamente son necesarias las columnas de la 1 a la j de la matriz Q para formar la columna j de la matriz A . Por lo tanto la matriz A se puede representar

a partir de las primeras n columnas de la matriz Q . De tal manera se asume que la matriz A tiene rango n , es decir, que todas las n columnas de A son linealmente independientes. Si el rango de la matriz A es igual a k (

rango( A ) = k ) y k < n , entonces solo son necesarias las primeras k columnas de la matriz Q para formar la matriz A [14]. Las últimas m − n renglones de R deben de ser cero, debido a que las columnas de Q son ortogonales, y por lo tanto linealmente independientes. Como la matriz Q es de (m, m) y las columnas de Q son ortogonales, forman una base para el espacio columna de A ( rango( A ) ). Las últimas m − n columnas de Q forman una base para el subespacio de ℝ m que yace fuera del rango de la matriz A . A este espacio se le conoce como complemento ortogonal del rango de A y se define como rango( A ) ⊥ , (┴

49

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

ɶ es la matriz que se compone de las n columnas de Q , y Q complemento ortogonal). Si Q c es una matriz que se compone de las columnas restantes de Q , donde Q c se refiere al hecho de que las columnas forman una base de rango( A) ⊥ , podemos escribir la matriz Q en forma matricial de la siguiente manera ɶ Q  Q = Q c

(3.32)

Por lo tanto la representación matricial de R es ɶ R R=  0

(3.33)

ɶ es matriz consistente de (n, n) correspondiente a la parte triangular superior de la donde R matriz R con elementos diferentes de cero y 0 se refiere a la matriz de ceros de dimensión ((m − n), n) . Con el conocimiento de que solamente son necesarios las primeras n columnas de la matriz

Q para crear la matriz A usando los coeficientes en R , se puede reducir el esfuerzo computacional y el uso de memoria durante el proceso para realizar la factorización QR. La factorización QR de la matriz A conocida como “economy-size” se define como

ɶ ɶ A = QR

(3.34)

La diferencia entre la factorización QR completa (full QR factorization) y la factorización QR “economy-size” radica

en que la factorización QR completa contiene las m − n

columnas adicionales de Q que forman la base para rango( A) ⊥ . Algunos de los métodos para obtener la factorización QR son: 1. Gram-Schmidt. 2. Reflexiones Householder. 3. Rotación de Givens. El proceso de Gram-Schmidt no se emplea en la práctica para hallar bases ortonormales porque existen otras técnicas mas eficientes para efectuar la ortogonalización mediante computadora .Uno de estos procedimientos es mediante las reflexiones de Householder. Éste método ya se ha incluido en MATLAB®.

50

Capítulo 3:Factoriozación QR

3.6.1 Gram-Schmidt El método de Gram-Schmidt es uno de los métodos de ortogonalización más antiguos que se utiliza para obtener la factorización QR. Éste método se deriva partiendo de la siguiente expresión. (3.35)

A = QR donde

Q ∈ ℝ mxn y R ∈ ℝ nxn . Si a j y q j definen las j-ésimas columnas de A y Q

respectivamente, tenemos que j

a j = ∑ rkj q k

(3.36)

k =1

Premultiplicando por qiT y tomando en cuenta que Q tiene columnas ortonormales, es decir

qiT a j = rij , i = 1: j − 1 . Por consiguiente, q j = q ' j / r jj

(3.37)

donde j −1

a' j = a j − ∑ rkj q k , r jj = q ' j k =1

2

(3.38)

por lo tanto se puede obtener una columna de Q y una columna de R al mismo tiempo. Para asegurar que r jj > 0 se requiere que la matriz A tenga rango completo [14].

3.6.2 Rotaciones de Givens

Otra forma de obtener la factorización QR es mediante las rotaciones de Givens. Las rotaciones de Givens se utilizan para eliminar los elementos de debajo de la diagonal de forma sistemática. Dada una matriz A o = A ∈ ℝ ( m,n ) , se busca una secuencia de matrices ortogonales Ω1 , Ω 2 ,..., Ω k ∈ ℝ ( m ,n ) tales que la matriz A i = Ωi Ai −1 tenga más elementos nulos debajo

de la diagonal de A i −1 , para i = 1, 2,..., k , de forma que la matriz R = A k final es triangular superior; por lo tanto, se tiene: Ω k , Ω k −1 ,..., Ω1A = R

51

(3.39)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

y en consecuencia, A = QR , siendo Q = Ω1T , Ω2T ,..., Ωk T ortogonal y R triangular superior. Se denomina matriz de Givens,  c s Ω[ m,n ] =    −s c 

(3.40)

a cualquier matriz ortogonal que coincide con la unidad, excepto en cuatro elementos, que se construyen de la siguiente forma: c = Ωmm[ m,n ] = Ωnn[ m,n ] = cos θ s = Ωmn[ m ,n ] = senθ

(3.41)

−s = Ω nm[ m ,n ] = − senθ

para algún θ ∈ [−π , π ] . Algebraicamente,  xk , k ≠ i, j  yk =  cxi + sx j , k =i  − sx + cx , k = j j  i

(3.42)

y y j = 0 si

s=

xj xi 2 + x j 2

, c=

xi xi 2 + x j 2

Figura 3.1 Rotación de Givens.

52

.

(3.43)

Capítulo 3:Factoriozación QR

3.6.3 Reflexiones de Householder El método de Householder es más eficiente que el método utilizado por Givens. Los elementos apropiados de una columna entera se eliminan mediante una reflexión (o matriz Householder). La matriz de transformación que describe la reflexión tiene la forma H =I−

2 vvT , 0 ≠ v ∈ ℝ n vT v

(3.44)

La matriz de transformación H es ortogonal ya que HT H = I

(3.45)

H = HT = H −1

(3.46)

Y también porque H es simétrica

Tomando en cuenta las propiedades de ortogonalidad y simetría de la matriz H , la aplicación de H a un vector cualquiera x produce

 2 vT x  Hx = x −  T  v  v v 

(3.47)

La Figura 3.2 muestra porqué H es algunas veces llamado reflector Householder: Éste refleja el vector x alrededor del hiperplano con span( v ) ⊥ [2].

Figura 3.2 Matriz Householder H , veces el vector

x.

Las matrices Householder son una herramienta muy poderosa para insertar ceros en los vectores. ¿Es posible encontrar una matriz Householder H tal que Hx = y , dado los vectores x y y ?. Si H es ortogonal, requerimos que x 2 = y 2 . Por lo tanto  vT x  Hx = y ⇔ x − 2  T  v v v

53

(3.48)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

la ecuación anterior se puede escribir de la siguiente manera α v = x − y para ciertos valores de α . Pero la matriz H es independiente del escalamiento que tenga el vector v así que podemos establecer a α = 1 . Si α = 1 entonces v = x − y y resulta que v T v = xT x + y T y − 2 xT y ,

(3.49)

1 T v v, 2

(3.50)

ya que xT x = y T y v T x = xT x − y T x = por consiguiente Hx = x − v = y Se concluye que proporcionando x 2 = y

2

(3.51)

es posible encontrar una matriz Householder

H tal que Hx = y . Excluyendo el caso en el que x = y , lo cual implicaría que v = 0 , por lo

que la matriz Householder será indefinida. Normalmente se escoge y tal que tenga un patrón especial de ceros. La elección usual es y = σ e1 donde σ = ± x 2 , la cual produce un número máximo de ceros en y . Entonces v = x − y = x − σ e1

(3.52)

En la ecuación anterior existe la libertad de elegir el signo de σ , por lo que se escoge de la siguiente manera para evitar la cancelación en la expresión para v1 . sin g (σ ) = − sin( x1 )

(3.53)

3.7 Aplicación de la Factorización QR a mínimos cuadrados El cálculo de AT A y la solución de este sistema de ecuaciones a menudo produce errores inaceptables cuando los cálculos se hacen en computadoras con aritmética de punto flotante. La descomposición QR normalizada, que puede calcularse con exactitud y eficiencia, proporciona un camino para eliminar la dificultad. Suponga que A = QR es una descomposición QR normalizada de A , de tal modo que, si A es de (m, n) tiene rango k ; Q es de (m, k ) y tiene columnas ortonormales, mientras que R es de (k , n) , triangular superior y tiene rango k .

54

Capítulo 3:Factoriozación QR

Considérense las ecuaciones AT Ax = AT y cuando A = QR

(3.54)

(QR )T (QR )x = (QR )T y

(3.55)

RT QT QRx = RT QT y

(3.56)

al sustituir A resulta

esto es

como Q tiene columnas ortonormales, así que QT Q = I k , por lo tanto, las ecuaciones vienen a ser RT Rx = RT QT y

(3.57)

Teorema 3.11 Inversas izquierdas y rango.

Sea A una matriz de (m, n) y sea k el rango de A . Entonces, A tiene una inversa izquierda L si solo si k = n y n ≤ m . Teorema 3.12 No singularidad y rango.

Sea A una matriz cuadrada de (m, m) . A es no singular si y sólo si el rango de A es igual m.

Teorema 3.13 Espacio renglón/ columna y rango.

Si A tiene rango k , entonces a) El espacio columna de A tiene dimensión k , y una base de él consiste en el conjunto de columnas de A correspondientes a las columnas delanteras de cualquier forma reducida de Gauss de A . b) El espacio renglón de A tiene dimensión k , y una base de él consiste en los renglones no cero en cualquier forma reducida de Gauss de A . c) La matriz A , de (m, m) , es no singular si y sólo si sus columnas forman un conjunto linealmente independiente, en cuyo caso el espacio columna tiene dimensión m . d) La matriz A , de (m, m) , es no singular si y sólo si sus renglones forman un conjunto linealmente independiente, en cuyo caso el espacio renglón tiene dimensión m .

55

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Corolario 3.2 Rango de renglón igual a rango de columna.

Los rangos de A , AT y A H son iguales. De acuerdo al Corolario 3.2 los rangos de R y de R T son iguales; por lo tanto, la matriz RT (n, k ) , tiene rango k . Por el Teorema 3.11, RT tiene una inversa izquierda L tal que LRT = I k . Multiplicando ambos miembros de la ecuación (3.57) por L y usando LRT = I k , Rx = QT y

(3.58)

Si RT Rx = RT QT y ; es la inversa, es claro que si Rx = QT y , entonces RT Rx = RT QT y . Por lo tanto, x resuelve el problema de mínimos cuadrados si y solo si Rx = QT y . Ya que R es de (k , n) , es triangular superior y tiene rango k , las ecuaciones Rx = QT y pueden resolver inmediatamente para x por sustitución en reversa. Es decir, que una vez que se ha encontrado la descomposición QR normalizada de A , es posible resolver inmediatamente el problema de mínimos cuadrados. La prueba de este desarrollo se puede ver más a detalle en [02].

3.8 Función QR en MATLAB® . La función qr que contiene MATLAB® también se encuentra implícita en el paquete LAPACK. Esta función es útil tanto para matrices cuadradas como matrices rectangulares, la cual expresa una matriz A determinada como un producto de una matriz real ortonormal o unitaria compleja y una matriz triangular superior. La sintaxis que se obtuvo de [34], es la siguiente:

[Q,R] = qr(A), produce una matriz triangular superior R de la misma dimensión que A y una matriz unitaria Q , tal que A = QR . Para matrices “sparse”, Q es frecuentemente cercana a una matriz de rango completo. Si [m n] = size(A), entonces Q es de mxm y R es de mxn

[Q,R] = qr(A,0) , produce una descomposición “economy-size”. Si [m n] = size(A), y

m > n , entonces qr calcula solamente las primeras n columnas de Q y R es de nxn . Ahora bien, si m n entonces solo se calculan las primeras

n columnas de Q . La sintaxis para obtener esta factorización es

[Q, R ] = qr (A,0)

(3.63)

Para ejemplificar esto se aumenta un renglón a la matriz (3.59), por lo que la matriz a factorizar es 1 2 3 4 5 6  A= 7 8 0   11 12 13

usando (3.63) obtenemos las matrices

(3.64)

Qy R

 −0.0731 −0.8104 0.1667   −13.6748 −15.3567 −12.4316   −0.2925 −0.4694 0.2473  , R= 0 −1.0822 −1.0081  Q=   −0.5119 −0.1285 −0.8494   0 0 7.6445     −0.8044 0.3261 0.4354 

(3.65)

Comprobando que A = QR 1 2 3 4 5 6  A = QR =  7 8 0   11 12 13

(3.66)

como (3.66) y (3.64) son idénticas la factorización es correcta. Ahora se obtiene la factorización QR de (3.60)

58

Capítulo 3:Factoriozación QR

 −0.0731 −0.8104 0.1667 −0.5569   −13.6748 −15.3567 −12.4316   −0.2925 −0.4694 0.2473 0.7956   0 −1.0822 −1.0081  ,R = Q=  −0.5119 −0.1285 −0.8494 0.0000   0 0 7.6445      0 0 0  −0.8044 0.3261 0.4354 −0.2387   

(3.67)

Comprobando A = QR 1 2 3 4 5 6  A = QR =  7 8 0   11 12 13

como se puede observar la matriz la matriz

Q

Q

(3.68)

en (3.67) realiza todas las operaciones en contraste con

que se obtiene en (3.65) que solo calcula las primeras n columnas de

Q.

Tomando en cuenta la matriz de permutación cando existe pivoteo [Q, R, P ] = qr ( A) , donde E es la matriz de permutación, tenemos  −0.2074 −0.3994 −0.8930   −9.6437 −3.7330 −8.0882  Q =  −0.5185 −0.7293 0.4465  , R =  0 −5.5736 0.5730  0 0.5023   −0.8296 0.5556 −0.0558  0

(3.69)

0 0 1  P = 1 0 0  0 1 0 

(3.70)

2 3 1 2 3 1 AP =  5 6 4  = QR =  5 6 4  8 0 7  8 0 7 

(3.71)

Entonces comprobando que AP = QR

La primera forma de obtener la factorización QR, es el mediante el proceso de ortogonalización de Gram-Schmidt, que se puede llevar acabo de dos formas: el proceso de Gram-Schmidt clásico y el proceso de Gram-Schmidt modificado. Para esto se tomaron los programas de [28], los cuales son [Q, R ] = gs _ m( A) ,

[Q, R ] = gs _ c( A )

y

®

puesto que MATLAB no cuenta con esta factorización se considera que

puede ser útil para el caso en el que la factorización QR de la matriz propiedad de que la matriz

Q

A

no cumpla con la

resultante sea ortogonal.

Obteniendo la factorización QR de (3.59) mediante el proceso de Gram-Schmidt resulta  0.1231 0.9045 −0.4082  8.1240 9.6011 3.3235    Q = 0.4924 0.3015 0.8165  , R =  0 0.9045 4.5227  0.8616 −0.3015 −0.4082   0 0 3.6742 

59

(3.72)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

entonces comprobando que A = QR , tenemos que  1 2 3   A = QR =  4 5 6   7 8 0  

(3.73)

Ahora realizando la factorización QR mediante el proceso de Gram-Schmidt modificado  0.1231 0.9045 −0.4082  8.1240 9.6011 3.3235  Q = 0.4924 0.3015 0.8165  , R =  0 0.9045 4.5227  0 3.6742  0.8616 −0.3015 −0.4082   0

(3.74)

Se obtiene el mismo resultado que en (3.72). La segunda forma de obtener la factorización QR es mediante rotaciones de Givens, para esto se tomó de [53] el algoritmo correspondiente, debido a que MATLAB® no cuenta con alguna función que permita obtener dicha factorización. Primero se ejecuta la función

[ R , thetac, thetas] = qrgivens ( A ) ,

la cual realiza el cálculo de la

factorización QR de una matriz A sin pivoteo en columnas, mediante rotaciones de Givens.

Donde

R

8.1240 9.6011 3.3235  R =  0 0.9045 4.5227   0 0 3.6742 

(3.75)

0  0   0  0 thetac =  0.1231 0  , thetas = 0.9924 0   −0.4961 0.9115 0.8682 −0.4114 

(3.76)

es la una matriz triangular superior,

valores de

cy s

thetac y thetas proporciona

las matrices de los

que son los elementos de la matriz de Givens (3.40) necesarias para

realizar la rotación. Posteriormente se usa la función ortogonal

Q = qrmakeqgiv(thetac, thetas) para

obtener finalmente la matriz

Q.

 0.1231 0.9045 −0.4082  Q = 0.4924 0.3015 0.8165  0.8616 −0.3015 −0.4082 

Las matrices

Q

y

R

(3.77)

que se obtienen son las mismas que en (3.74), por lo que la

comprobación de A = QR es válida. La Tabla 3.1 muestra los programas que proporciona [53] para obtener la factorización QR, así como una breve descripción de ellos.

60

Capítulo 3:Factoriozación QR

Tabla 3.1Algoritmos para obtener la factorización QR mediante rotaciones de Givens.

Función

Descripción

[ R , thetac, thetas] = qrgivens ( A )

Recibe como argumento la matriz a factorizar A y regresa una matriz triangular superior R y las matrices que contienen los valores de c y s necesarias para la matriz de Givens

Q = qrmakeqgiv(thetac, thetas)

Recibe como argumento las matrices de ellas encontrar la matriz ortogonal

[B] = qrqtbgiv(B, thetac, thetas)

thetac y thetas

para que a partir

Q.

Recibe como argumentos thetac , thetas y B . Y regresa la matriz B rotada, la cual es una matriz triangular superior, útil para resolver el problema de mínimos cuadrados en el cálculo de parámetros para sistemas Ax = b y A = B = Q H .

[c, s] = qrtheta ( x, y )

Recibe como argumentos las coordenadas

x,y.

Y regresa las

matrices de valores para los elementos de la matriz de Givens

c

y s.

La tercer forma de encontrar la factorización QR es utilizando reflexiones o transformaciones Householder, en la Tabla 3.2 presentan las programas que presenta [53] junto con una descripción breve para obtener dicha factorización. Tabla 3.2 Algoritmos para obtener la factorización QR mediante transformaciones householder.

Función [ V , R ] = qrhouse( A )

Descripción Recibe como argumento la matriz

A que

se va a factorizar y regresa

un vector v householder y una matriz triangular superior Q = qrmakeq ( V )

R.

Recibe como argumento un vector householder v y regresa una matriz ortogonal Q .

H = housholder ( A , k )

Recibe una matriz A y una constante k que permite la elección de (3.80) y regresa una matriz Householder H .

[ vl, vr ] = makehouse( A )

Genera un vector Householder tal que HA tenga ceros en sus elementos excepto el primer componente, en este caso vl o vr para

61

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

aplicarlos en houseleft o houserigth según sea el caso. B = houseright ( A, v )

Recibe como argumentos la matriz A y el vector Householder v . Aplica la transformación Householder con base a v hacia la derecha de A , regresando la matriz resultante en B .

B = houseleft ( A, v )

Recibe como argumentos A y v , y aplica la transformación Householder con base a v hacia la izquierda de A , regresando la matriz resultante en B .

MATLAB® cuenta con la función house dentro de la galería de matrices de Higham, el cual en [36] presenta un análisis de precisión y estabilidad de algoritmos numéricos. Por lo que se utiliza nuevamente la matriz  1 2 3   A =  4 5 6  7 8 0  

(3.78)

para encontrar la factorización. Primero se analiza si la función house de MATLAB® proporciona correctamente la transformación Householder, escribiendo

[ v, beta, S] = gallery(' house ', A(:,1), k )

(3.79)

ésta función toma un vector A de (m,1) y regresa v y beta tal que HA = Se1 , donde Se1 es la primer columna de eye(n) ; abs(S) = norm( A) y P = I − beta * v * v ' es una matriz Householder y el parámetro k determina el signo de S k = 0 sign(S) = − sign( A(1)) k = 1 sign(S) = sign( A(1))

(3.80)

donde k = 0 es la opción por defecto. Para el caso especial en (3.79) toma cada columna de la matriz A y encuentra el vector v necesario para determinar la matriz householder utilizando (3.44) tal que multiplicando H por un vector x reducirá todos los componentes del vector a cero, excepto el elemento j-ésimo, que para este caso queda el elemento A(1,1) . De [28] se tomó un programa codificado en MATLAB® para encontrar la transformación householder dada una matriz A .

[Q, R ] = qr _ factor (A) Tomando nuevamente la matriz A en (3.78) obtenemos la matriz householder de A

62

(3.81)

Capítulo 3:Factoriozación QR

 −0.1231 −0.9045 −0.4082  Q =  −0.4924 −0.3015 0.8165  ,  −0.8616 0.3015 −0.4082 

 −8.1240 −9.6011 −3.3235  R =  0 −0.9045 −4.5227  0 3.6742   0

(3.82)

y comprobando que A = QR  1 2 3   A = QR =  4 5 6   7 8 0  

(3.83)

Concluimos que la factorización QR mediante transformaciones householder es correcta, tanto para (3.79) como para (3.81). Ahora obtenemos la factorización QR utilizando las funciones Q = qrmakeq (V ) ,

[ V , R ] = qrhouse( A )

y

tenemos

 −8.1240  −0.1231 0.9045 0.4082   0 Q =  −0.4924 0.3015 −0.8165 , R =   −0.8616 −0.3015 0.4082   0

−9.6011 −3.3235 0.9045 4.5227  0 −3.6742 

(3.84)

tomando en cuenta (3.84) comprobamos A = QR  1 2 3   A = QR =  4 5 6   7 8 0  

(3.85)

Con lo que se determina que la factorización es correcta. Las otras funciones de la Tabla 3.2 son un caso particular de las dos primeras.

3.10 Conclusiones En este capítulo se presentaron los conceptos y definiciones necesarias para la comprensión de la factorización QR dada una matriz A ya sea cuadrada o rectangular. La factorización QR produce una matriz triangular superior R y una matriz ortogonal Q tal que QQ = I . Algunos de los métodos para encontrar la factorización QR son: •

Gram-Schmidt

•

Rotaciones de Givens

•

Reflexiones o transformaciones Householder

63

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

En la práctica la factorización QR es muy útil, ya que generalmente se emplea para obtener el rango de una matriz, que es, la cantidad de elementos diferentes de cero en la diagonal de la matriz R . En la etapa de evaluación, se comprobó que la programación de las diferentes formas de obtener la factorización QR utilizando las funciones de MATLAB® y los algoritmos codificados que se obtuvieron de la literatura, es correcta; con ello en el Capítulo 6 se validarán aplicados a la identificación de procesos reales.

64

Capítulo 4 Factorización AUDI 4.1 Introducción La identificación UD Aumentada (AUDI) es una familia de nuevos algoritmos de identificación que se basan en descomposiciones matriciales que han sido estudiadas ampliamente y se presentaron en los capítulos anteriores. Este algoritmo comparado con los métodos tradicionales de mínimos cuadrados son conceptualmente más eficientes, numéricamente más robustos y con aplicaciones más completas [48]. Además de que es más simple, más flexible en su implementación y produce más información [49]. Este capítulo se organiza de la siguiente forma: En la Sección 4.2 se establecen los conceptos y definiciones que se necesitan para la comprensión y el desarrollo de esta factorización. En la Sección 4.3 se presenta la factorización UD Aumentada (AUDI), la cual se basa en el algoritmo de mínimos cuadrados recursivos y en el algoritmo UD de Bierman. En la Sección 4.4 se presenta el algoritmo AUDI-LS, el cual es útil para modelos tipo ARX. En la Sección 4.5 se aborda el algoritmo AUDI-ELS, el cual se utiliza para modelos tipo ARMAX; donde el ruido que se le agrega al sistema es modelado. En la Sección 4.6 se realiza la evaluación del algoritmo AUDI para un modelo tipo ARX con la finalidad de cotejar si el algoritmo es capaz de determinar los coeficientes de un sistema simple. Y para finalizar en la Sección 4.7 se presentan las conclusiones del capítulo.

4.2 Definiciones y conceptos generales Para este apartado se tomaron de [47], [48] y [50] las siguientes definiciones y conceptos.

65

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Definición 4.1 Factorización de mínimos cuadrados hacia delante Si se rescribe la matriz de datos aumentada A en columnas (4.1)

T

A = [ x1 , x2 ,..., xn , − y ]

Considerando el siguiente conjunto de ecuaciones simultáneas u12 a1   − a2 + e1  u a + u a  −a + e  23 2  13 1   3 2  ⋮  = ⋮      u a + u a + ... + u a − a + e 2n 2 n −1 n −1  1n 1   n n −1  θ a + θ a + ... + θ a + θ a   y + e   1 1 2 2 n −1 n −1 n n  n 

(4.2)

La ecuación anterior representa el problema original de mínimos cuadrados que se expresa en (1.1) ó (1.2). Lo demás renglones son ecuaciones adicionales. La ecuación (4.2) se puede escribir en forma matricial de la siguiente manera

[ a1

a2

a3 ⋯ an

1 u12 u13  1 u23   1 − y]     

⋯ u1n θ1  ⋯ u2 n θ 2  ⋯ u3n θ3  =E ⋱ ⋮ ⋮  1 θn   1 

(4.3)

En forma compacta AU = E

(4.4)

Las soluciones de mínimos cuadrados para todos los conjuntos de n ecuaciones lineales se pueden obtener simultáneamente mediante una simple descomposición. De acuerdo a lo anterior se plantea el siguiente teorema.

Teorema 4.2 Factorización de mínimos cuadrados hacia delante Asumiendo que la matriz de coeficientes aumentada A como se define como en (4.4) y además es no-singular, se define la siguiente factorización LDLT. S = AT A = LDLT

(4.5)

donde L es una matriz triangular inferior unitaria y D es diagonal. Entonces L y D simultáneamente proveen todas las soluciones de mínimos cuadrados y funciones de costo para los n conjuntos de ecuaciones lineales que se definieron mediante (4.2) ó (4.4) en la forma

66

Capítulo 4:Factoriozación AUDI

U = L−T

(4.6)

D=D

donde el superíndice

" −T "

es la matriz transpuesta e inversa,

U

y D son las matrices solución

y función de costo respectivamente.

Definición 4.2 Factorización de mínimos cuadrados hacia atrás Si la matriz de datos aumentada A se formula como sigue (4.7)

A = [ − y, x1 , x2 ,..., xn ]

Y si se considera los siguientes

n

conjuntos de ecuaciones

 y − e1  u12 a1   y − e  u a + u a  2    13 1 23 2  ⋮  = ⋮       y − en −1  u1n a1 + u2 n a2 + ... + un −1, n an −1   y − e  θ a + θ a + ... + θ a + θ a   n   1 1 2 2 n −1 n −1 n n

(4.8)

esto es, cada conjunto de ecuaciones usa los vectores x para ajustar el vector de observación y . La diferencia de la ecuación anterior con respecto a (4.2) radica en el número de vectores x (ordenados de 1 a n) que se utilizan para ajustar y . Si se rescribe la ecuación (4.8) en forma matricial, tenemos

[− y

a1 a2 ⋯ an −1

1 1 1 1  u u ⋯ 12 13   ⋱ an ]     

1 u1n u2 n ⋮ un −1, n

1



θ1  θ2 

=E  θ n −1   θ n 

⋮

(4.9)

Donde E es la matriz de residuos de dimensión apropiada. En forma más compacta AU = E

(4.10)

La solución de la matriz U se obtiene mediante el siguiente teorema

Teorema 4.3 Factorización hacia atrás de mínimos cuadrados Para el problema que se define en (1.2), la matriz producto de datos aumentada (no simétrica) S se define como T

S = [ − y , A ] [ A, − y ]

67

(4.11)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

y S se descompone con la factorización LDU en S = LDU , la matriz solución de mínimos cuadrados U y la matriz correspondiente a la función de costo D para los n conjuntos de ecuaciones lineales en (4.8) o (4.10) se determinan por U = Lɶ −1

(4.12)

D = diag (ET E)

donde la matriz de residuos E se calcula mediante (4.10), y Lɶ −1 es la matriz normalizada cuya inversa toma la forma de una matriz triangular superior como en (4.9). La prueba de este teorema se puede ver en [48]. Para llevar acabo la factorización LDU es necesario que las matrices L y U estén bien condicionadas [22].

Definición 4.3 Modelos múltiples de mínimos cuadrados Para las ecuaciones múltiples que se definen en (4.34), la estimación de mínimos cuadrados de la matriz solución se define como U = arg min AΘ θ

2

(4.13)

donde A ∈ ℝ mxn es la matriz de datos aumentada y puede tomar cualquier forma. Θ es la matriz solución con una estructura especial que especifica el usuario. La Figura 4.1 muestra la equivalencia del método convencional de mínimos cuadrados y el modelo múltiple de mínimos cuadrados, donde se puede observar que el método convencional es un subconjunto del modelo múltiple del problema de mínimos cuadrados.

Figura 4.1 Mínimos cuadrados vs Modelos múltiples de mínimos cuadrados [48].

Definición 4.4 Estimación de parámetros con mínimos cuadrados La aplicación del estimador de mínimos cuadrados para la estimación de parámetros se puede describir brevemente de la siguiente manera. Suponiendo que los datos de entrada/salida de un sistema están disponibles hasta un tiempo t

68

Capítulo 4:Factoriozación AUDI

(4.14)

{u (1), y (1)} ,{u (2), y(2)} ,...,{u (t ), y (t )}

donde

u(•) y y (•) son la entrada y salida del sistema respectivamente. Por lo tanto el

problema general de mínimos cuadrados para la estimación de parámetros es encontrar un modelo lineal que minimice el siguiente criterio cuadrático 2

N

(4.15)

J (t ) = ∑ [ y (t ) − yˆ (t ) ] t =1

donde se sabe que yˆ (t ) es la salida estimada del sistema. Si se supone que el sistema en cuestión se describe como un modelo de ecuación en diferencias lineal tipo ARX y (t ) + a1 y (t − 1) + ... + an y (t − n) = b1u (t − 1) + ... + bnu (t − n) + e(t )

(4.16)

o en forma compacta (4.17)

y (t ) = jT (t )θ + e(t )

donde n es el orden del modelo. e(t ) es la parte del sistema no modelada y se supone que es una señal de tipo ruido blanco con media cero y varianza σ v 2 . El vector de datos es j(t ) y el vector de parámetros θ es T

j(t ) = [ − y (t − 1),..., −y (t − n), u(t − 1),..., u (t − n)] T

θ = [ a1 ,..., an , b1 ,...bn ]

(4.18)

De tal suerte que la solución de mínimos cuadrados para este caso es equivalente a encontrar la mejor estimación θˆ (t ) para j(t )θ(t ) = y (t )

(4.19)

La ecuación anterior tiene el mismo formato que el modelo que se expresó en (criterio de minimización normal). El criterio cuadrático para (4.19) se expresa de la siguiente manera T

J (t ) = [ y (t ) − j(t )θ(t )] [ y (t ) − j(t )θ(t ) ]

(4.20)

Por consiguiente los parámetros estimados se obtienen mediante −1 θˆ (t ) =  jT (t ) j(t )  j(t )y (t )

(4.21)

La ecuación anterior representa el mismo estimador que en (ecuación normal) y de la misma forma la estimación de mínimos cuadrados para (4.21) depende de que la matriz jT (t ) j(t ) sea invertible. Con base al principio del modelo múltiple de mínimos cuadrados se establece el estimador múltiple de mínimos cuadrados, el cual tiene como principal característica el mejoramiento 69

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

del desempeño numérico debido a la utilización de técnicas de factorización y la estructura del modelo múltiple resulta del reacomodo y aumento del vector de regresión.

Definición 4.5 Método de mínimos cuadrados de modelos múltiples (MMLS). El método MMLS (por sus siglas en inglés “Multiple Model Least-Squares”) se encuentra implícito en el método de mínimos cuadrados recursivos, el cual también es un método muy antiguo al igual que el método convencional de mínimos cuadrados estándar. y el método de MMLS integra un gran número de técnicas bien establecidas, con características y propiedades favorables en una estructura muy confiable. Esta estructura proporciona más integridad de información y comodidad para la implementación en la estimación de parámetros y procesamiento de señales. Si se define que el vector de datos aumentado con un orden máximo posible n como T

j(t ) = [ − y (t − n), u (t − 1),..., −y (t − 1), u(t − 1), − y (t )]

(4.22)

Es notable que el acomodo de los datos de entrada/salida en el vector de regresión para el método convencional en (4.18) es en pares de datos entrada/salida. La salida actual del sistema en cuestión se incluye en el vector de regresión. Esta estructura es la base del método MMLS y además también es la principal diferencia entre la formulación del método MMLS y la formulación del método convencional. La matriz de datos aumentada se representa por  ϕ T (1)   T  ϕ (2)  Φ(t ) =   ⋮   T   ϕ (t ) 

(4.23)

y la matriz producto de datos análogamente a (4.5) es S(t ) = ΦT (t )Φ(t ) = ∑ j =1 j( j ) j( j )T t

(4.24)

A la matriz S(t ) se le conoce como matriz de información aumentada (AIM, por sus siglas en inglés “Augmented Information Matrix”) dentro del contexto de estimación de parámetros e identificación de sistemas. Factorizando la AIM de acuerdo al Teorema 4.2, es decir S(t ) = L(t )D(t )LT (t ) , las matrices resultantes U(t ) = L−T (t ) y D(t ) = L(t ) son la matriz de parámetros y la matriz de función de costo respectivamente para todas las ecuaciones (modelos) definidas por jT (t )U (t ) = eT (t )

70

(4.25)

Capítulo 4:Factoriozación AUDI

La cual es análogo al modelo de mínimos cuadrados en (4.17). La matriz de parámetros U (t ) tiene la forma 1 αˆ1(1) θˆ1(1) αˆ1(2) θˆ1(2)  1 θˆ2 (1) αˆ 2( 2) θˆ2 ( 2)   1 αˆ 3( 2) θˆ3(2)   1 θˆ4 ( 2) U (t ) =   1      

⋯ αˆ1( n )

θˆ1( n ) 

(n )

θˆ2( n ) 

⋯ αˆ3( n )

θˆ3( n )   θˆ4( n ) 

⋯ αˆ 2

⋯ αˆ 4 ( n ) ( n)

⋯ αˆ 5 ⋱ ⋮

1

 

(4.26)

θˆ5  ⋮   θˆ2 n ( n )  1  (n )

y la matriz de la función de costo D(t ) de la forma (4.27)

D(t ) =  J f (0) (t ), J b (1) (t ), J f (1) (t ), J b (2) (t ), J f (2) (t ),..., J b ( n ) (t ), J f ( n ) (t ) 

El superíndice

"( i )"

(i = 0,1, 2,.., n) representa el orden del modelo. Por ejemplo si θˆ 1( n )

significa que es el primer parámetro estimado del vector de parámetros para un modelo de orden nth . Los subíndices

"f"

y

"b "

son los modelos que se obtienen mediante la

factorización hacia delante (modelos hacia delante) y hacia atrás respectivamente (modelos hacia atrás). Los modelos hacia delante (forward models) usan las entradas y salidas pasadas para ajustar la salida futura. Y los modelos hacia atrás usan las entradas y salidas pasadas para ajustar las entradas futuras. Si existe retroalimentación a la salida del sistema, los modelos hacia atrás (backward models) corresponden al modelo retroalimentado y por consiguiente éstos modelos se conocen como modelos retroalimentados (feedback models). La estructura MMLS utilizar para obtener los mismos parámetros y la misma función de costo, en la Tabla 4.1 se muestran alguno de ellos. Tabla 4.1 Métodos de factorización para MMLS [48].

Método

Fórmula

Matriz de parámetros Matriz de función de costo

QR

Φ = QR

U = R −1 diag (R )

D = [ diag (R ) ]

LU

S = LU

U = L−T

D = diag (U )

LDLT

S = LDLT

U = L−T

D=D

Cholesky

S = GG T

U = G −T diag (G )

D = [ diag (G )]

71

2

2

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

UDUT

C = UDUT

D = D −1

U =U

4.3 Factorización UD Aumentada (AUDI) En la Sección 1.3.1 se presentó el problema principal que el método convencional de mínimos cuadrados tiene. Si se utilizan métodos de factorización apropiados que se apliquen a la matriz de coeficientes A , toda la información que contiene se puede extraer y expresarla explícitamente en formas útiles (factorizaciones). En la Tabla 4.2 se muestran algunos de los métodos de factorización como LU/LDLT, Cholesky, QR y UDUT que se utilizan para mejorar el desempeño numérico con respecto al método de mínimos cuadrados. Donde G es una matriz triangular inferior, L es una matriz triangular unitaria inferior, U es una matriz triangular superior unitaria, D es una matriz diagonal y Q es una matriz ortogonal. Tabla 4.2 Métodos de descomposición para la solucionar mínimos cuadrados[48] .

Método

Descomposición

Solución por mínimos cuadrados

AT A = LU

θˆ = U −1L−1 ( AT y )

Factorización LDLT

AT A = LDLT

θˆ = L−T D−1L−1 ( AT y )

Factorización Cholesky

AT A = GG T

θˆ = G −1T G −1 ( AT y )

A = QR

θˆ = R −1Q −T ( AT y )

Factorización LU

Factorización QR

El problema de mínimos cuadrados se puede presentar de una forma diferente θ  =e 1 

[ A, − y ] 

(4.28)

que es equivalente a escribir Aθ = e

(4.29)

A la forma (4.28) ó (4.29) se conoce como sistema aumentado. A es la matriz de datos aumentada y θ es el vector de solución aumentado. La matriz A contiene mucho más información que simplemente el vector de parámetros estimados θ , y esta información se puede extraer sin necesidad de realizar rutinas de cálculo extra. 72

Capítulo 4:Factoriozación AUDI

En la Tabla 4.3 muestra la relación que existe entre la factorización QR y las demás factorizaciones Tabla 4.3 Relación entre la factorización QR y las factorizaciones LDLT, LU y Cholesky[48].

Factorización LDLT

S = LDLT

Factorización LU

S = LU

Factorización Cholesky

S = GG T

Factorización QR

A = QR

U = DLT , G = LD1/ 2 , R = G T

Utilizando las ecuaciones aumentadas (4.29) el problema de mínimos cuadrados convencional se puede redefinir como sigue [48] θˆ = arg min Aθ θ

2

sujeto a θn +1 ≡ 1

(4.30)

La ecuación anterior es una forma más general de expresar el problema de mínimos cuadrados. Los siguientes tres problemas de mínimos cuadrados se pueden definir con la misma función de costo (4.30) de la siguiente manera •

Predicción de mínimos cuadrados hacia delante θˆ = arg min Aθ θ

•

sujeto a θn +1 ≡ 1

(4.31)

Predicción de mínimos cuadrados hacia atrás θˆ = arg min Aθ θ

•

2

2

sujeto a θ1 ≡ 1

(4.32)

2

sujeto a θ2 ≡ 1

(4.33)

Mínimos cuadrados completo θˆ = arg min Aθ θ

El problema convencional de mínimos cuadrados se extiende en forma de modelos múltiples sin necesidad esfuerzo computacional. Si se considera la siguiente extensión del problema de mínimos cuadrados como en la ecuación (ecuación normal).

AΘ = E

73

(4.34)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

donde A ∈ ℝ mx ( n +1) que en (4.4) se define y Θ es una matriz solución cuadrada de dimensión ((n + 1), (n − 1)). Mientras que el método convencional de mínimos cuadrados que se define en (4.4) resuelve para un conjunto de ecuaciones, la ecuación (4.34) define n conjuntos de ecuaciones y se conoce como modelos múltiples (ver la Definición 4.3). La implementación recursiva del estimador múltiple de mínimos cuadrados MMLS (ver Definición 4.4) tiene como base el algoritmo recursivo de mínimos cuadrados (RLS) que se muestra en la Tabla 4.4 el cual se basa en el estimador de mínimos cuadrados convencional. El método de Bierman es una alternativa de solución al algoritmo de mínimos cuadrados recursivos (Recursive Least-Squares), el cual se basa en el estimador de mínimos cuadrados (LS). En la Tabla 4.4 se presenta el procedimiento paso a paso del algoritmo de RLS. Tabla 4.4 Algoritmo de mínimos cuadrados recursivos (RLS)[48].

j(t ) = [− z (t − 1), ⋯ −z (t − n), u(t − 1) ⋯ u(t − n)]T zɶ (t ) = z (t ) − jT (t )θˆ (t − 1)

Vector de datos Innovación

C(t ) = C(t − 1) − C(t − 1) j(t )s (t ) j (t )C(t − 1) k (t ) = C(t ) j(t ) θˆ (t ) = θˆ (t − 1) + k (t ) zɶ (t )

Factor de escalamiento Covarianza Ganancia de Kalman Parámetros

ε (t ) = z (t ) − jT (t )θˆ (t ) J (t ) = λ (t ) J (t − 1) + zɶ (t )

Residuo Función costo

s(t ) = λ (t ) + j (t )P (t − 1) j(t ) T

−1

T

Las condiciones iniciales son C(0) = σ 2 I y θˆ (0) = 0 , donde σ es un número entero grande. A pesar de que en la teoría y en la práctica el algoritmo RLS es ampliamente aceptado su pobre propiedad numérica es muy notoria debido a que se basa en el estimador de LS. Tomando de la Tabla 4.4 el factor de escalamiento y sustituyéndolo en la ecuación de la covarianza tenemos que la matriz de covarianza se actualiza mediante la siguiente ecuación C(t ) = C(t − 1) −

C(t − 1) j(t ) jT (t )C(t − 1) λ (t ) + jT (t )C(t − 1) j(t )

(4.35)

Al mismo tiempo que se actualiza la matriz de covarianza, ésta converge a una matriz de magnitud muy pequilla, lo cual implica que la actualización de la matriz de covarianza incluya la substracción de por lo menos dos matrices iguales con magnitudes muy pequeñas. De esta manera se resuelve fácilmente el problema del pobre desempeño numérico cuando dicho algoritmo es implementado en computadores digitales con longitud de palabra finita y errores de redondeo [48].

74

Capítulo 4:Factoriozación AUDI

Existen algunos algoritmos derivados del RLS para mejorar el desempeño numérico, entre ellos se encuentra la factorización UD de Bierman, el cual es uno de los intentos más exitosos. Matemáticamente el método de Bierman es equivalente al método de RLS. Sin embargo mediante una reformulación diferente al algoritmo de Bierman, se obtiene la factorización UD que presenta propiedades mucho más estables que el algoritmo RLS. El método de Bierman además de actualizar directamente la matriz de covarianza C(t ) mediante el uso de la ecuación (4.35), se obtiene una forma factorizada de la matriz C(t ) que tiene la forma UDUT , por lo tanto C(t ) = U(t )D(t )U(t )T . A cada instante de tiempo se actualiza U(t ) y D(t ) , en lugar de la matriz C(t ) . En la Tabla 4.5 se presenta el procedimiento paso a paso de la factorización UD de Bierman. Tabla 4.5 Algoritmo para obtener la factorización UD de Bierman[48]. T

h(t ) = [ − z (t − 1), ⋯ , − z (t − n), u (t − 1), ⋯ , u (t − n) ] zɶ (t ) = z (t ) − hT (t )θˆ(t − 1) f = U T (t − 1)h(t ), g = D (t − 1), β 0 = λ (t ) X 1 ,… , X n

Para j = 1: d , hacer β j = β j −1 + f j g j D jj (t ) = D jj (t − 1) β j −1 / β j / λ (t )

µ j = − f j / β j −1 , v j = g j

Para j = 1: j − 1, hacer U ij (t ) = U ij (t − 1) + vi µi vi = vi + U ij (t − 1)v j  v1  k (t ) =  ⋮  , k (t ) = k (t ) / β d  vd 

θˆ(t ) = θˆ(t − 1) + k (t ) zɶ (t ) ε (t ) = z (t ) − hT (t )θˆ(t ) J (t ) = J (t − 1) + zɶ (t )ε (t )

A pesar de que con la factorización UD se logra mejorar las propiedades numéricas del método de mínimos cuadrados recursivos dicha factorización no se ha utilizado ampliamente. Esto ese debe principalmente al grado de complejidad para la interpretación e implementación. Además de que las matrices se utilizan como matrices intermedias durante el proceso y carecen de significado físico. Queda claro que la factorización UD se utiliza para reemplazar la actualización de la matriz de covarianza (4.35) para tener un desempeño numérico más estable. La factorización UD no solo es una técnica que mejora el desempeño numérico del método RLS, sino que 75

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

además contribuye a la creación de un nuevo concepto y un nuevo enfoque para la identificación: La factorización UD aumentada (AUDI). AUDI es un enfoque mejor que el método RLS tanto en concepto como en la implementación. Por otra parte el enfoque AUDI forma una familia completa de algoritmos correspondiente a los algoritmos de la familia original de mínimos cuadrados., y proporciona una alternativa más versátil y natural para todas las aplicaciones del RLS convencional.

4.4 Algoritmo AUDI-LS El algoritmo AUDI-LS se basa en el algoritmo UD de Bierman para llevar a cabo la estimación de parámetros en modelos de tipo ARX. Si j (t ) es el vector de datos aumentado de la forma

[ −z(t − n),

T

u(t − n), ⋯ − z (t − 1), u(t − 1), − z (t )]

(4.36)

La matriz de información aumentada (Augmented Information Matrix “AIM”) esta definida por  t  S(t ) =  ∑ j ( j ) j T ( j )   j =1  dxd

(4.37)

donde d = 2n + 1 . Para realizar identificación recursiva, además de tomar las ventajas que tiene el algoritmo de Bierman, factorización UD, la factorización UD Aumentada se centra en la inversa de la matriz de información aumentada C(t ) = S −1 (t )

(4.38)

Que se define como matriz de covarianza aumentada (Augmented Covariance Matrix “ACM”). Si denotamos la descomposición LDLT de la AIM como S(t ) = L s (t )Ds (t )LT s (t )

(4.39)

y la factorización UDUT de la ACM como C(t ) = U c (t )Dc (t )UT c (t )

(4.40)

donde los subíndices “s” y “c” indican que la descomposiciones son para AIM S(t ) y ACM C(t ) , respectivamente. a partir de C(t ) = S −1 (t ) es fácil encontrar que

76

Capítulo 4:Factoriozación AUDI

T  U(t ) = U c (t ) = L s (t )  −1  D(t ) = Dc (t ) = Ds (t )

(4.41)

La ecuación anterior demuestra que la factorización UDUT de la ACM produce la misma información U(t ) y D(t ) como la descomposición LDLT de la AIM. Ya que ACM contiene toda la información que se necesita para la identificación recursiva, la actualización recursiva de la matriz C(t ) es necesaria. Esto es (4.42)

C−1 (t ) = C−1 (t − 1) + j(t ) jT (t )

A partir del lema de inversión de matrices [48], tenemos C(t ) = C(t − 1) −

C(t − 1) j (t ) j T (t )C(t − 1) λ (t ) + j T (t )C(t − 1) j (t )

(4.43)

Haciendo una comparación entre las ecuaciones (4.43) con (4.35) se puede observar que tienen la misma estructura excepto que la matriz de covarianza C(t ) es sustituida por la matriz de covarianza aumentada C(t ) y el vector de datos j(t ) se sustituye por el vector de datos aumentado j (t ) . En la Tabla 4.6 se presenta el procedimiento paso a paso del algoritmo de identificación UD Aumentado (AUDI). Tabla 4.6 Algoritmo AUDI-LS recursivo[48]. T

ϕ (t ) = [ − z (t − n), u (t − n), ⋯ , − z (t − 1), u (t − 1), − z (t )] f = U T (t − 1)ϕ (t ), g = D (t − 1) f , β 0 = λ (t )

Para j = 1: d hacer β j = β j −1 + f j g j S D jj (t ) = D jj (t − 1) β j −1 / β j / λ (t )

µ j = − f j / β j −1 , v j = g j

Para i = 1: j − 1 , hacer (saltar j = 1 ) U ij (t ) = U ij (t − 1) + vi µi vi = vi + U ij (t − 1)v j U (t ) = U (t ), D (t ) = D −1 (t )

77

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

El algoritmo se inicializa con C(0) = U(0)D(0) UT (0) = σ 2 I donde σ es un número entero grande. Lo cual es equivalente a inicializar C(0) = σ 2 I y θˆ (0) = 0 en el algoritmo básico de RLS. Las matrices U(t ) y D(t ) son matrices temporales, matrices intermedias durante el proceso y carentes de significado físico. Para la estimación de parámetros todos los pasos presentados en la Tabla 4.4 se deben llevar a cabo. En el algoritmo AUDI, sin embargo, las matrices U(t ) y D(t ) contienen explícitamente toda la información de los parámetros estimados y la

función costo. Algunas de las ventajas del algoritmo AUDI es que omite los últimos cuatro pasos del algoritmo de Bierman porque esencialmente estos pasos repiten procedimiento previo hasta llegar al orden mas grande del modelo. Además el algoritmo AUDI es automáticamente mejorado con la inserción de los datos de salida actuales en el vector de datos aumentado. Como resultado tenemos que el algoritmo AUDI es más simple en interpretación e implementación.

4.5 Algoritmo AUDI-ELS Los algoritmos AUDI y AUDI-LS que se presentaron anteriormente se basan en la suposición de que el ruido del sistema es blanco con media cero. En el caso de sistemas con ruido coloreado2, el método AUDI-LS producirá estimación con bias, por lo tanto no puede se puede usar. Existen varias modificaciones en el principio básico del algoritmo de mínimos cuadrados que permiten remover el bias en las estimaciones realizadas. El método de mínimos cuadrados extendidos (ELS, por sus siglas en inglés) y el método de variables instrumentales (IV), son los dos métodos que se utilizan comúnmente. El método de ELS introduce una estructura pseudo lineal para estimar de forma simultánea el modelo del ruido que pueda estar agregado a la salida del sistema, y de esta manera produce una estimación libre de bias. Esta estructura se conoce como ARMAX. Suponiendo el siguiente modelo en ecuación en diferencias lineal z (t ) + a1 z (t − 1) + ... + an z (t − n) = b1u(t − 1) + ... + bn u(t − n) − e(t ) e(t ) = v(t ) + d1 v(t − 1) + ... + d n v (t − n)

2

Ruido coloreado: Alto nivel de ruido. En contraste ruido blanco. 78

(4.44)

Capítulo 4:Factoriozación AUDI

Se define el vector de datos aumentado para AUDI-ELS de la misma forma que (4.36) en AUDI-LS T

.

j(t ) = [ − z (t − n), vˆ (t − n), u(t − n),..., − z (t − 1), vˆ (t − 1), u(t − 1), − z (t )]

(4.45)

Donde vˆ ( ) se introduce en el vector de regresión, el cual es un dato estimado que permite modelar el ruido blanco agregado al sistema; éste se puede reemplazar ya sea por los residuos o por los errores de predicción [53]. La diferencia entre el vector de regresión en (4.36) con respecto a (4.45), es que en éste último el ordenamiento de los datos son en forma de tercias { z, vˆ, u} , mientras que en AUDILS el acomodo se realiza en forma de pares { z, u} . Cabe señalar que el término vˆ (t ) de ruido se determina de una manera diferente y más eficiente que en el método convencional de mínimos cuadrados extendidos, esto es (4.46)

vˆ(t ) = − µ d

Tabla 4.7 Algoritmo recursivo AUDI-ELS [48].

ϕ (t ) = [ − z (t − n), vˆ(t − n), u (t − n), ⋯ , − z (t − 1), vˆ(t − 1), u (t − 1), − z (t ) f = U T (t − 1)ϕ (t ), g = D(t − 1) f , β 0 = λ (t )

Para j = 1: d hacer β j = β j −1 + f j g j D jj (t ) = D jj (t − 1) β j −1 / β j

µ j = − f j / β j −1 , v j = g j

Para i = 1: j − 1 , hacer U ij (t ) = U ij (t − 1) + vi µi vi = vi + U ij (t − 1)v j vˆ(t ) = − µ d U (t ) = U (t ), D (t ) = D −1 (t )

79

T

]

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

4.6 Evaluación del algoritmo AUDI Para ejemplificar la aplicación de la factorización AUDI en la estimación de parámetros consideremos el siguiente sistema lineal de ecuaciones de la forma Ax = b

(4.47)

1 2 3 10    A =  4 5 6  [ x1 , x2 , x3 ] , y =  28 7 8 0  37 

(4.48)

con

Aplicando la descomposición LDLT estándar a la matriz de datos aumentada T

S = [ A, −b ] [ A, −b ] = LDLT

(4.49)

donde 78 27  66  78 93 36 S=  27 36 45   −381 −456 −198

−381 −456  −198   2253 

(4.50)

Con base a la factorización lu de la matriz A 0 0 0   1  1.1818 1.0000  0 0  L=  0.4091 5.0000 1.0000 0     −5.7727 −7.0000 −1.0000 1.0000 

(4.51)

66.0000 78.0000 27.0000 −381.0000   0 0.8182 4.0909 −5.7273  U=  0 0 13.5000 −13.5000    0 0 0  0 

(4.52)

Por lo tanto se obtienen la matriz

Uy D

1 −1.1818 5.5000 3.0000  0 1 −5.0000 2.0000  U = L−T =  0 0 1 1.0000    0 0 1  0

(4.53)

D = D = diag (66.000, 0.8182,13.500, 0.0000)

(4.54)

Como el número de parámetros a estimar son 3, por lo tanto la última columna proporciona el mejor ajuste linear de y para x1 , x2 , x3 . El error de ajuste (“loss function”) se proporciona por el último elemento de D igual a cero, ya que la solución xˆ coincide con la solución

80

Capítulo 4:Factoriozación AUDI

T

verdadera [3, 2,1] . Para evaluar la factorización cholesky es necesario replantear una matriz semidefinida positiva debido a que uno de los valores propios de (4.50) es negativo, lo cual no cumple con la propiedad por definición de matriz semidefinida positiva. Para el mismo sistema lineal (4.47) con A y b como en (4.48) obtenemos su factorización QR ya que trabaja directamente con la matriz de datos aumentada en lugar de la matriz producto de datos como lo hace LDU/LDLT , que son un caso particular de la factorización LU, y puede perder precisión a la hora de formar la matriz producto de datos debido a la longitud de palabra finita y errores de redondeo. Por lo que se recomienda la factorización QR en la implementación práctica. La factorización qr de A es  −1.6693e − 001 −2.5956e − 001 5.5205e − 001  −1.9728e − 001 −93281e − 002 −8.2808e − 001 Q=  −6.8288e − 002 9.6120e − 001 6.9007e − 002   9.6362e − 001 4.0557e − 003 −6.9007e − 002

7.7460e − 001 5.1640e − 001 2.5820e − 001  2.5820e − 001

(4.55)

 −3.9539e + 002 −4.7323e + 002 −2.0548e + 002 2.3381e + 003   0 3.8326e + 000 3.2085e + 001 −3.9750e + 001 R=  0 0 1.8632e + 000 1.8632e + 000    0 0 0 −1.6320e − 013 

Se toma la matriz

R

para obtener las matrices

Se extraen los elementos de la diagonal de renglón de

R

U

R

y

D mediante

(4.56)

el siguiente procedimiento:

para formar la matriz

D

y se divide cada

por su elemento diagonal para formar la matriz triangular superior unitaria

U.

Por lo tanto 1 −1.1969 9.5000  1 −8.3714 U=  1  

3.0000  2.0000  1.0000   1.0000 

(4.57)

D = diag (R ) = (156330,14.689,3.4714, 2.6635e − 26)

Si observamos el último elemento de la matriz diagonal

D es

(4.58)

prácticamente cero debido a

que la última columna de (4.57) proporciona los parámetros reales, por consiguiente el error de ajuste se reduce a cero. Por lo tanto se puede concluir que obteniendo la factorización UD ya sea mediante LU/LDLT o QR se puede obtener la solución a sistemas de ecuaciones lineales.

81

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

4.7 Conclusiones En este capítulo se realizó un análisis de la factorización UD Aumentada (AUDI), así como las posibles factorizaciones con las cuales se puede obtener esta factorización. La factorización AUDI es una familia de modelos estimados de manera simultánea desde un modelo de primer orden hasta un modelo de orden n , este valor es especificado por el usuario. Este algoritmo se basa en el principio del algoritmo MMLS el cual es una reformulación fundamental y una implementación del método básico de mínimos cuadrados. La ventaja de este algoritmo es que los n + 1 modelos se obtienen sin necesidad de cálculo extra, lo cual lo convierte en un algoritmo versátil y confiable. En la Figura 4.2 se muestra que el estimador múltiple de mínimos cuadrados es la base de muchas aplicaciones potenciales incluyendo identificación de sistemas.

Figura 4.2 Árbol de aplicación del estimador MMLS/AUDI [48].

El método MMLS supera ampliamente la inferioridad numérica del método convencional de mínimos cuadrados y al mismo tiempo produce más información. En la sección de evaluación se tomó por simplicidad un modelo tipo ARX para corroborar que el algoritmo estime los parámetros del modelo propuesto, teniendo un ajuste exacto en los parámetros estimados. En el Capítulo 6 se validará este algoritmo de la misma forma que la factorización LU y QR.

82

Capítulo 5 Norma L1 5.1 Introducción En este capítulo se aborda la utilización de la norma L1 como criterio de minimización aplicado a identificación de sistemas. Para ello el capítulo se organiza como sigue: En la Sección 5.2 se presentan algunos conceptos que sirven de fundamento para las siguientes secciones. En la Sección 5.3 se muestra la aplicación de este criterio a identificación de sistemas. En la Sección 5.4 se presenta el factor de predicción de error, el cual generalmente se basa en el criterio de Akaike (norma L2), utilizando la norma L1. En la Sección 5.5 Se lleva a cabo la etapa de evaluación del algoritmo estimando parámetros de un modelo sencillo tipo ARX. Por último las conclusiones se presentan en la Sección 5.6 .

5.2 Definiciones y conceptos generales Definición 5.2 Norma (o norma vectorial) Una norma en γ es una función que asigna a cada vector v en γ un número real no negativo, llamado la norma de v y está simbolizada por v (o algunas veces por v γ ) y satisface las siguientes condiciones[45]:

a)

v > 0 para v ≠ 0 y 0 = 0 , elemento no negativo.

b) α v = α v para todo escalar α y todo vector v . c)

u + v ≤ u + v (la desigualdad del triángulo) para todo vector u y v .

83

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Existen tres formas de utilizar las normas tanto en ℝ n como en ℂ n frecuentemente en las aplicaciones [02].

Definición 5.3 Norma • 1 , • 2 y •

∞

, norma 1, norma 2 y norma infinito

respectivamente, para vectores x = [ x1

x2 ⋯ x p ]T en ℝ n y ℂ n , de la siguiente

forma x 1 = x1 + x2 + ... xn 2

2

2

x 2 = ( x1 + x2 + ... xn )1/ 2 x

∞

(5.1)

= máx{ x1 + x2 + ... xn }

Para el caso de identificación la norma 1 consiste en encontrar una solución que minimice la suma del valor absoluto de la diferencia de los datos observados, es decir θˆ = arg min y − Aθ θ

(5.2)

La rutina de optimización que se utiliza para la estimación de parámetros se basa en el Método Simplex de Nelder-Mead, por lo tanto es importante conoces la siguiente definición [18].

Definición 5.4 Simplex. Un “simplex” en ℝ n es un conjunto de n + 1 puntos x1 ,..., xn en ℝ n tal que el conjunto de vectores { xi − xo :: i = 1...n} es linealmente independiente en ℝ n .

5.3 Aplicación de la norma L1 a identificación La norma L1 se basa en el valor absoluto de la diferencia de los residuos. Una de las grandes ventajas de utilizar la norma L1 es que es mucho más robusto que la norma L2 en la estimación estadística [40]. Este criterio consiste en encontrar una solución que minimice la suma del valor absoluto de la diferencia de los datos observados, es decir θˆ = arg min VN (θ, Z N ) θ

1 N ∑ ε(t , θ) N t =1 ε(t , θ) = y (t ) − yˆ (t ) VN (θ, Z N ) =

donde: θ

es el vector de parámetros a estimar.

84

(5.3)

Capítulo 5: Norma L1

ZN

es un conjunto de datos.

N

es el número de datos de estimación.

VN

es la función de costo.

ε(t , θ) es el error entre la salida real del sistema y la salida estimada. θˆ

es el vector de parámetros estimados.

A continuación se realiza una comparación entre el criterio tradicional basado en la norma L1 y el nuevo criterio que se genera de la aplicación de la norma L1 como función objetivo.

5.4 Error de predicción final (fpe) para L1 Comúnmente el criterio de error de predicción final en L2 se utiliza ampliamente para validación de modelos. Akaike describió este criterio en 1969, el cual consisten en cómo modificar la función de costo para obtener una estimación confiable de validación y del criterio de comparación únicamente a partir de los datos estimados. En [21] se presenta el desarrollo y análisis más a detalle. Este criterio se describe de la siguiente forma d 1 N ≈ d 1− N N 1+

J L2

N

∑V

N

(θˆ N , Z N )

t =1

1 VN (θˆN , Z N ) = N

N

∑ε

2

(5.4)

(t , θˆ N )

t =1

donde: d

es el número de parámetros estimados.

ZN

es un conjunto de datos.

N

es el número de datos de estimación.

VN

es la función de costo.

ε (t , θˆN )

es el error entre la salida real del sistema y la salida de modelo estimado.

θˆN

es el modelo obtenido a partir de los parámetros estimados.

Tomando el criterio (5.3) se puede escribir el error de predicción final como sigue:

85

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

d J L1 ≈ VN (θˆ N , Z N ) + N N 1 N VN (θˆ N , Z ) = ∑ ε(t , θˆ N ) N t =1

(5.5)

Donde los parámetros d , N y ε(t , θˆ N ) son los mismos para ambos criterios. Con respecto al criterio (5.5) se sabe en la comunidad científica que la estimación en L1 presenta una robustez numérica mucho más fuerte que la estimación hecha en L2, por lo que es posible el tratamiento con datos ruidosos [20].

5.5 Evaluación de la norma 1 para estimación de parámetros Para realizar la evaluación de este criterio de minimización del error que se expresa en (5.5) se considera el problema de estimación de parámetros de un sistema lineal representado por la siguiente ecuación en diferencias y (t ) − 1.5(t − 1) + 0.7(t − 2) = u (t − 1) + 0.5(t − 2) − v(t )

(5.6)

donde y (t ) y u(t ) son la estrada y salida del sistema respectivamente; v (t ) es ruido blanco con media cero y varianza σ v 2 . La señal de entrada es una secuencia binaria aleatoria de longitud N = 500 . Agregando un nivel de ruido σ v = 0.5 , se simula el sistema y se recopilan los datos de entrada/salida. Para llevar acabo este método se utiliza el toolbox de optimización con el que cuenta MATLAB®, apoyándonos principalmente de la función fminsearch, la cual utiliza el método simplex de Nelder-Mead. Esta función recibe como argumentos una función que se le conoce como “function handle”3, que para este caso se llama MinL1; El vector de regresión A,

el número de renglones con que cuenta A que se denota por k , y el vector de

observación Out _ temp de la misma longitud de A . [teta,fval,exitflag]=fminsearch(@MinL1,teta,optimset('Display','off','MaxFunEvals',7000), Out_temp,A,k); Los argumentos de salida, que proporciona la función fminsearch, de nuestro interés, son los parámetros estimados “theta,” el valor de la función objetivo al momento de terminar la minimización “fval” y la condición de salida de la función fminsearch “exitflag”, la cual puede tomar uno de tres valores 1(converge a la solución teta), 0(número máximo de evaluaciones en la función o iteraciones alcanzada) y -1(Algoritmo terminado por la función 3

Function hanndle: Es un valor que proporciona MATLAB a una función para llamarla indirectamente 86

Capítulo 5: Norma L1

de salida). Una vez que se lleva acabo la minimización los parámetros estimados para el sistema (5.6) que arroja la función fminseach son T

teta = [ −1.4821,0.6792,1.0803,0.59149]

(5.7)

Donde a1 = −1.4821, a2 = 0.6792, b1 = 1.0803, b2 = 0.59149 . Por consiguiente los parámetros estimados comparados con los parámetros reales se pude concluir que la estimación es buena.

5.6 Conclusiones En este apartado se presentó el uso de la norma L1 como criterio de minimización del error existente entre el vector de datos observados de un sistema y el vector de datos estimados. El uso de este criterio generalmente es en el área de control robusto, ya sea para la sintonización o diseño de controladores debido su propiedad de robustez numérica. Por ello se plantea el uso de la norma L1 a la identificación de sistemas, ya que los datos con los que se trabajan siempre presentan ciertas anomalías, principalmente ruido. En la sección de evaluación del algoritmo, éste obtuvo un buen ajuste con respecto a los parámetros del modelo propuesto. En el siguiente capítulo se validarán todos los algoritmos que se expusieron anteriormente con datos reales recolectados de forma experimental, que se tomaron de la base de datos Daisy.

87

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

88

Capítulo 6 Validación de algoritmos 6.1 Introducción En este capítulo, el proceso de identificación se lleva a cabo con base en los puntos que se expusieron en la Sección 1.9.1. De acuerdo a los criterios para la selección y validación del modelo, que se presentaron en la Sección 1.9.2, nos basaremos en el análisis de complejidad del modelo que tiene como principio básico el criterio de Akaike (FPE), tanto para L2 como para L1. Para llevar acabo la validación nos basaremos en el criterio de coherencia en el comportamiento de entrada-salida mediante las gráficas de Bode y el criterio de simulación, esto para los casos de estudio 2-4. También se utilizará el criterio del menor tiempo de ejecución en la estimación de los parámetros, tomando en cuenta el tiempo de cálculo de los mismos. Para el primer caso de estudio el criterio de validación se basará en la comprobación de parámetros fijos, ya que de este proceso se tiene conocimiento de los parámetros reales. En el Apéndice A se describen brevemente los procesos que se van a utilizar para validar los algoritmos LU, QR, AUDI y la aplicación de la norma L1 en la estimación de parámetros. Estos procesos son: 1. Sistema de tercer orden 2. Banco de nivel RCN100 que se encuentra en el laboratorio de control en CENIDET. 3. Secadora de pelo. 4. Sistema de calentamiento.

89

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Este capítulo se divide en tres secciones. En la Sección 6.2 se describen las dos estructuras que se utilizarán para la identificación del modelo. La Sección 6.3 consiste en la etapa de identificación de cada uno de los procesos seleccionados. Por último en la Sección 6.4 se presentan las conclusiones del capítulo.

6.2 Estructura del modelo Dentro del contexto de identificación de sistemas, en este caso, sistemas SISO, lineales determinísticos, existen estructuras bien definidas que algunas de ellas se nombran en la Tabla 6.1. Aquí nos enfocaremos a utilizar dos tipos de estructuras ARX y ARMAX. Partiendo de la que la estructura a bloques del modelo general es:

Figura 6.1 Esquema de bloques del modelo general

La Figura 6.1 muestra la relación entre el modelo general y los casos especiales que existen conocidos como caja negra cuya estructura representativa es A( q ) y (t ) =

B(q) C (q) u (t ) + e( t ) F (q) D(q )

(6.1)

Tabla 6.1 Casos especiales más comunes de (6.1) [53].

Polinomios utilizados B AB ABC AC ABD ABCD BF BFCD

Nombre de la estructura del modelo FIR ARX ARMAX ARMA ARARX ARARMAX OE BJ

La diferencia entre los modelos que se describen la tabla anterior radica en tipo de modelo con el cual puede ser representada la perturbación que se agrega al sistema. Todos partiendo del modelo ARX.

90

Capítulo 6: Validación de algoritmos

6.2.1 Estructura ARX. El acrónimo ARX se refiere a: AR, parte auto regresiva A(q)y (t ) y X a la entrada exógena B (q )u(t ) .

El caso especial donde na = 0 , z (t ) se modela como una respuesta al impulso finita

(FIR)[53]. El modelo auto regresivo se denota por (6.2)

A(q −1 )y (t ) = B (q −1 )u(t − nk ) + e(t )

La representación en ecuación en diferencias de (6.2) es y (t ) + a1 y (t − 1) + ... + ana y (t − na ) = b1u (t − 1) + ... + bnb u (t − nb ) + e(t )

(6.3)

El término del ruido blanco e(t ) se introduce como un error directo en la ecuación en diferencia, el modelo (6.3) se conoce como modelo de ecuación de error. Los parámetros que definen esta estructura están representados por θ =  a1

a2 ⋯ ana

b1 ⋯ bnb 

T

(6.4)

Los polinomios A y B se describen como A(q) = 1 + a1q −1 + ... + ana q − na

(6.5)

B(q) = b1q −1 + ... + bnb q nb

(6.6)

Lo que permite obtener el predictor (6.7)

yˆ (t | θ) = B (q)u(t ) + [1 − A( q)] y (t )

Introduciendo el vector de regresión T

j(t ) = [ − y (t − 1)... − y (t − na ) u (t − 1)...u (t − nb )]

(6.8)

Con lo que finalmente se puede reescribir (6.2) como yˆ (t | θ) = θT j(t ) = j(t )T θ

(6.9)

El predictor (6.9), el cual es equivalente a (1.7), es un producto escalar entre un vector de datos conocido j(t ) y el vector de parámetros θ . En estadística, los modelos estimados con este tipo de estructura se conoce como regresión lineal, y el vector j(t ) se conoce como vector de regresión [53]. Para este caso algunos de los coeficientes de los polinomios A y B se conocen, por lo que el modelo (6.9) en la forma de regresión lineal es yˆ (t | θ) = jT (t )θ + µ (t )

91

(6.10)

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

donde µ (t ) también es un término conocido.

6.2.2 Estructura ARMAX. Una de las desventajas que presenta la estructura (6.3) es que no describe de manera adecuada el término de la perturbación e(t ) . Esto se puede solucionar describiendo la ecuación de error como un ruido blanco de media móvil [53], lo cual proporciona un modelo como y (t ) + a1 y (t − 1) + ... + ana y (t − na ) = b1u (t − 1) + ...

+bnb u (t − nb ) + e(t ) + c1e(t − 1) + .. + cnc e(t − nc )

(6.11)

donde (6.12)

C (q) = 1 + c1q −1 + ... + cnc q nc e(t )

Ahora el vector de parámetros a estimar es θ =  a1 ⋯ ana

b1 ⋯ bnb

c1 ⋯ cnc 

T

(6.13)

el predictor obtenido para esta estructura se representa por yˆ(t | θ ) =

 A(q )  B(q ) u (t ) + 1 −  y (t ) C ( q)  C (q) 

(6.14)

Si se introduce el error de predicción ε (t , θ) = y (t ) − yˆ (t | θ)

(6.15)

j(t ,θ ) = [ − y (t − 1)... − y (t − na ) u (t − 1)...u (t − nb ) ε (t − 1,θ )...ε (t − nc ,θ ) ]

(6.16)

y el vector de regresión

el modelo estimado mediante la estructura ARMAX se obtiene mediante yˆ (t | θ) = jT (t , θ)θ

(6.17)

El modelo anterior es muy similar al modelo (6.9) de regresión lineal, pero el modelo (6.17) no es una regresión lineal debido al efecto no lineal de θ en el vector j(t ,θ ) por lo que a este tipo de modelos se le considera como regresión pseudo lineal.

92

Capítulo 6: Validación de algoritmos

6.3 Identificación de los procesos La razón por la cual los procesos que a continuación se analizan, a excepción del sistema de tercer orden, es porque los datos proporcionados se obtuvieron de forma experimental por lo que se considera un buen criterio para validar los algoritmos propuestos aplicados a identificación de sistemas. En esta sección se validan los algoritmos propuestos aplicados a la identificación de sistemas. Para realizar la validación, se cuenta con datos experimentales de los procesos: RCN100, secadora de pelo y sistema de calentamiento. Primero se identificará un proceso sintético que se tomó de [25]. Después se identificarán procesos cuyos datos se obtuvieron de [10]. La selección del periodo de muestreo para este primer caso de estudio se describe en el Apéndice B. El proceso de identificación se basa en el procedimiento que mostró en la Sección 1.9.1.

6.3.1 Sistema de tercer orden El primer caso de estudio es un sistema de tercer orden donde se conocen perfectamente sus parámetros. El objetivo de utilizarlo para realizar el proceso de identificación con los algoritmos propuestos es para verificar que éstos determinan los valores de los parámetros del sistema. El criterio que se utilizará para validar los algoritmos, en este primer caso, se basará en la comparación de parámetros fijos.

6.3.1.1 Estructura ARX La estructura en el dominio del tiempo del proceso es G (s) =

s + 10 s 3 + 6 s 2 + 9s + 10

(6.18)

Para discretizar el sistema es necesario determinar algunos parámetros con base a la respuesta transitoria. En el Apéndice B se determina el tiempo de muestreo que para este caso es Ts = .246s . Utilizando un retenedor de orden cero y con el tiempo de muestreo determinado se procede a discretizar el sistema (6.18) mediante la función C2D de MATLAB® que permite convertir la función transferencia de un sistema en tiempo continuo a su representación correspondiente en tiempo discreto. La sintaxis es la siguiente Gz = C2D(Gs,Ts,'zoh') Por lo tanto se tiene 93

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

G( z) =

0.03029 z 2 + 0.03402 z − 0.003445 z 3 − 1.978 z 2 + 1.297 z − 0.2577

(6.19)

La señal de excitación que se utiliza como señal entrada es una señal binaria pseudo aleatoria (SBPA), con una longitud N = 500 que se genera utilizando la función idinput, la cual sirve para generar señales de entrada para tareas de identificación. La sintaxis es la siguiente sbpa=idinput(N,'prbs'); donde el argumento prbs indica que es una señal binaria pseudo aleatoria. Una vez que se tiene la función transferencia en tiempo discreto y la señal de entrada, se aplica ésta última al sistema (6.19) mediante el comando filter, con la siguiente sintaxis yplant=filter(numplant,denplant,sbpa); donde num y den.son los parámetros del polinomio del numerador y denominador de (6.19) respectivamente. Entonces la señal de entrada sbpa se filtra a través de un filtro que se forma con los vectores num y den para crear la señal de salida yplant.. En la Figura 6.2 se muestra las señales de entrada/salida obtenidas Datos de Entrada/Salida

Magnitud

1 0.5 0 -0.5 -1 0

50

100

150

200 250 300 Numero de muestras

350

400

450

500

50

100

150

200 250 300 Numero de muestras

350

400

450

500

Magnitud

0.5 0 -0.5 -1 0

Figura 6.2 Datos de entrada/salida para el sistema de tercer orden.

Con los datos de entrada/salida y el tiempo de muestreo correspondiente se genera un modelo no paramétrico con la función iddata MODELO=iddata(Out,In,Ts); La Figura 6.3 muestra la respuesta transitoria del MODELO, que corresponden a la respuesta al escalón e impulso del modelo no paramétrico. La respuesta en frecuencia (bode) se muestra en la Figura 6.4.

94

Capítulo 6: Validación de algoritmos

A estas gráficas se conocen como estimaciones no paramétricas las cuales nos sirven como referencias para determinar el modelo óptimo durante la estimación paramétrica. Respuesta al escalon del modelo no paramétrico (REAL)

Amplitud

Amplitud

10

1

0.5 REAL 0 0

1

2

3

4 Tiempo (s)

5

6

-2

10

10

8

0

0.5

1

1.5

2

2.5

0

1

SPA Fase (grados)

Amplitud

-1

10

-3

7

Respuesta al impulso del modelo no parampetrico (REAL)

0.5

0

-0.5 0

BODE

0

1.5

1

2

3

4 Tiempo (s)

5

6

7

8

Figura 6.3 Respuesta al escalón y al impulso del sistema de tercer orden.

-100 -200 -300 -400 0

0.5

1

1.5 Frecuencia (Hz)

2

2.5

Figura 6.4 Respuesta en frecuencia del sistema de tercer orden.

La Tabla 6.2 muestra el porcentaje de desviación de los parámetros estimados con respecto a los parámetros reales del sistema, este porcentaje de desviación se determina para cada algoritmo. Tabla 6.2 Porcentaje de desviación de los parámetros estimados para el sistema de tercer orden con estructura ARX.

a1

ARX [3 3 0] LS [3 3 0] LU [3 3 0] QR [3 3 0] AUDI [3 3 0] 3.994% 0.0% 0.0% 0.0% 0.0%

L1 [3 3 0] 20.020%

a2

10.409%

0.0%

0.0%

0.0%

0.0%

52.660%

a3

24.253%

0.0%

0.0%

0.0%

0.0%

76.193%

b1

99.970%

0.0%

0.0%

0.0%

0.0%

89.993%

b2

11.022%

0.0%

0.0%

0.0%

0.0%

-34.891%

b3

-5.660%

0.0%

0.0%

0.0%

0.0% -278.519%

Es claro que el modelo estimado con la función arx de MATLAB® no encuentra los verdaderos parámetros, no obstante los modelos obtenidos con LS, LU, QR y AUDI tienen una aproximación exacta. Por lo tanto para un modelo con estructura tipo ARX los algoritmos más precisos son LS, LU, QR y AUDI. La Tabla 6.3 muestra los tiempos de ejecución que cada uno de los algoritmos se tomó para estimar los parámetros que se muestran en la Tabla 6.2. El tiempo de ejecución se obtuvo utilizando la función etime de MATLAB® mediante la siguiente sintaxis 95

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

t0=clock; operación, tf=etime(clock,t0); donde t0 es el tiempo inicial el cual se inicializa según el reloj del sistema, operación es la rutina o subrutina que se va a desarrollar y tf es el tiempo final, el cual se calcula mediante la resta del tiempo transcurrido y el tiempo inicial. Es importante mencionar que el tiempo de ejecución no se mantiene constante, aunque idealmente se desearía eso; pero las funciones con las que cuenta MATLAB® se basan en el tiempo que ocupa la Unidad Central de Proceso (CPU, por sus siglas en inglés Central Processing Unit) en llevar acabo la operación, existe la posibilidad de que si se tienen procesos ejecutándose en segundo plano afecte el resultado que arroje la función etime. Este tiempo de ejecución (TE) se determinó desactivando los servicios de Microsoft y los programas que se ejecutaban en segundo plano en ese momento para que la medición fuera lo más confiable posible. A pesar de que la rutina arx es pre-compilada y que los programas elaborados para cada algoritmo son interpretados, éstos últimos tienen un tiempo de ejecución menor. LU1 se refiere al algoritmo LU utilizando la función lu de MATLAB®, LU2 es la factorización LU mediante el método Cholesky (2.13), LU3 es la factorización LU utilizando el método Crout (2.12) y LU4 es la factorización LU con la forma Doolittle (2.11). Tabla 6.3 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el sistema de tercer orden con estructura ARX.

ARX LS LU1 LU2 LU3 LU4 AUDI L1 TE 0.4210s 0.0700s 0.0300s 0.0300s 0.0300s 0.0300s 0.4210s 0.3110s En la Tabla 6.4 se presentan los tiempos de ejecución para el algoritmo QR donde QR1 es la factorización QR utilizando la función de MATLAB®, QR2 es la factorización QR utilizando el proceso de Gram-Schmidt clásico, QR3 es la factorización QR utilizando el proceso de Gram-Schmidt modificado, QR4 es la factorización QR utilizando rotaciones de Givens y QR5 es la factorización QR utilizando transformaciones o reflexiones de Householder. Los cuatro métodos de obtener la factorización QR proporcionan los mismos parámetros, por lo que se muestra en la Tabla 6.2 el tiempo de ejecución (TE) que lleva cada algoritmo para estimar los parámetros correspondientes incluyendo la factorización LU, AUDI y L1. Tabla 6.4 Tiempo de ejecución de QR vs ARX, LS, AUDI y L1 para el sistema de tercer orden con estructura ARX.

ARX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 0.4210s 0.0700s 0.0100s 0.0100s 0.2110s 0.3510s 0.2500s 0.0300s 0.3110s

96

Capítulo 6: Validación de algoritmos

En la Figura 6.5 se presentan las estimaciones paramétricas para cada algoritmo, primeramente se muestran las gráficas correspondientes a la respuesta escalón e impulso. Enseguida se tienen las gráficas de BODE y por último una gráfica comparativa entre las salidas de los modelos estimados y la salida del proceso real.

Respuesta al escalon de los modelos obtenidos

Respuesta al impulso de los modelos obtenidos

1.2

0.8 0.7

REAL

0.6

LS ARX

1

LU

Amplitud

0.6 REAL LS

0.4

L1

0.1 0

L1

0

QR AUDI

0.4

0.2

LU QR AUDI

0.2

0.5

0.3

ARX

-0.1

1

2

3

4 Tiempo (s)

5

6

7

8

-0.2 0

a) Respuesta al escalón.

1

2

3

4 Tiempo (s)

BODE

2

REAL LS ARX

0

10

LU QR

-2

10

-4

10

0

0.5

1

1.5

2

2.5

0.5

1

1.5 Frequencia (Hz)

2

2.5

0 -100 -200 -300 -400 0

5

b) Respuesta al impulso.

10

Amplitud

-0.2 0

Fase (grados)

Amplitud

0.8

c) Respuesta en frecuencia 1.

97

6

7

8

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. BODE

2

10

0

Amplitud

10

-2

10

-4

10

0

0.5

1

1.5

2

2.5

0 Fase (grados)

REAL -100

AUDI L1

-200 -300 -400 0

0.5

1

1.5 Frecuencia (Hz)

2

2.5

d) Respuesta en frecuencia 2. Salida del sistema real y salida de los modelos estimados 1 Measured Output IDLS Fit: 100% IDARX Fit: 98.89%

Magnitud

0.5

IDLU Fit: 100% IDQR Fit: 100% IDAUDI Fit: 100% IDL1 Fit: 94.97%

0

-0.5

-1

-1.5 0

20

40

60 Tiempo (s)

80

100

120

e) Comparación de las salidas de los modelos estimados. Figura 6.5 Estimaciones paramétricas para el sistema de tercer orden para un modelo tipo ARX.

Con respecto a las gráficas anteriores, los modelos LS, LU, QR y AUDI presentan un mejor ajuste en los parámetros del sistema (6.19). Los modelos que se obtuvieron con ARX y L1 son los que tienen mayor porcentaje de desviación de los parámetros reales. Todos los algoritmos estimaron los parámetros del sistema real con un tiempo de ejecución mucho menor que ARX. Sin embargo todos los modelos estimados son los suficientemente representativos de la dinámica del sistema real por lo que cualquiera de ellos puede ser buena opción.

98

Capítulo 6: Validación de algoritmos

6.3.1.2 Estructura ARMAX Para estimar los modelos tipo ARMAX correspondiente a cada algoritmo se realiza el mismo procedimiento de identificación. La Tabla 6.5 muestra los parámetros obtenidos para los modelos estimados. Donde los modelos LS, LU, QR son los que obtienen los parámetros exactamente igual al sistema real. Y los modelos ARMAX, AUDI y L1 presentan una desviación notable en la mayoría de los parámetros estimados. Tabla 6.5 Porcentaje de desviación de los parámetros estimados para el sistema de tercer orden con estructura ARMAX. ARMAX [3 3 1 0] LS[3 3 1 0] LU [3 3 1 0] QR [3 3 1 0] AUDI [3 3 1 0] L1 [3 3 3 0]

a1

4.500%

0.0%

0.0%

0.0%

5.763%

25.278%

a2

11.950%

0.0%

0.0%

0.0%

15.111%

63.778%

a3

28.211%

0.0%

0.0%

0.0%

35.351%

60.962%

b1

99.994%

0.0%

0.0%

0.0%

0.660%

-25.619%

b2

11.229%

0.0%

0.0%

0.0%

-10.111%

-43.357%

b3

-967.343%

0.0%

0.0%

0.0%

-967.343%

-416.98%

En la Tabla 6.6 se muestran los tiempos de ejecución de ARMAX, LS, los diferentes métodos de obtener LU, AUDI y L1. Siendo mucho muy superior los algoritmos LS, LU, AUDI y L1 con respecto a la función ARMAX. Tabla 6.6 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el sistema de tercer orden con estructura ARXMAX. ARMAX LS LU1 LU2 LU3 LU4 AUDI L1 TE 1.5120s 0.0600s 0.0600s 0.0300s 0.0300s 0.0300s 0.8420s 0.4500s

La Tabla 6.7 presenta los tiempos de ejecución para los algoritmos ARMAX, LS, los diferentes métodos de obtener QR, AUDI y L1. Se puede observar que los tiempos de ejecución siguen siendo menor que ARMAX al igual que la familia de algoritmos LU. Tabla 6.7 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para el sistema de tercer orden con estructura ARXMAX. ARMAX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 1.5120s 0.0600s 0.0600s 0.0600s 0.0600s 0.4410s 0.3210s 0.8420s 0.4500s

En la Figura 6.6 se reúnen las gráficas correspondientes a los modelos estimados con los algoritmos LU, QR, AUDI y L1 con estructura ARMAX. El modelo que presenta menos ajuste al modelo real es el que se obtuvo con L1 en el análisis temporal (Figura 6.6.a), b)). En la respuesta en frecuencia Figura 6.6 c) y d) L1 sigue siendo el modelo, que a comparación de los restantes, que menos ajusta al modelo del sistema real. Simulando los modelos estimados con respecto al sistema real con la misma señal de entrada ambos, se obtiene el porcentaje de ajuste de la señal de salida de cada modelo

99

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

(Figura 6.6 e)), siendo los modelos LU y QR los modelos más exactos representando en un 100% la dinámica del modelo real. Respuesta al escalon de los modelos obtenidos

Respuesta al impulso de los modelos obtenidos

1.2

0.8 0.7

REAL

0.6

LS ARMAX

1

LU

0.4

Amplitud

0.6

REAL

0.5

QR AUDI

0.4

L1

0.3 0.2

LS ARMAX LU QR

0.2

0

1

2

3

4 Tiempo (s)

0

5

6

-0.1

7

-0.2 0

8

1

a) Respuesta al escalón.

2

3

4 Tiempo (s)

BODE

2

0

10

-2

10

-4

10

0

0.5

1

1.5

2

2.5

0 REAL LS ARMAX LU

-100 -200

QR -300 -400 0

0.5

5

6

b) Respuesta al impulso.

10

Amplitud

-0.2 0

0.1

AUDI L1

Fase (grados)

Amplitud

0.8

1

1.5 Frecuencia (Hz)

2

c) Respuesta en frecuencia 1.

100

2.5

7

8

Capítulo 6: Validación de algoritmos

BODE

5

Amplitud

10

0

10

-5

10

0

0.5

1

1.5

2

2.5

0 Fase (grados)

REAL -100

AUDI L1

-200 -300 -400 0

0.5

1

1.5 Frecuencia (Hz)

2

2.5

d) Respuesta en frecuencia 2. Salida del sistema real y salida de los modelos estimados 1 Measured Output IDLS Fit: 100% IDARMAX Fit: 98.23%

Magnitud

0.5

IDLU Fit: 100% IDQR Fit: 100% IDAUDI Fit: 98.33% IDL1 Fit: 90.38%

0

-0.5

-1

-1.5 0

20

40

60 Tiempo (s)

80

100

120

e) Comparación de las salidas de los modelos estimados. Figura 6.6 Estimaciones paramétricas para el sistema de tercer orden para un modelo tipo ARMAX.

6.3.2

Banco de nivel (RCN 100)

El segundo caso de estudio es un banco de nivel RCN 100 que se encuentra en el laboratorio de control automático en CENIDET. En el Apéndice A.2 se realiza una descripción breve del proceso. En la sección 6.3.2.1 se identifica el proceso tomando como premisa que su estructura es de tipo ARX. La característica principal de este proceso es que tiene una dinámica rápida. El criterio que se utilizará a partir de este caso de estudio se basa principalmente en la coherencia en el comportamiento de los datos de entrada/salida porque no se conoce la estructura real del sistema y por ende los valores reales de los parámetros a estimar.

101

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Para poder determinar los modelos para cada algoritmo primeramente se realizó la identificación paramétrica con las funciones que proporciona el Toolbox de identificación de MATLAB®: arx y armax. ¿Cómo se realizó esto? Para poder determinar el orden óptimo de un modelo es necesario realizar una barrido en los órdenes que definen la estructura ARX, na, nb y nk. Primero se realiza un barrido en el orden del retardo, dejando fijos na y nb con un valor elevado, por ejemplo 10. Posteriormente efectúa el barrio de 1-10 en nk y se determinan los valores mínimos de cada modelo estimado, extrayendo la “Loss Function” o la “FPE” de cada modelo. Teniendo así, para este caso, 10 candidatos a ser el orden del retardo. Con base a la magnitud de los mínimos obtenidos para cada modelo se determina que el orden que tenga menor magnitud sea el primer candidato. Una vez elegido el orden óptimo para nk, se deja fijo el valor de nb y se efectúa el barrido en na. De la misma forma que se eligió el orden para nk se obtiene el orden para nb y una vez teniendo el orden de nb se determina el orden de na. Finalmente se obtiene el modelo óptimo con los órdenes elegidos y se generan las respuestas al escalón e impulso, para el análisis temporal, y la respuesta en frecuencia del modelo estimado. Las respuestas del modelo estimado se comparan con las respuestas que se tienen del modelo real para determinar si caracteriza lo suficientemente bien la dinámica del proceso. En caso de que el modelo no sea bueno se vuelve a realizar el cálculo de los mínimos para cada orden hasta que se obtenga el modelo óptimo, como se muestra en el diagrama de flujo de la Figura 1.1. Una vez que se obtuvieron los órdenes óptimos para cada proceso, éstos se utilizaron para validar los algoritmos alternativos.

102

Capítulo 6: Validación de algoritmos

6.3.2.1 Estructura ARX La Figura 6.7 muestra los datos de entrada/salida. Datos de Entrada/Salida

Magnitud

1

0.5

0 0

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

Magnitud

15

10

5 0

Figura 6.7 Datos de entrada/salida del proceso banco de nivel RCN100.

Las perturbaciones en baja frecuencia, los desplazamientos (offsets), las tendencias, etc., son comunes en los datos que se obtienen de forma experimental. Existen dos enfoques para tratar estos problemas: Remover las perturbaciones por medio de pretratamiento de datos ó dejar que la parte encargada de modelar el ruido se encargue de estas perturbaciones en el caso de estructuras como ARMAX. Para este caso se hace un pretratamiento a los datos de entrada/salida utilizando el comando dtrend, del Toolbox de identificación de MATLAB®, para eliminar tendencias; el cual consiste en eliminar el mejor ajuste en línea directa de un vector determinado, para este caso tanto para el vector de entrada como el de salida. En la Figura 6.8 se muestran los datos de entrada/salida después de eliminar tendencias en los vectores. Datos de Entrada/Salida

Magnitud

0.5

0

-0.5 0

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

Magnitud

10

5

0

-5 0

Figura 6.8 Datos de entrada/salida después de eliminar tendencias del proceso banco de nivel RCN100.

103

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

En la Figura 6.9 se muestran las estimaciones no paramétricas del proceso tanto la respuesta al escalón como la respuesta al impulso, y en la Figura 6.10 se muestran las gráficas de Bode. De la misma forma que el proceso anterior éstas gráficas se utilizarán para la validación de los modelos estimados para cada estructura: ARX y ARMAX. A diferencia del proceso anterior, los datos se obtuvieron de forma experimental. Respuesta al escalon del modelo no paramétrico (REAL)

10

30

Amplitud

Amplitud

40

20 REAL

10 0 0

50

100

150 Tiempo (s)

200

250

10

10

300

Respuesta al impulso del modelo no parampetrico (REAL)

1

0

-1

0

0.4

0

0.3

-50

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

SPA Fase (grados)

Amplitud

10

BODE

2

0.2 0.1 0 0

50

100

150 Tiempo (s)

200

250

-100 -150 -200 0

300

Figura 6.9 Respuesta al escalón y al impulso del banco de nivel RCN100.

0.01

0.02

0.03

0.04 0.05 0.06 Frecuencia (Hz)

0.07

0.08

0.09

0.1

Figura 6.10 Respuesta en frecuencia del banco de nivel RCN100.

Existen diferentes alternativas de medir la diferencia entre dos señales. A partir de este proceso se utilizaran dos índices de desempeño para medir la precisión de los modelos estimados para cada algoritmo. Estos dos criterios son: Error cuadrático medio, el cual representa la desviación entre dos señales, definido por [42] n

∑(y ECM =

i

− yˆ i )

i =1

n(n − 1)

2

(6.20)

donde: yi

es la salida del proceso real

yˆ i

es la salida del modelo estimado

n

número de muestras

El error cuadrático medio nos da la medida de las diferencias en promedio entre los valores reales y los estimados. Otra medida que proporciona información similar es el error absoluto medio, definido por [38]: 104

Capítulo 6: Validación de algoritmos

n

∑y EAM =

i

− yˆ i

i =1

n(n − 1)

(6.21)

En la Tabla 6.8 se muestran los criterios ECM y EAM para las respuestas al escalón (Esc,) e impulso (Imp.) de cada modelo estimado. También se muestran los órdenes óptimos para cada algoritmo. Tabla 6.8 ECM y EAM de los modelos estimados para el proceso banco de nivel RCN100 con estructura ARX. ARX [2 1 1] LS [2 1 1] LU [2 1 1] QR [2 1 1] AUDI [2 2 1] L1 [2 1 1] Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. ECM 0.0968 0.0011 0.1001 0.0011 0.0968 0.0011 0.0968 0.0011 0.1002 0.0010 0.0950 0.0011 EAM 0.0122 0.0001 0.0126 0.0001 0.0122 0.0001 0.0122 0.0001 0.0133 0.0001 0.0118 0.0001

El tiempo de ejecución que consumió la familia de algoritmos LU, y los algoritmos ARX, LS, AUDI y L1 se presenta en la Tabla 6.9. Los tiempos de los algoritmos alternativos siguen siendo menor que el algoritmo ARX.

TE

Tabla 6.9 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el banco de nivel RCN100 con estructura ARX. ARX LS LU1 LU2 LU4 LU4 AUDI L1 0.4400s 0.0400s 0.0800s 0.0600s 0.0400s 0.0400s 0.0300s 0.3700s

Por otro lado la Tabla 6.10 se muestran los tiempos de ejecución para la familia de algoritmos QR y los algoritmos ARX, LS, AUDI y L1. Aquí se puede observar que los algoritmos QR4 y QR5 son los únicos que superan el tiempo de ejecución de ARX. Esto se debe al número de parámetros a estimar así como la longitud de los vectores entrada-salida.

TE

Tabla 6.10 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para el banco de nivel RCN100 con estructura ARX. ARX LS LU QR1 QR2 QR3 QR4 QR5 AUDI L1 0.4400s 0.0400s 0.0800s 0.0100s 0.0100s 0.0100s 1.3020s 1.1310s 0.0300s 0.3700s

En la Figura 6.11 se muestran las estimaciones paramétricas para cada algoritmo. Los modelos estimados presentan buen ajuste tanto en la respuesta temporal como en la respuesta en frecuencia, ver Figura 6.11 a)-d). En la Figura 6.11 e) se muestra un bajo porcentaje de ajuste con respecto a la salida del sistema real, a pesar de eso los modelos que se obtuvieron se consideran representativos de la dinámica del proceso real.

105

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. Respuesta al impulso de los modelos obtenidos

Respuesta al escalon de los modelos obtenidos 35

0.35

30

0.3

25

0.25

20

0.2

REAL LS

Amplitud

15

REAL LS

LU QR AUDI L1

0.15

ARX

10

0.1

LU QR 5

0.05

AUDI L1 50

100

150 Tiempo (s)

200

250

0 0

300

50

a) Respuesta al escalón.

100

150 Tiempo (s)

BODE

2

Amplitud

1

10

0

10

-1

10

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Fase (grados)

0 -50

REAL LS

-100

ARX LU QR

-150 -200 0

0.01

0.02

0.03

0.04 0.05 0.06 Frecuencia (Hz)

0.07

0.08

0.09

0.1

0.09

0.1

c) Respuesta en frecuencia 1. BODE

2

10

1

10

0

10

-1

10

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0 REAL -50

AUDI L1

-100 -150 -200 0

0.01

0.02

200

b) Respuesta al impulso.

10

Amplitud

0 0

Fase (grados)

Amplitud

ARX

0.03

0.04 0.05 0.06 Frecuencia (Hz)

0.07

0.08

d) Respuesta en frecuencia 2.

106

0.09

0.1

250

300

Capítulo 6: Validación de algoritmos

Salida del sistema real y salida de los modelos estimados 12 Measured Output IDLS Fit: 71.72% IDARX Fit: 71.7% IDLU Fit: 71.72% IDQR Fit: 71.72% IDAUDI Fit: 72.01% IDL1 Fit: 71.64%

10 8

Magnitud

6 4 2 0 -2 -4 -6 0

100

200

300

400 500 600 Tiempo (s)

700

800

900

1000

e) Comparación de las salidas de los modelos estimados. Figura 6.11 Estimaciones paramétricas del proceso banco de nivel RCN100 para un modelo tipo ARX.

6.3.2.2 Estructura ARMAX En la Tabla 6.11 se presentan los índices de desempeño ECMA y EAM de los modelos estimados para las respuestas al escalón e impulso, para una estructura ARMAX, con su respectivo orden.

ECM EAM

Tabla 6.11 ECM y EAM de los modelos estimados para el proceso banco de nivel con estructura ARMAX. ARMAX [2 1 1 1] LS [2 1 1 1] LU [2 1 1 1] QR [2 1 1 1] AUDI [2 2 1 1] L1 [2 1 1 1] Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. 0.2800 0.0025 0.2407 0.0022 0.2800 0.0025 0.2800 0.0025 0.6671 0.0057 0.2779 0.0026 0.0322 0.0003 0.0278 0.0003 0.0322 0.0003 0.0322 0.0003 0.0882 0.0006 0.0320 0.0003

La Tabla 6.13 muestra el tiempo TE para la familia de algoritmos LU y para los algoritmos ARMAX, LS, AUDI y L1. Se puede observar que para este caso de estudio, el que presentó mayor tiempo de ejecución es el modelo AUDI, siendo el modelo de QR, en todas sus variantes, el que menos tiempo de ejecución tomó para el cálculo de los parámetros. Tabla 6.12 Tiempo de ejecución de ARMAX, LS, LU, AUDI y L1 para el banco de nivel RCN100 con estructura ARMAX. ARMAX LS LU1 LU2 LU3 LU4 AUDI L1 TE 1.3020s 0.1310s 0.1400 0.0400 0.0400 0.0400 1.5720s 0.4500s

En la Tabla 6.13 se presentan los tiempo TE para la familia de algoritmos QR, y los algoritmos ARMAX, LS, AUDI y L1. Tabla 6.13 Tiempo de ejecución de ARMAX, LS, QR, AUDI y L1 para el banco de nivel RCN100 con estructura ARMAX. ARMAX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 1.3020s 0.1310s 0.1000s 0.1100s 0.1110s 1.3220s 1.3620s 1.5720s 0.4500s

107

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Las estimaciones paramétricas correspondientes a los modelos para una estructura ARMAX se muestran en la Figura 6.12. Para los modelos estimados con estructura ARMAX el modelo menos representativo de la dinámica del proceso real es el que se obtiene con el algoritmo AUDI. Respuesta al impulso de los modelos obtenidos

Respuesta al escalon de los modelos obtenidos 20

REAL

18

0.25

LS ARMAX

16

LU

0.2

14

QR AUDI

Amplitud

10 8

REAL

L1

0.15

0.1

LS

6 4

ARMAX LU QR

2

AUDI

0.05

0

L1 0 50

100

150 Tiempo (s)

200

250

50

a) Respuesta al escalón.

100

150 Tiempo (s)

b) Respuesta al impulso.

BODE

2

Amplitud

10

1

10

0

10

-1

10

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0 REAL Fase (grados)

Amplitud

12

-50

LS ARMAX LU

-100

QR

-150 -200 0

0.01

0.02

0.03

0.04 0.05 0.06 Frecuencia (Hz)

0.07

0.08

c) Respuesta en frecuencia 1.

108

0.09

0.1

200

250

Capítulo 6: Validación de algoritmos

BODE

2

Amplitud

10

1

10

0

10

-1

10

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0 Fase (grados)

REAL -50

AUDI L1

-100 -150 -200 0

0.01

0.02

0.03

0.04 0.05 0.06 Frecuencia (Hz)

0.07

0.08

0.09

0.1

d) Respuesta en frecuencia 2. Salida del sistema real y salida de los modelos estimados 12 Measured Output IDLS Fit: 61.3%

10

Magnitud

IDARMAX Fit: 48.32% 8

IDLU Fit: 61.3% IDQR Fit: 61.3%

6

IDAUDI Fit: 30.5% IDL1 Fit: 47.11%

4 2 0 -2 -4 -6 0

200

400 600 Tiempo (s)

800

1000

e) Comparación de las salidas de los modelos estimados. Figura 6.12 Estimaciones paramétricas del proceso banco de nivel RCN100 para un modelo tipo ARMAX.

6.3.3

Secadora de pelo (Dryer)

El tercer caso de estudio es un experimento de laboratorio que emula el funcionamiento de una secadora de pelo (ver descripción del proceso en el Apéndice A.3). En este proceso la variable de interés es la temperatura del aire a la salida del tubo. Por lo tanto su dinámica de es lenta presentando un retardo durante el transitorio. 6.3.3.1 Estructura ARX Los datos de entrada/salida se muestran en la Figura 6.13.

109

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. Datos de Entrada/Salida

Magnitud

6

5

4 0

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

Magnitud

6

5

4

0

Figura 6.13 Datos de entrada/salida del proceso secadora de pelo.

Después de eliminar tendencias a los datos de entrada/salida, los nuevos vectores de entrada/salida que se utilizan para el proceso de identificación se muestran en la siguiente figura. Datos de Entrada/Salida

Magnitud

1

0

-1 0

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

100

200

300

400 500 600 Numero de muestras

700

800

900

1000

Magnitud

1

0

-1 0

Figura 6.14 Datos de entrada/salida después de eliminar tendencias del proceso secadora de pelo.

Las estimaciones no paramétricas para este caso de estudio se muestran en la Figura 6.15 y la Figura 6.16.

110

Capítulo 6: Validación de algoritmos

Respuesta al escalon del modelo no paramétrico (REAL)

REAL

Amplitud

Amplitud

0.5

0 0

-1

10

-2

10

-3

0.5

1

1.5 Tiempo (s)

2

10

2.5

Respuesta al impulso del modelo no parampetrico (REAL)

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0 Fase (grados)

1.5

Amplitud

BODE

0

10

1

1 0.5

REAL -200

-400

0 0.5

1

1.5 Tiempo (s)

2

2.5

Figura 6.15 Respuesta al escalón y al impulso de la secadora de pelo.

-600 0

0.5

1

1.5

2

2.5 3 Frecuencia (Hz)

3.5

4

4.5

5

Figura 6.16 Respuesta en frecuencia de la secadora de pelo.

Los índices de desempeño para los modelos estimados con estructura ARX para cada algoritmo se muestran en la Tabla 6.14.

ECM EAM

Tabla 6.14 ECM y EAM de modelos estimados para el proceso secadora de pelo con estructura ARX. ARX [6 5 1] LS [6 5 1] LU [6 5 1] QR [6 5 1] AUDI [6 6 1] L1 [6 5 1] Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. 0.0028 0.0021 0.0028 0.0021 0.0028 0.0021 0.0028 0.0021 0.0028 0.0021 0.0068 0.0553 0.0005 0.0003 0.0005 0.0003 0.0005 0.0003 0.0005 0.0003 0.0005 0.0003 0.0011 0.0076

En la siguiente tabla se presentan los tiempos de ejecución para la familia de algoritmos LU, y los algoritmos ARX, LS, AUDI y L1. El algoritmo L1 es el que tiene mayor tiempo de ejecución. Tabla 6.15 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para la secadora de pelo con estructura ARX.

ARX LS LU1 LU2 LU3 LU4 AUDI L1 TE 0.4710s 0.0700s 0.1200s 0.0900s 0.0600s 0.0800s 0.0300s 1.1420s Para la familia de algoritmos QR, y los algoritmos ARX, LS, AUDI y L1 en la Tabla 6.16 se muestran los tiempos de ejecución. En este caso los algoritmos que consumieron mayor tiempo de ejecución son QR4, QR5 y L1. Tabla 6.16 Tiempo de ejecución de ARX, QR, LU, AUDI y L1 para la secadora de pelo con estructura ARX.

ARX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 0.4710s 0.0700s 0.0700s 0.0800s 0.0900s 1.8120s 1.6220s 0.0300s 1.1420s En la Figura 6.17 se presentan las estimaciones paramétricas correspondientes a los modelos estimados. El algoritmo menos representativo a la dinámica del proceso que se está identificando es el que se obtiene con L1. Los algoritmos restantes presentan una buena aproximación.

111

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. Respuesta al escalon de los modelos obtenidos

Respuesta al impulso de los modelos obtenidos

1 REAL LS

2 0.8

ARX LU QR

1.5

AUDI

REAL LS

1

0.5

ARX

0.2

L1

Amplitud

0.4

LU QR 0

AUDI 0

L1 0.5

1

1.5 Tiempo (s)

2

2.5

0.5

a) Respuesta al escalón.

1

1.5 Tiempo (s)

b) Respuesta al impulso. BODE

0

Amplitud

10

-1

10

-2

10

-3

10

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Fase (grados)

0 REAL LS ARX LU QR

-200

-400

-600 0

0.5

1

1.5

2

2.5 3 Frecuencia (Hz)

3.5

4

4.5

5

4.5

5

c) Respuesta en frecuencia 1.

Amplitud

10

10

10

10

BODE

0

-1

-2

-3

0

0.5

1

1.5

2

2.5

3

3.5

4

0 REAL Fase (grados)

Amplitud

0.6

AUDI

-200

L1

-400

-600 0

0.5

1

1.5

2

2.5 3 Frecuencia (Hz)

3.5

4

d) Respuesta en frecuencia 2.

112

4.5

5

2

2.5

Capítulo 6: Validación de algoritmos

Salida del sistema real y salida de los modelos estimados 1.5 Measured Output IDLS Fit: 88.78% IDARX Fit: 88.78%

1

IDLU Fit: 88.78% IDQR Fit: 88.78% IDAUDI Fit: 88.78%

0.5

Magnitud

IDL1 Fit: 80.46% 0

-0.5

-1

-1.5

-2 0

10

20

30

40

50 60 Tiempo (s)

70

80

90

100

d) Comparación de las salidas de los modelos estimados. Figura 6.17 Estimaciones paramétricas del proceso secadora de pelo para un modelo tipo ARX.

6.3.3.2 Estructura ARMAX Para la estructura ARMAX la desviación del error de los modelos estimados se muestra en la Tabla 6.17.

ECM EAM

Tabla 6.17 Desviación del error en los modelos estimados para el proceso secadora de pelo con estructura ARMAX. ARMAX [6 5 1 1] LS [6 5 1 1] LU [6 5 1 1] QR [6 5 1 1] AUDI [3 3 2 3] L1 [10 7 3 2] Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. 0.0026 0.0021 0.0064 0.0045 0.0026 0.0021 0.0026 0.0021 0.0334 0.0322 0.0029 0.0239 0.0005 0.0003 0.0011 0.0009 0.0005 0.0003 0.0005 0.0003 0.0060 0.0043 0.0005 0.0040

En la Tabla 6.19 se concentran los tiempos de ejecución para cada algoritmo. Se puede observar que hasta el momento los algoritmos que tienen un menor tiempo de ejecución son LU y QR en sus 4 métodos distintos. A pesar de que L1 es una función de optimización que requiere muchas iteraciones, la diferencia de ejecución es menor a un segundo. Tabla 6.18 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para la secadora de pelo con estructura ARX. ARMAX LS LU QR1 QR2 QR3 AUDI L1 TE 1.4220s 0.1910 0.1700s 0.0600s 0.0900s 0.0900s 2.2940s 2.2230s

En la Tabla 6.19 se presentan los tiempos de ejecución para la familia de algoritmos QR y los algoritmos ARMAX, LS AUDI y L1. Tabla 6.19 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para la secadora de pelo con estructura ARX. ARMAX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 1.4220s 0.1910 0.0600s 0.0900s 0.0900s 0.0700s 0.0800s 2.2940s 2.2230s

La Figura 6.18 presenta las estimaciones no paramétricas para la estructura ARMAX del proceso secadora de pelo. Para este caso el mejor modelo se obtiene con LU, QR y ARMAX (toolbox de identificación), y el peor modelo es el que se obtiene con L1.

113

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Respuesta al impulso de los modelos obtenidos

Respuesta al escalon de los modelos obtenidos 1.5

1

REAL LS ARMAX

0.8

LU QR

1

AUDI L1

Amplitud

0.4 REAL

0.5

LS ARMAX

0.2

LU QR AUDI

0

0

L1 0.5

1 Tiempo (s)

1.5

0

2

0.5

Amplitud

10

10

10

10

1

1.5 Tiempo (s)

a) Respuesta al escalón.

b) Respuesta al impulso. BODE

0

-1

-2

-3

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Fase (grados)

0 REAL LS ARMAX LU QR

-200

-400

-600 0

0.5

1

1.5

2

2.5 3 Frecuencia (Hz)

3.5

4

4.5

5

c) Respuesta en frecuencia 1. BODE

0

10

Amplitud

0

-1

10

-2

10

-3

10

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0 REAL Fase (grados)

Amplitud

0.6

-200

AUDI L1

-400 -600 -800 0

0.5

1

1.5

2

2.5 3 Frecuencia (Hz)

3.5

4

d) Respuesta en frecuencia 2.

114

4.5

5

2

Capítulo 6: Validación de algoritmos

Salida del sistema real y salida de los modelos estimados 1.5 Measured Output IDLS Fit: 88.78% IDARMAX Fit: 88.44%

1

IDLU Fit: 88.78% IDQR Fit: 88.78% IDAUDI Fit: 73.09%

0.5

Magnitud

IDL1 Fit: 87.28% 0

-0.5

-1

-1.5

-2 0

10

20

30

40 50 60 Tiempo (s)

70

80

90

100

e) Comparación de las salidas de los modelos estimados. Figura 6.18 Estimaciones paramétricas del proceso secadora de pelo para un modelo tipo ARMAX

6.3.4

Sistema de calentamiento (Heating System)

El último caso de estudio es un sistema de calentamiento cuya descripción se presenta en el Apéndice A.4. En este proceso también la variable de interés es la temperatura, por lo que presenta un retardo durante el transitorio al igual que el proceso anterior. El experimento originalmente se llevó a cabo para el diseño de un controlador robusto. Por las características propias de diseño, la señal de excitación no es una señal de excitación persistente. 6.3.4.1 Estructura ARX En la Figura 6.19 se muestran as gráficas correspondientes a los datos de entrada/salida del sistema a identificar. Datos de Entrada/Salida

Magnitud

8

6

4 0

100

200

300 400 500 Numero de muestras

600

700

800

100

200

300 400 500 Numero de muestras

600

700

800

Magnitud

180 160 140 120 100 80 0

Figura 6.19 Datos de entrada-salida del sistema de calentamiento.

115

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Una vez eliminadas las tendencias de los datos entrada/salida en la siguiente figura se muestran los nuevos datos con los que se llevará a cabo la identificación de los modelos. Datos de Entrada/Salida

Magnitud

2

0

-2 0

100

200

300 400 500 Numero de muestras

600

700

800

100

200

300 400 500 Numero de muestras

600

700

800

Magnitud

50

0

-50 0

Figura 6.20 Datos de entrada/salida después de eliminar tendencias del proceso sistema de calentamiento.

Las estimaciones no paramétricas se muestran en la Figura 6.21 y la Figura 6.22. Para obtener la respuesta en frecuencia de este proceso es necesario aumentar la resolución de frecuencias γ mediante el comando spa(MODELO, γ ). Donde γ = 100 para este caso. La selección del tamaño de γ es un compromiso entre resolución de frecuencia y varianza (variabilidad) [30]. Un buen valor de γ para sistemas sin resonancias muy pronunciadas es de γ = 20 − 30 . Para un espectro con picos de resonancia muy estrechos es necesario seleccionar un valor de γ elevado. El tamaño de ventana por defecto se calcula mediante min(length(DATA)/10,30), y se puede modificar usando SPA(DATA,M). Donde M = γ . Respuesta al escalon del modelo no paramétrico (REAL) 10

BODE

2

10 5

REAL

0

Amplitud

Amplitud

15

-5 10

20

30

40

50 60 Tiempo (s)

70

80

90

10

100

Respuesta al impulso del modelo no parampetrico (REAL)

0

-2

0

0.05

0.1

0.15

0.2

0.25

0.05

0.1 0.15 Frecuencia (Hz)

0.2

0.25

0 Fase (grados)

0.4 Amplitud

10

0.2 0

-200 -400 -600 REAL

0

10

20

30

40

50 60 Tiempo (s)

70

80

90

-800 0

100

Figura 6.21 Respuesta al escalón y al impulso de la secadora de pelo.

Figura 6.22 Respuesta en frecuencia de la secadora de pelo.

116

Capítulo 6: Validación de algoritmos

Para un modelo tipo ARX los índices de desempeño ECM y EAM de los modelos estimados se muestran en la Tabla 6.20. También es importante remarcar que con el algoritmo L1 se estimó un modelo de mucho menor orden que los restantes algoritmos.

Tabla 6.20 ECM y EAM de los modelos estimados para el proceso sistema de calentamiento con estructura ARX. ARX [8 7 1] LS [8 7 1] LU [8 7 1] QR [8 7 1] AUDI [8 8 1] L1 [1 2 2] Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. ECM 0.0201 0.0019 0.0201 0.0019 0.0201 0.0019 0.0201 0.0019 0.0269 0.0018 0.1245 0.0061 EAM 0.0020 0.0002 0.0020 0.0002 0.0020 0.0002 0.0020 0.0002 0.0035 0.0002 0.0147 0.0007

En la Tabla 6.21 se presentan los tiempos de ejecución para la familia de algoritmos LU, y los algoritmos ARX, LS, AUDI y L1, predominando con un menor tiempo los algoritmos alternativos. Tabla 6.21 Tiempo de ejecución de ARX, LS, LU, AUDI y L1 para el sistema de calentamiento con estructura ARX. ARX LS LU1 LU2 LU3 LU4 AUDI L1 TE 0.4810s 0.0700s 0.0800s 0.0800s 0.0800s 0.0800s 0.0300s 0.2100s

Para la familia de algoritmos QR, y los algoritmos ARX, LS, AUDI y L1 el tiempo de ejecución se presentan en la siguiente tabla. Tabla 6.22 Tiempo de ejecución de ARX, LS, QR, AUDI y L1 para el sistema de calentamiento con estructura ARX. ARX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 0.4810s 0.0700s 0.0600s 0.0800s 0.0700s 1.4920s 1.1820s 0.0300s 0.2100s

Los modelos estimados con los algoritmos QR4 y QR5 consumen más tiempo de ejecución y como se mencionó anteriormente, se debe a que el número de parámetros a estimar son considerables. Esto afecta directamente a los algoritmos de Givens y Householder, QR4 y QR5 respectivamente; ya que para encontrar la factorización QR se basa en el cálculo de matrices de forma sucesiva hasta que la matriz final Q sea ortogonal y la matriz R sea triangular superior. En la Figura 6.23 se muestran las gráficas correspondientes a la estimación paramétrica. En el análisis temporal el único algoritmo que presenta una notoria diferencia entre el modelo estimado y el modelo real es el algoritmo L1 (Figura 6.23 a) y b)). En la respuesta en frecuencia el algoritmo que tiene un mejor desempeño es L1, ya que debido a su propiedad de robustez mantiene una trayectoria promedia a pesar de las variaciones que se presentan en el rango de frecuencias de 0.1-0.25 Hz.

117

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. Respuesta al escalon de los modelos obtenidos

Respuesta al impulso de los modelos obtenidos 0.5

16

REAL LS

14

ARX 0.4

LU QR AUDI

Amplitud

10 8 6

REAL LS

4

ARX

2

LU QR

0.3

L1

0.2

0.1

AUDI 0

L1

0 20

30

40

50 60 Tiempo (s)

70

80

90

100

0

10

a) Respuesta al escalón.

20

30

40

50 60 Tiempo (s)

BODE

2

Amplitud

0

10

-2

10

0

0.05

0.1

0.15

0.2

0.25

0.05

0.1 0.15 Frecuencia (Hz)

0.2

0.25

0 -200 REAL LS ARX LU

-400 -600

QR -800 0

c) Respuesta en frecuencia 1. 10

10

10

10

BODE

2

0

-2

-4

0

0.05

0.1

0.15

0.2

0.25

0.05

0.1 0.15 Frecuencia (Hz)

0.2

0.25

0

-500

REAL AUDI L1

-1000 0

70

80

b) Respuesta al impulso.

10

Fase (grados)

10

Amplitud

0

Fase (grados)

Amplitud

12

d) Respuesta en frecuencia 2.

118

90

100

Capítulo 6: Validación de algoritmos

Salida del sistema real y salida de los modelos estimados 60 Measured Output IDLS Fit: 89% IDARX Fit: 88.99%

40

IDLU Fit: 89% IDQR Fit: 89% IDAUDI Fit: 89.69%

20

Magnitud

IDL1 Fit: 82.13% 0

-20

-40

-60

-80 0

200

400

600

800 1000 Tiempo (s)

1200

1400

1600

1800

e) Comparación de las salidas de los modelos estimados. Figura 6.23 Estimaciones paramétricas del proceso sistema de calentamiento para un modelo tipo ARX.

6.3.4.2 Estructura ARMAX Para el caso de la estimación de parámetros para una estructura tipo ARMAX, los coeficientes ECM y EAM se resumen en la Tabla 6.23.

ECM EAM

Tabla 6.23 Desviación del error en los modelos estimados para el proceso sistema de calentamiento con estructura ARMAX. ARMAX [8 7 3 1] LS [8 7 3 1] LU [8 7 3 1] QR [8 7 3 1] AUDI [4 4 2 2] L1 [3 3 1 1] Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. Esc. Imp. 0.0509 0.0036 0.0575 0.0024 0.0509 0.0036 0.0509 0.0036 0.0581 0.0091 0.0278 0.0180 0.0086 0.0006 0.0103 0.0004 0.0086 0.0006 0.0086 0.0006 0.0079 0.0011 0.0493 0.0026

La Tabla 6.24 muestra el tiempo de ejecución para la familia de algoritmos LU y para los algoritmos ARMAX, LS, AUDI y L1. Se sigue manteniendo el patrón de que los algoritmos alternativos consumen menor tiempo de ejecución que la función ARMAX del Toolbox de Identificación. Tabla 6.24 Tiempo de ejecución de ARMAX, LS, LU, AUDI y L1 para el sistema de calentamiento con estructura ARMAX. ARMAX LS LU1 LU2 LU3 LU4 AUDI L1 TE 1.5620s 0.1700s 0.1610s 0.1300s 0.1510s 0.1500s 1.3320s 0.3800s

En la Tabla 6.25 se presentan los tiempos de ejecución para la familia de algoritmos QR, y los algoritmos ARMAX, LS, AUDI y L1. Para este caso el algoritmo que consume mayor tiempo de ejecución en el cálculo de los parámetros es QR4, y el que consume menor tiempo es L1, además de estimar un modelo de mucho menor orden que los otros algoritmos. Tabla 6.25 Tiempo de ejecución de ARMAX, LS, QR, AUDI y L1 para el sistema de calentamiento con estructura ARMAX. ARMAX LS QR1 QR2 QR3 QR4 QR5 AUDI L1 TE 1.5620s 0.1700s 0.1300s 0.1510s 0.1500s 1.7620s 1.3920s 1.3320s 0.3800s

119

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Para este caso a pesar de que L1 presenta una deferencia importante en el análisis transitorio, en la respuesta en frecuencia sigue una trayectoria aceptable con respecto al sistema real. Respuesta al escalon de los modelos obtenidos

Respuesta al impulso de los modelos obtenidos

16

0.5 REAL LS

0.45

14

0.4

ARMAX

12

LU QR

0.35

AUDI

Amplitud

0.3

8 REAL LS

6

L1 0.25 0.2 0.15

ARMAX

4

0.1

LU QR

2

0.05

AUDI L1

0 20

30

40

50 60 Tiempo (s)

70

80

90

0 100

0

a) Respuesta al escalón.

10

20

30

40

50 60 Tiempo (s)

BODE

2

0

10

-2

10

-4

10

0

0.05

0.1

0.15

0.2

0.25

0.1 0.15 Frecuencia (Hz)

0.2

0.25

0 -200 -400

REAL LS

-600

ARMAX LU QR

-800 0

0.05

70

b) Respuesta al impulso.

10

Amplitud

10

Fase (grados)

Amplitud

10

c) Respuesta en frecuencia 1.

120

80

90

100

Capítulo 6: Validación de algoritmos

BODE

2

Amplitud

10

0

10

-2

10

0

0.05

0.1

0.15

0.2

0.25

0.1 0.15 Frecuencia (Hz)

0.2

0.25

Fase (grados)

0 -200 -400 REAL AUDI L1

-600 -800 0

0.05

d) Respuesta en frecuencia 2. Salida del sistema real y salida de los modelos estimados 80 Measured Output IDLS Fit: 84.46% IDARMAX Fit: 91.63% IDLU Fit: 84.46% IDQR Fit: 84.46% IDAUDI Fit: 89.59% IDL1 Fit: 84.05%

60 40

Magnitud

20 0 -20 -40 -60 -80 0

200

400

600

800 1000 Tiempo (s)

1200

1400

1600

1800

e) Comparación de las salidas de los modelos estimados. Figura 6.24 Estimaciones paramétricas del proceso sistema de calentamiento para un modelo tipo ARMAX

6.4 Conclusiones Este capítulo se enfocó a la identificación de los cuatro procesos seleccionados: a) Sistema de tercer orden, b) banco de nivel RCN100, c) secadora de pelo (Dryer) y d) sistema de calentamiento (heating system). Se abordaron dos tipos de estructuras de las existentes para estimación paramétrica, ARX y ARMAX.

Con respecto al ajuste de la dinámica de los modelos estimados. Para el primer caso de estudio, el cual es un sistema muy simple que su dinámica esencial es la magnitud del sobre pico, y de acuerdo con las gráficas que se mostraron en la Sección

121

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

6.3.1 para ambas estructuras, se tiene que los modelos estimados se consideran buenos debido a que caracterizan lo suficientemente bien la dinámica del sistema. Para el segundo caso de estudio, el banco de nivel RCN100, es un sistema cuya respuesta transitoria ante una entrada de tipo escalón, su dinámica se asemeja a un sistema de primer orden. Los modelos que se obtuvieron para este proceso en su análisis temporal, para la estructura tipo ARX, tiene un mejor ajuste que los modelos obtenidos de estructura tipo ARMAX. Y si observamos la respuesta en frecuencia de ambas estructuras, el modelo que no es tan aceptable es el que se obtiene con el algoritmo AUDI para un modelo con estructura de tipo ARMAX. Los siguientes dos procesos son procesos de dinámica lenta debido a que la variable de interés en ambos es la temperatura, por lo tanto presentan un retardo en su respuesta transitoria. Para el proceso: secadora de pelo, el modelo que no es lo suficientemente representativo de la dinámica del proceso real es el modelo que se obtuvo con la norma L1, esto tanto para la estructura ARX y ARMAX. Por otro lado los modelos estimados con los otros algoritmos se consideran buenos modelos. Por último para el proceso: sistema de calentamiento, el modelo estimado con mejor ajuste a la dinámica real es el que se obtuvo con la norma L1 ya que debido a su propiedad de robustez presenta una trayectoria promedio en el rango de frecuencias (en la respuesta en frecuencia); mientras que los modelos restantes tienden a seguir los cambios que se presentan específicamente en alta frecuencia. En la Tabla 6.26 se resumen los órdenes óptimos que se obtuvieron para la familia de algoritmos LU, LU clásico, Cholesky, Crout y Doolittle que se aplicaron a identificación para una estructura tipo ARX. Tabla 6.26 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos LU para una estructura tipo ARX.

LU1 LU2

LU3

LU4

Tercer Orden

3

3

3

3

RCN100

2

2

2

2

Dryer

6

6

6

6

122

Capítulo 6: Validación de algoritmos

Heating System

8

8

8

8

Los órdenes óptimos para los modelos estimados con estructura ARMAX se presentan en la siguiente tabla. Tabla 6.27 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos LU para una estructura tipo ARMAX.

LU1 LU2

LU3

LU4

Tercer Orden

3

3

3

3

RCN100

2

2

2

2

Dryer

6

6

6

6

Heating System

8

8

8

8

NOTA: Modelado:

bueno,

aceptable,

regular

Los órdenes óptimos para la familia de algoritmos QR, función QR de MATLAB®, GramSchmidt clásico, Gram-Schmidt modificado, rotaciones de Givens y reflexiones de Householder, que se aplicaron a identificación para una estructura tipo ARX se muestran en la Tabla 6.28 Tabla 6.28 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos QR para una estructura tipo ARX.

QR1 QR2

QR3

QR4

QR5

Tercer Orden

3

3

3

3

3

RCN100

2

2

2

2

2

Dryer

6

6

6

6

6

Heating System

8

8

8

8

8

Los órdenes correspondientes de la familia de algoritmos QR para una estructura tipo ARMAX se muestran en la Tabla 6.29. Tabla 6.29 Ordenes óptimos para los algoritmos ARX, LS, AUDI, L1 y la familia de algoritmos QR para una estructura tipo ARMAX.

123

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

QR1 QR2

QR3

QR4

QR5

Tercer Orden

3

3

3

3

3

RCN100

2

2

2

2

2

Dryer

6

6

6

6

6

Heating System

8

8

8

8

8

NOTA: Modelado:

bueno,

aceptable,

regular

Es importante mencionar que cada una de las diferentes formas de obtener la factorización LU y QR se utilizaron para identificar cada caso de estudio obteniendo buenos modelos a excepción del último caso de estudio, donde estas familias obtienen modelos de orden elevado mientras que el algoritmo L1 obtiene un orden menor.

Con respecto al tiempo de cálculo transcurrido en la obtención de los parámetros. El tiempo de cálculo en términos generales en la mayoría de los casos fue menor que el tiempo de cálculo que consumen las funciones que contiene el Toolbox de Identificación de MATLAB®. Es importante resaltar que las funciones que contiene MATLAB® para realizar las factorizaciones son rutinas pre-compiladas, es decir, el tiempo de ejecución es más rápido que un archivo que se genera en el editor de texto. Mientras que las rutinas programadas son programas interpretados que se generan en el editor m-file de MATLAB®. Esta reducción de tiempo es benéfica en tareas de identificación durante el proceso de obtención del modelo paramétrico. Para que el tiempo de ejecución fuera los más confiable posible se tuvieron que desactivar los servicios de Microsoft que se ejecutan en segundo plano así como programas de monitoreo. Esto con la finalidad de obtener un grado de precisión lo más cercano al tiempo real de ejecución. Los programas codificados se ejecutaron en una computadora HP Pavilion ze5427LA con un procesador Intel Pentium(R) 4 de 2.40 GHz y 512 MB de memoria RAM.

124

Capítulo 6: Validación de algoritmos

Es posible que si se ejecutan los programas codificados para la estimación de parámetros en otra computadora de mejores características, los tiempos de ejecución reduzcan pero el umbral de diferencia en los tiempos se mantendrá. Una alternativa para medir el tiempo de ejecución es utilizando un lenguaje de bajo nivel, como ensamblador, para tener una medición exacta de los tiempos que ocupa cada operación e instrucción.

125

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

126

Capítulo 7 Conclusiones y trabajos futuros 7.1 Conclusiones El algoritmo de mínimos cuadrados ha sido una de las técnicas de estimación de parámetros más utilizadas durante muchos años. La aplicación de ésta técnica en la estimación de parámetros se han estudiado demasiado, siendo coincidentes todos los trabajos en que la matriz A puede estar mal condicionada cuando es de gran dimensión, que es el caso de la identificación de sistemas, con lo cual genera errores numéricos importantes en la solución de (1.5). Por tal motivo se presentaron diferentes alternativas para resolver esta ecuación, dentro de las cuales se pueden nombrar a: factorización LU, factorización QR, factorización AUDI y la norma L1 como criterio de minimización del error entre la salida estimada y la salida real. En la Tabla 7.1 se resumen algunos de los criterios que se utilizan para la selección y validación de los modelos. Para cada caso de estudio los criterios que se utilizaron están marcados con “X”.

127

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Tabla 7.1 Criterios para la selección y validación de modelos.

PROCESOS

SELECCIÓN VALIDACION

COMPARACIÓN APD

FPE L1

X

X

RCN100

X

X

X

X

X

X

DRYER

X

X

X

X

X

X

HS

X

X

X

X

X

X

AP: Validación en base a la aplicación del modelo. CPF: Comprobación de parámetros fijos. CIO: Coherencia en el comportamiento entrada-salida. RM: Reducción de modelo. IFP: Intervalos de fiabilidad de parámetros. S: Simulación

APD: Análisis preliminares de los datos. AE: Análisis espectral. RMC: Evaluar rango de la matriz de covarianza. CO: Correlación. LFCN: Función de perdida (Loss Function). FPE: Error de predicción final (Final predective error). CP-Cs: Cancelación polos-ceros. PR: Prueba de residuo.

NOTA: Para el caso del algoritmo utilizando la norma L1 se generó la función FPE para dicha norma. Los criterios utilizados para la selección y validación de los modelos estimados con cada uno de los algoritmos se encuentran marcados con una X en la Tabla 7.1. Para el primer caso de estudio se utilizó el criterio de análisis espectral para la etapa de selección durante el paso de análisis preliminares de los datos. Esto porque se generó la señal de excitación cumpliendo con base a la Definición 1.1. Otro criterio durante la etapa de selección fue el criterio de información de Akaike, por consiguiente se empleó el factor de predicción de error para L2, que se describió en la Sección 5.4. Para la etapa de validación del modelo se utilizó el criterio de comprobación de parámetros fijos ya que se conocen los valores de los parámetros a estimar. 128

PR

FPE L2

X

C P-Cs

LFCN

X

CO

X

RMC

AE

X

IFP

X

Tercer Orden

RM

CIO

S

AP

CPF

F-TEST

FCN. ERROR

Capítulo 7: Conclusiones y trabajos futuros

Además del criterio de simulación donde se comparan las salidas tanto del modelo no paramétrico y el modelo estimado ante señales de entrada distinta con respecto a la original. Para los procesos banco de nivel, secadora de pelo y sistema de calentamiento, en la etapa de selección de modelos se utilizaron los criterios: Análisis espectral y el criterio de información de Akaike (5.4) y el error de predicción final para L1 (5.5). Para la etapa de validación se utilizó el criterio de coherencia en el comportamiento de entrada-salida mediante el apoyo de las gráficas de BODE. En la Tabla 7.2 se muestran los órdenes de los modelos estimados de cada caso de estudio, para una estructura ARX. Para la mayoría de los procesos todos los algoritmos presentan buen ajuste con respecto al sistema real. El modelo que se obtuvo con L1 para la secadora de pelo no es un modelo que represente la dinámica del proceso real. Sin embargo para el último caso de estudio con el algoritmo L1 se obtuvo un modelo de menor orden que los algoritmos restantes. Tabla 7.2 Orden máximo de cada algoritmo para cada proceso con una estructura ARX.

ARX

LS

LU

QR

AUDI

L1

Tercer Orden

3

3

3

3

3

3

RCN100

2

2

2

2

3

2

Dryer

6

6

6

6

7

6

8

8

8

8

9

4

Heating System

Para el caso de los modelos estimados con estructura ARMAX la Tabla 7.3 muestra los órdenes óptimos. En general se obtuvieron bueno modelos con los algoritmos alternativos; al igual que en el caso de una estructura de tipo ARX, para el caso de estudio del proceso secadora de pelo, el peor modelo se obtuvo con el algoritmo L1. También para la estructura ARMAX el modelo estimado con el algoritmo L1 es un buen modelo y de un orden menor que los algoritmos restantes. La información que arrojan estas dos tablas es importante, ya que nos dice que los algoritmos alternativos si son capaces de reproducir la dinámica de procesos reales y en un tiempo mucho menor que las funciones tradicionales que contiene el Toolbox de

129

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

identificación como se muestra en la Figura 7.1 para la estructura tipo ARX y la Figura 7.3 para la estructura tipo ARMAX. Tabla 7.3 Orden máximo de cada algoritmo para cada proceso con una estructura ARMAX.

ARMAX

LS

LU

Tercer Orden

3

3

3

3

3

3

RCN100

2

2

2

2

3

2

Dryer

6

6

6

6

6

10

8

8

8

8

6

4

Heating System

NOTA: Modelado:

bueno,

QR AUDI

aceptable,

L1

regular

En la Figura 7.1 se muestran los tiempos de ejecución de la familia de algoritmos LU y los algoritmos ARX, LS, AUDI y L1 para cada uno de los casos de estudio teniendo como estructura un modelo tipo ARX. Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARX 2 1.8

Tiempo de Ejecución (s)

1.6 1.4 1.2 1

ARX LS LU CHOL CROUT DOOLITTLE AUDI L1

0.8 0.6 0.4 0.2 0 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.1 Tiempos de ejecución de los algoritmos ARX, LU, LS, AUDI y L1 para una estructura ARX.

Se puede observar que los algoritmos que consumen mayor tiempo de ejecución son la función arx de MATLAB® y el algoritmo L1. En la Figura 7.3 se hace un acercamiento para observar de manera más clara la diferencia entre los algoritmos LS, LU y AUDI. Es claro que la familia de algoritmos LU tiene una diferencia en el tiempo de ejecución relativamente pequeña, siendo el algoritmo AUDI el que consume menor tiempo para todos los procesos.

130

Capítulo 7: Conclusiones y trabajos futuros

Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARX ARX LS LU CHOL CROUT DOOLITTLE AUDI L1

Tiempo de Ejecución (s)

0.15

0.1

0.05

0

-0.05

Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.2 Acercamiento de los tiempos de ejecución de los algoritmos ARX, LU, LS, AUDI y L1 para una estructura ARX.

Para una estructura tipo ARMAX la Figura 7.3 muestra los tiempos de ejecución de la familia de algoritmos LU y los algoritmos ARMAX, LS, AUDI y L1 para cada uno de los casos de estudio. En este caso el algoritmo que consume mayor tiempo de ejecución son: AUDI, ARMAX y

Tiempo de Ejecución (s)

L1. Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARMAX 2.5 ARMAX LS LU CHOL 2 CROUT DOOLITTLE AUDI L1 1.5

1

0.5

0 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.3 Tiempos de ejecución de los algoritmos ARMAX, LS, LU, AUDI y L1 para una estructura ARMAX.

En la Figura 7.4, haciendo un acercamiento, se puede observar que la familia de algoritmos LU es la que utiliza menor tiempo de ejecución durante la estimación de los parámetros.

131

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARMAX 0.3 ARMAX LS 0.25 LU

Tiempo de Ejecución (s)

0.2

0.15

CHOL CROUT DOOLITTLE AUDI L1

0.1

0.05

0

-0.05 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.4 Acercamiento de los tiempos de ejecución de los algoritmos ARMAX, LS, LU, AUDI y L1 para una estructura ARMAX.

Ahora bien, la Figura 7.5 muestra los tiempos de ejecución par la familia de algoritmos QR,

Tiempo de Ejecución (s)

y los algoritmos ARX, LS, AUDI y L1. Donde los algoritmos que consumen mayor tiempo de ejecución son QR mediante rotaciones de Givens y reflexiones Householder. Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARX 2 ARX LS 1.8 QR GMC GM 1.6 M GIVENS HOUSEHOLDER 1.4 AUDI L1 1.2 1 0.8 0.6 0.4 0.2 0 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.5 Tiempos de ejecución de los algoritmos ARX, LS, QR, AUDI y L1 para una estructura ARX.

Haciendo un acercamiento en la Figura 7.6 se observa que el algoritmo AUDI es el que mantiene un tiempo de ejecución constante en cada uno de los procesos, no obstante los algoritmos Gram-Schmidth clásico, Gram-Schmidth modificado, mínimos cuadrados y ARX presentan una diferencia pequeña en el tiempo de ejecución.

132

Capítulo 7: Conclusiones y trabajos futuros

Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARX

0.2

ARX LS

Tiempo de Ejecución (s)

QR GMC 0.15

0.1

GM

M

GIVENS HOUSEHOLDER AUDI L1

0.05

0 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.6 Acercamiento de los tiempos de ejecución de los algoritmos ARX, LS, QR, AUDI y L1 para una estructura ARX.

Tiempo de Ejecución (s)

Los tiempos de ejecución de la familia de algoritmos QR, y los algoritmos ARMAX, LS, , AUDI y L1 para una estructura tipo ARMAX se muestran en la siguiente figura. Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARMAX 2.5 ARMAX LS QR GMC GMM 2 GIVENS HOUSEHOLDER AUDI L1 1.5

1

0.5

0 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.7 Tiempos de ejecución de los algoritmos ARMAX, LS, QR, AUDI y L1 para una estructura ARMAX.

Los algoritmos que utilizan menos tiempo de ejecución son QR mediante Gram-Schmidth clásico y modificado, LS y QR de MATLAB®. Esto se puede ver más claro en la Figura 7.8.

133

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas. Tiempos de ejecución de los algoritmos para cada caso de estudio para la estructura ARMAX 0.25

Tiempo de Ejecución (s)

0.2

0.15

0.1 ARMAX LS QR GMC

0.05

GMM GIVENS HOUSEHOLDER AUDI L1

0

-0.05 Tercer orden

RCN100

Sistema de calentamiento Procesos

Secadora de pelo

Figura 7.8 Acercamiento de los tiempos de ejecución de los algoritmos ARMAX, LS, QR, AUDI y L1 para una estructura ARMAX

A pesar de que las rutinas con las que cuenta el toolbox de Matlab son pre-compiladas en lenguaje C, es decir, son rutinas que permiten que el tiempo de ejecución sea menor que un cualquier archivo m-file, en la mayoría de los procesos, los algoritmos alternativos tienen un tiempo de ejecución menor. Si se toma en cuenta que el proceso de identificación es iterativo, y dependiendo del tipo de estructura y el orden que se desea estimar, el número de ciclos que se necesita para obtener el orden óptimo son muchos, estos algoritmos alternativos pueden agilizar ese proceso que puede resultar muy tedioso cuando el orden del modelo estimado es muy grande. La elección del mejor modelo para cada proceso es un tanto complicada ya que se pueden tomar en cuentas muchos aspectos y en este trabajo solo se estimaron los modelos con la finalidad de validar los algoritmos propuestos, ya que ningún modelo estimado tiene un fin aplicativo. En la práctica, por ejemplo, si el modelo estimado se utiliza para el ajuste de algún controlador y éste satisface los requerimientos de control se puede decir que el modelo es válido para esa específica aplicación. Para tener una idea de una aplicación real, los modelos estimados con LU, QR, AUDI se pueden utilizar para tareas de predicción o simulación. Mientras que los modelos estimados con la norma L1 puede utilizarse en diseño y/o sintonización de controladores robustos. Por otro lado, la elección de la mejor estructura del modelo es un compromiso entre flexibilidad y ahorro [53]. Entiéndase por flexibilidad que la estructura del modelo tiene la capacidad para describir diferentes sistemas y una forma de obtenerlo es utilizando muchos parámetros. Por otro lado el ahorro requiere la utilización de pocos parámetros, que implica por supuesto menor tiempo de ejecución en la estimación de los mismos. 134

Capítulo 7: Conclusiones y trabajos futuros

7.2 Trabajos futuros Dentro del contexto de sistemas lineales sería interesante implementar estos algoritmos para las estructuras lineales que no se abordaron: OE, BJ, ARARX y ARARMAX. Además se puede hacer el esfuerzo por integrar estás técnicas de estimación de parámetros a sistemas no lineales, con estructuras NARX, NARMAX. Medir el tiempo de ejecución utilizando de una manera distinta para tener un grado de error mucho menor. Compilar las rutinas programadas para reducir el tiempo de ejecución con la finalidad de que la comparación, entre los algoritmos alternativos y las rutinas de estimación paramétrica (ARX, ARMAX, OE, BJ, etc.), sea bajo las mismas condiciones.

135

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

136

Referencias [01] Alan Jennings and J.J. McKeown, Matrix Computation, Second Edition, Wiley, 1992. [02] Ben Noble y James W. Daniel. Algebra Lineal Aplicada, Tercera Edición, Prentice Hall, 1989. [03] Bikal P. Singh, Carlos N. Bouza, Luís C. Martínez y Sira M. Allende, “El Ajuste de la Ecuación de Regresión en Problemas Económicos: Mínimos Cuadráticos vs. Mínima Desviación Absoluta Media,” Documentos de Trabajo en Análisis Económico, Vol. 3, No. 16, 2004. [04] Carlos Mario Vélez, “Curso de Identificación de Sistemas Dinámicos”, Maestría en Matemáticas aplicadas, Universidad EAFIT, Grupo de investigación en sistema de control digital, Colombia, 2004. http//www.control-system.net. [05] Christakis Charalambous, Adib Nashashibi, A computational Method for The LU Decomposition of Rectangular Matrices and its Implications to The Realization of 2-D Digital Filters. [06] Christian H. Bischof, “A Parallel QR Factorization Algorithm using Local Pivoting,” Conference on High Performance Networking and Proceedings of the 1988 ACM/IEEE conference on Supercomputing.

Computing

[07] David C. Yu and Huili Wang, “A new Parallel LU Decomposition Method,” IEEE Transactions on Power System, Vol. 5, No. 1, February 1990. [08] Daniel W. Apley and Jianjun Shi, “The inverse QR descomposition in order recursive calculation of least squares coefficients”, Dept. of Mechanical Engineering and Applied Mechanics, The University of Michigan, Ann Arbor, MI. 1995. [09] David C. Yu and Huili Wang, “A new parallel LU descomposition mtehod”, IEEE Transactions on Power Systems, Vol. 5, No. 1, February 1990. [10] Daysi (Database for the Identification of Systems), Developed, maintained, and hosted by SISTA. http//homes.esat.kuleuven.be/~smc/daisy. [11] Erik Elmorth and Fred Gustavson, “High-Performance Library software for QR Factorization,” Springer-Verlag, Berlin Heidelberg, 2001. 137

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

[12] Escuela Universitaria politécnica de Manresa, “Apuntes de Identificación de sistemas”, Departamento de ingeniería en sistemas automáticos e informática industrial, Julio 2004. http//www.eupm.upc.es/~esaii/assign/ident. [13] Gerald and Wheatley, Applied Numerical Analysis, Addison Wesley, Sixth Edition 1999. [14] Gerald Recktenwald, Numerical Methods with Matlab®, Implementation and Application, Prentice Hall, 2000. [15] Grewal M. S. and Andrews A. P., Kalman Filtering-Theory and Practice, Prentice Hall, Englewood Cliffs, New Jersey, 1993. [16] Hideo Yamashita and Eihachiro Nakamae, “A Pivot Ordering Algorithm Aimed at Minimizing Computation Time,” IEE Transactions on Circuits and Systems, Vol. CAS-25, No. 8, August 1978. [17] Jean Claude Carmona and Víctor Alvarado, “Active Noise control of a Duch using robust control theory”, IEEE Transactions on Control Systems Technology 8 (2000), No. 6, pag. 930-938. [18] Jeffrey C. Lagarias, James A. Reeds, Margaret H. Wright, Paul E. Wright, “Convergence Properties of the Nelder--Mead Simplex Method in Low Dimensions,” SIAM Journal Optimization, Vol. 9, No.1, 1998. [19] Jhon H. Mathews and Kurtis D. Fink, Métodos Numéricos con MATLAB, Tercera Edición, Prentice Hall, Madrid, 2000. [20] Jiandong Wang, “Identification of Multirate Systems and Closed-Loop Systems,”

Master’s Thesis, Department of Electrical and Computer Engineering, Edmonton, Alberta Canada, 2003. [21] J.C. Carmona and V.M. Alvarado, “L1 Predictor Error Approach in System

Identification,” Proceedings of the American Control Conference, Anchorage, AK, May 8-10, 2002. [22] J. M. Peña, “LDU Decompositions with L and U well conditioned,” Electronic

Transactions on Numerical Analysis. Vol. 18, pp. 198-208, 2004. [23] Karl Johan Ǻström and Björn Wittenmark, Adaptive Control, Addison Wesley. [24] Karl J. Ǻström and Björn Wittenmark, Computer-Controlled Systems Thoery and

Design, Third Edition, Prentice Hall.

138

Referencias

[25] Katsuhiko Ogata, Ingeniería de control Moderna, Prentice Hall, Segunda Edición,

1998. [26] K. N. Balasubramanya Murthy and C. Siva Ram Murthy, “A new Parallel Algorithm

for solving sparse linear systems”, Department of computer Science and Engineering Indian Institute of Tehcnology, Madras, India, 1995. [27] K. N. Balasubramanya Murthy and C. Siva Ram Murthy, “Parallel algorithm for solving sparse linear systems,” ELECTRONICS LETTERS, Vol. 2, No. 19, September 1996. [28] Laurence V. Fausett, Appiled Numerical Analysis using MATLAB®, Prentice Hall, 1999. [29] Lennart Ljung, “Issues in System identification”, IEEE Control systems, 02721708/91/0100-0029, 1991. [30] Lennart Ljung, “System Identification,” Department of Electrical Engineering, Linköping University S-581 83 Linköping, Sweden, May 29, 1995. [31] Lennart Ljung, System Identification Theory for the user, Prentice Hall: Englewood Cliffs, NJ. 1999. [32] Martin Golubitsky y Michael Dellnitz, Álgebra lineal y ecuaciones diferenciales, con uso de Matlab®, Thomson Learning, 2001. [33] MATLAB® 7.0 (R14), The Language of Technical Computing, The Mathworks Inc, 2004. [34] MATLAB®, User´s Guide. [35] Mustafa Khammash, Murti V. Salpaka, and Tim Van Voorhis, “Robust Synthesis in L1: A Globally Optima Solution,” IEEE Transactions on Automatic Control, Vol. 46, No. 11, November 2001. [36] Nicholas J. Higham, “Accuracy and Stability of Numerical Algorithms”, SIAM, 1996. [37] Oliver Y. K. Chen, “A Comparison of Pivoting strategies for the Direct LU Factorization,” Department of Mathematics, Montana State University – Billings, billings, MT 59101-0298. [38] Oriol Jorba I Casellas, “Simulación de los campos de viento de la península Ibérica y el área geográfica de Catalunya con alta resolución especial para distintas situaciones meteorológicas típicas”, Tesis de Doctorado, Departamento de doctorado en ingeniería Ambiental, Universitat Politècnica de Catalunya, Barcelona, Marzo de 2005. 139

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

[39] Petko Yanev, Paolo Foschi and Erricos John Kontoghiorghes, “Algoritms for Computing the QR Decomposition of a Set of Matrices with Common Columns,” Algorithmica, Springer-Verlag, New york Inc., 2004. [40] Qifa Ke and Takeo Kanade, “Robust Subspace Computation Using L1 Norm”, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA 15213, CMU-CS03-172, August 2003. [41] Ramón Doallo, Juan Touriño y Amilio L. Zapata, “Sparse Householder QR Factorization on a Mesh,” Parallel and Distributed Processing, Proceedings of the Fourth Euromicro Workshop on, 1996. [42] Ricardo Verdeja Higareda, “Comparación experimental del desempeño de un filtro activo paralelo de potencia bajo la acción de controladores tipo: Histéresis, proporcional-integral-derivativo y basado en pasividad”, Tesis de Maestría, Departamento de Electrónica, Cenidet, Cuernavaca, Diciembre de 2004. [43] Seymour Lipschutz, Algebra Lineal, Serie Schaum, McGraw-hill, Segunda Edición, 1992. [44] Shaohua Niu and D. Grant Fisher, “A Multiple Model Least-Squares Estimation Method”, Department of Chemical Engineering, University of Alberta, Edmonton, Canada, T6G2G6. 1994. [45] Sigurd Skogestad and Ian Postlethwaite, Multivariable Feedback Control Analysis and Design, Wiley, 1996. [46] SISTA, Database for Identification of Systems, 2004. [47] Steve Niu, “Monitoring Parameter Identifiability With AUDI,” June 24, 1994 [48] Steve S. Niu, D. Grant Fisher, Lennart Ljung and L. shah, “A Tutorial On Multiple Model Least-Squares and Augmented UD Identification,” December 21, 1994. [49] Steve S. Niu, Lennart Ljung and Åke Björck,“Descomposition Methods for Solving Least-Squares Parameter Estimation”, IEEE Transactions on Signal Processing, Vol.44, No. 11, November 1996. [50] S. Niu, “Augmented ud identification for process control”, Ph.D. Thesis, University of Alberta, Edmonton, Canada, 1994. [51] Söderström, T. and Stoica, P., System Identification, Prentice Hall, Englewood Clifs, New Jersey, 1989. [52] The Mathworks, ‘System identification toolbox”, 2004. 140

Referencias

[53] Todd K. Moon and Wynn C. Stirling, Mathematical Methods and Algorithms, Prentice Hall, 2000. [54] Toolbox Identification, MATLAB® Inc. 2004. [55] Toolbox Optimización, MATLAB® Inc., 2005. [56] Universidad Nacional de Quilmas, “Apuntes de Identificación y Control Adaptivo”, Área de Automatización y Control Industrial, Buenos Aires, Argentina, Junio 2004. [57] Yu-fai Fung, Wai-leung Cheung, Michael G. Singh and Muhammet F. Ercan, “A Pc Based Parallel LU Decomposition Algorithm for Sparse Matrices,” Department of Electrical Engineering, The Hong Kong Polytechnic Unversity, Hung Hom, Hong Kong. 2003.

141

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

142

Apéndice A Descripción de Procesos A.1 Sistema de tercer orden. Este proceso es un sistema propuesto, tomado de [25], con la finalidad de cotejar los parámetros estimados con los reales. Este proceso tiene una respuesta semejante a la de un sistema de segundo orden por lo que la obtención de los parámetros de la respuesta transitoria se muestran en el Apéndice 0.

Figura A. 1 Sistema de tercer orden

A.2 Banco de nivel RCN 100. Este banco de pruebas se encuentra en el laboratorio de control de cenidet. Modo de operación: •

El proceso consiste en regular el nivel en el recipiente de acrílico (1). La válvula de control (2), actúa sobre la entrada del tanque, permitiendo la entrada de agua.

•

Cuenta con un sensor de nivel el cual es de tipo ultrasónico (3). Y varía de forma lineal desde 1.5 volts (tanque vacío) hasta 0.6 volts (tanque lleno hasta la altura de 18cms).

143

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

Figura A. 2 Banco de nivel RCN100.

Características Técnicas: •

Llave de bola de tres vías paso L en PVC DN15.

•

Caudalímetro (4) de sección variable 10 litros/min, presición ± 3%, de la escala completa.

•

Sensor de nivel de ultrasonido. Escala de 0.2-1 metro con una señal de salida de 05v. Frecuencia 200khz. Angulo de emisión 5º.

•

Micro válvula de regulación. Cuerpo, asiento y cierre de acero inoxidable 316. Cv 1.25. Escala de 40:1 normalmente cerrada.

•

Convertidor de corriente a presión. Señal de entrada 4-20mA. Señal de salida de 0.2 a 1 bario. Precisión de ± 0.54% de la escala completa.

A.3 Secadora de pelo (Dryer) Experimento de laboratorio que se comporta como una secadora de pelo. El aire se ventila a través de un tubo y se calienta a la entrada. La temperatura del aire se mide con un termopar a la salida del tubo. La entrada es el voltaje que se aplica al dispositivo de calentamiento (una malla resistiva).

144

Apéndice A

Figura A. 3 Secadora de pelo (dryer).

A.4 Sistema de calentamiento (Heating System) El experimento es un sistema de calentamiento tipo SISO. La entrada acciona una lámpara de Halógeno, la cual se encuentra suspendida varios centímetros sobre un plato delgado de acero. La salida es la medición del termopar ubicado en la parte posterior de plato.

Figura A. 4 Sistema de calentamiento (heating system).

145

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

146

Apéndice B Selección del periodo de muestreo B.1 Obtención del periodo de muestreo para el sistema de tercer orden. La Figura B. 1 muestra la respuesta al escalón del primer caso de estudio representado por (6.18), por lo tanto puede representarse mediante el modelo paramétrico en el dominio del tiempo siguiente G( s) =

s

2

ωn 2

+

A 2ξ s

ωn

e− Ls

(B.1)

+1

donde ξ es el coeficiente de amortiguamiento y ωn la frecuencia natural del sistema.

Figura B. 1 Respuesta al escalón.

Con base a la respuesta al escalón del sistema se obtiene el periodo adecuado y necesario para transformar el sistema a su representación correspondiente en tiempo discreto. La obtención de dicho periodo de muestreo depende de algunos parámetros con base a la respuesta transitoria, como son: El tiempo de establecimiento ts , tiempo de subida tr ,

147

Evaluación comparativa de otras alternativas al algoritmo clásico de mínimos cuadrados para la identificación de sistemas.

amortiguamiento ζ , sobrepaso máximo M p [25]. Generalmente se toman solo los parámetros de sobrepaso máximo y tiempo de establecimiento. Para obtener un periodo de muestreo adecuado se utiliza la siguiente relación descrita en [24]: Nr =

tr h

(B.2)

donde N r es el número de periodos de muestreo por cada tiempo de subida tr y tr es el tiempo de subida, en otras palabras, es el tiempo requerido para que la respuesta pase del 10 al 90%, del 5 al 95% o del 10 al 100% de su valor final. Un buen criterio es muestrear de 4 a 10 veces más rápido que la constante de tiempo. Para sistemas de primer orden el tiempo de establecimiento es igual a la constante de tiempo que tenga el sistema. En cambio para sistema de segundo orden el tiempo de establecimiento queda sujeto al coeficiente de amortiguamiento ζ y a la frecuencia natural del sistema ω0 . Por lo tanto la expresión que determina el tiempo de establecimiento ts para sistemas de segundo orden queda expresada como Tr = ω0 −1eϕ / tan ϕ

(B.3)

donde: ζ = cos ϕ y ϕ es la frecuencia amortiguada. Utilizando la función damp se obtiene la frecuencia natural y el coeficiente de amortiguamiento del sistema en cuestión, en este caso el sistema (6.18), y utilizando las ecuaciones (B.3) y (B.2) se obtiene el periodo de muestreo correspondiente.

148