Repositorio Institucional de Documentos

Proyecto Fin de Carrera Algoritmo de simulación de Dymola en Matlab para investigación de Modelos de Usuario deterministas en su aplicación a edificios Autor: Luis Jarque Catalán Director: Marcus Fuchs Codirectora: Sara Manzano Martínez Ponente: Carlos Monné Bailo Ingeniería Industrial Escuela de Ingeniería y Arquitectura Mayo 2013 1 / 2 Resumen ALGORITMO DE SIMULACIÓN DE DYMOLA EN MATLAB PARA INVESTIGACIÓN DE MODELOS DE USUARIO DETERMINISTAS EN SU APLICACIÓN A EDIFICIOS Los modelos de representación de edificios son una herramienta importante para simular su comportamiento dinámico. La creación de dicho modelo requiere una gran labor de esfuerzo y tiempo, pero el objetivo principal es el estudio de la influencia de los parámetros sobre el edificio para su control o mejora. Hasta el momento se han utilizado los perfiles de usuario aplicando lo establecido en la norma DIN 18599 y SIA 2014 para estimar las cargas internas del edificio debidas a personas, máquinas y luz y con ellas se simula el comportamiento energético del edificio. Este proyecto aborda la estimación de los parámetros debido a cargas internas a partir de medidas reales de energía consumida de calefacción. La optimización consiste en un proceso iterativo de calibración del resultado de energía de calefacción de la simulación del modelo con datos reales de consumo. Este proceso se realiza mediante tres métodos. El primer método corresponde a una función implementada dentro del programa de simulación Dymola. El segundo método es un programa desarrollado con Matlab que permite acoplar Dymola para la optimización de parámetros. El tercero consiste en el desarrollo de un algoritmo de Matlab que implementa el método matemático de Gauss-Newton por pasos. Tras el desarrollo y aplicación de los métodos correspondientes se estudia su comportamien- to sobre el modelo del edificio de referencia. El estudio sobre uno o varios parámetros del modelo de usuario permite obtener los resultados del proyecto y compararlos. Mediante el análisis de sensibilidad se estudia la dependencia entre parámetros del modelo de usuario y gracias a los resultados se concluye que los dos programas desarrollados con Matlab convergen con la solución óptima en menos tiempo y con menos iteraciones que Dymola. El proyecto se ha desarrollado en el Centro de Investigación E.ON para el Instituto de Eficiencia Energética en Edificios y Clima Interior de Aachen (Alemania) en el ámbito del desarrollo de una herramienta de planificación integral de energía para un campus de Alemania. El presente proyecto puede ayudar en la estimación de parámetros de los edificios estudiados. II Agradecimientos Agradezco al centro de investigación E.ON por haberme dado la oportunidad de desarrollar mi experiencia profesional en el ámbito de las energías. Mencionar a Marcus Fuchs por su seguimiento constante a lo largo del proyecto. Agradezco a mi codirectora Sara Manzano por su disponibilidad e interés en todo momento. En especial a mi familia, Lola y Jacinto, vuestro esfuerzo diario me ha ayudado en la lucha por mis objetivos. A ti Nieves, por todo lo que me has demostrado a lo largo de estos dos años. A mis amigos. Muchas gracias. III Índice general Resumen II Agradecimientos III Table of Contents IV List of Figures VII List of Tables IX 1. Introducción 1 1.1. Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3. Estructura de la tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2. Métodos de optimización 5 2.1. Visión general de las optimizaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2. Edificio de referencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2.1. Métodos elegidos para la estimación de parámetros . . . . . . . . . . . . . 8 2.3. Calibración del modelo con Dymola . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1. Configuración de la calibración . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.2. Ajuste de parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4. Herramienta de optimización de Matlab . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.1. Conexión entre Dymola y Matlab . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4.2. Algoritmo y optimización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.3. Precisión del algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5. Método de optimización matemática . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3. Desarrollo de optimizaciones 19 3.1. Modelo de edificio simplificado: J1615 como una zona térmica . . . . . . . . . . . 19 3.2. Estudio de parámetros fijos para el modelo de una zona . . . . . . . . . . . . . . . 20 IV Índice general 3.2.1. Un parámetro: flujo de calor fijo . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2.2. Tres parámetros: personas, máquinas y ratio de luz . . . . . . . . . . . . . . 22 3.2.3. Once parámetros: perfil de usuario de personas . . . . . . . . . . . . . . . . 23 3.3. Modelo de edificio: J1615 como seis zonas térmicas . . . . . . . . . . . . . . . . . . 24 3.4. Estudio de parámetros fijos para el modelo de seis zonas . . . . . . . . . . . . . . 26 3.4.1. Tres parámetros: personas, máquinas y ratio de luz . . . . . . . . . . . . . . 26 3.4.2. Seis zonas con flujo de potencia fijo . . . . . . . . . . . . . . . . . . . . . . . 27 4. Resultados 29 4.1. Análisis de sensibilidad y peso de los parámetros . . . . . . . . . . . . . . . . . . . 29 4.2. Comparación de métodos y simulaciones . . . . . . . . . . . . . . . . . . . . . . . . 31 5. Conclusión 34 Anexo A: Código del script de optimización 36 Anexo B: Código de la función del algoritmo matemático 40 Bibliografía 44 English version 45 A. Introduction I A.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I A.2. Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I A.3. Structure of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV B. Optimization methods V B.1. General optimization view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V B.2. Reference building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VI B.2.1. Methods chosen for parameter estimation . . . . . . . . . . . . . . . . . . . VIII B.3. Model Calibration with Dymola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X B.3.1. Calibration setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X B.3.2. Tune the parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X B.4. Matlab Optimization toolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XV B.4.1. Dymola to Matlab connection . . . . . . . . . . . . . . . . . . . . . . . . . . XVI B.4.2. Algorithm and solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVII B.4.3. Algorithm precision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIX V Índice general B.5. Mathematical optimization method . . . . . . . . . . . . . . . . . . . . . . . . . . . XX B.5.1. Gauss-Newton method (Levenberg-Marquardt) . . . . . . . . . . . . . . . . XXI C. Development of optimizations XXIV C.1. Simplified building model: J1615 as one thermal zone . . . . . . . . . . . . . . . . XXIV C.2. Study of fixed parameters in one thermal zone model . . . . . . . . . . . . . . . . XXV C.2.1. One parameter: fixed heat flow . . . . . . . . . . . . . . . . . . . . . . . . . . XXV C.2.2. Three parameters: people, machines and light ratio . . . . . . . . . . . . . . XXVI C.2.3. Eleven parameters: user profile people . . . . . . . . . . . . . . . . . . . . . XXVII C.3. Building model: J1615 as six thermal zones . . . . . . . . . . . . . . . . . . . . . . . XXVIII C.4. Study of fixed parameters in six thermal zone model . . . . . . . . . . . . . . . . . XXX C.4.1. Three parameters: people, machines and Light ratio . . . . . . . . . . . . . XXX C.4.2. Six zones with fixed heat flow . . . . . . . . . . . . . . . . . . . . . . . . . . . XXX D. Results XXXII D.1. Sensibility analysis and weight of parameters . . . . . . . . . . . . . . . . . . . . . XXXII D.2. Comparison of methods and simulations . . . . . . . . . . . . . . . . . . . . . . . . XXXIV E. Conclusion XXXVII Annex A: Code for optimization algorithm script XXXIX Annex B: Code for mathematical algorithm function XLIII Bibliografía XLVII VI Índice de figuras 1.1. Modelo de usuario del edificio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2. Contribucion de energia de las cargas internas durante una semana . . . . . . . . 4 2.1. Datos medidos filtrados por horas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2. Cargas internas correspondientes al modelo complejo de seis zonas . . . . . . . . 7 2.3. Diferentes zonas en las que se clasifica el edificio . . . . . . . . . . . . . . . . . . . 8 2.4. Interaction between parameters and demand . . . . . . . . . . . . . . . . . . . . . 9 2.5. Resultados de validación para una semana . . . . . . . . . . . . . . . . . . . . . . . 11 2.6. Resultados de calibración para una semana . . . . . . . . . . . . . . . . . . . . . . 12 2.7. Basis of pattern search method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.8. Conexión entre Dymola - Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.9. Evaluación de las opciones para mejora del algoritmo . . . . . . . . . . . . . . . . 14 2.10.Programa desarrollado en Matlab para estimación de parámetros . . . . . . . . . 15 2.11.Valor de la función y evaluación en función de las iteraciones con Pattern Search 16 2.12.Método matemático desarrollado como una función de Matlab . . . . . . . . . . 18 3.1. Edificio de referencia J1615 reducido a una zona . . . . . . . . . . . . . . . . . . . 20 3.2. Simplificación de parámetros en el edificio de referencia . . . . . . . . . . . . . . 21 3.3. Resultados dinámicos de potencia térmica, valores de los parámetros y potencia de las cargas internas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4. User profile for one zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.5. Reparto de energía total aportada por las cargas internas . . . . . . . . . . . . . . 25 3.6. Edificio de referencia J1615 con seis zonas . . . . . . . . . . . . . . . . . . . . . . . 26 3.7. Fallo de calibracion con la función de Dymola . . . . . . . . . . . . . . . . . . . . . 28 4.1. Solution after sweep two parameter on one zone building . . . . . . . . . . . . . . 30 4.2. Variable parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3. Simulación de la energía acumulada para 3 casos durante un año . . . . . . . . . 33 A.1. User model in building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II VII Índice de figuras A.2. Energy contribution of inner loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV B.1. Measurement data collected in hours . . . . . . . . . . . . . . . . . . . . . . . . . . VII B.2. Inner loads corresponding to complex model of six zones . . . . . . . . . . . . . . VII B.3. Zones of building . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIII B.4. Interaction between parameters and demand . . . . . . . . . . . . . . . . . . . . . IX B.5. Model specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XI B.6. Case of studio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XII B.7. Details of measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XII B.8. Time interval for simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIII B.9. Results of Validation for a week . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIII B.10.Tuner parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIV B.11.Results of calibration for a week . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XV B.12.Basis of pattern search method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVI B.13.Dymola - Matlab conecction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVII B.14.Options evaluation for algorithm in a simple one week study . . . . . . . . . . . . XVII B.15.Algorithm running with Optimization toolbox . . . . . . . . . . . . . . . . . . . . . XVIII B.16.Function value and evaluation with pattern search . . . . . . . . . . . . . . . . . . XIX B.17.Mathematical method as function in Matlab . . . . . . . . . . . . . . . . . . . . . . XXIII C.1. J1615 one zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXIV C.2. J1615 people heat contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXV C.3. Dynamic results of heating power, parameter values and inner loads heat power XXVII C.4. User profile for one zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXVIII C.5. Total inner loads energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXIX C.6. J1615 model as six zones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXIX C.7. Failure in calibration with Dymola for six parameters . . . . . . . . . . . . . . . . . XXXI D.1. Solution after sweep two parameter on one zone building . . . . . . . . . . . . . . XXXIII D.2. Variable parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXXIV D.3. Heating energy simulation for first four cases during one year . . . . . . . . . . . . XXXVI VIII Índice de cuadros 2.1. Estimación de parámetros de cargas internas . . . . . . . . . . . . . . . . . . . . . 8 2.2. Presencia de personas en tanto por uno de forma horaria . . . . . . . . . . . . . . 17 3.1. Resultados de calibración de una zona y un parámetro . . . . . . . . . . . . . . . . 22 3.2. Resultados de calibración de un año en el modelo de una zona y tres parámetros 22 3.3. Calibration for eleven parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4. Calibración tres parámetros en el modelo de seis zonas para un año . . . . . . . . 27 B.1. Building zones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIII B.2. Percentage of people hourly on one zone building . . . . . . . . . . . . . . . . . . XIX C.1. Results of one zone calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXVI C.2. Calibration of one zone with three parameters for one year . . . . . . . . . . . . . XXVII C.3. Calibration for eleven parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . XXVIII C.4. Calibration of six zone with three parameters for one year . . . . . . . . . . . . . . XXX IX Métodos de optimización 2.2 Edificio de referencia Figura 2.1: Datos medidos filtrados por horas (a) Personas (b) Máquinas (c) Luces Figura 2.2: Cargas internas correspondientes al modelo complejo de seis zonas 7 Métodos de optimización 2.2 Edificio de referencia (a) Planta baja (b) Primera planta Figura 2.3: Diferentes zonas en las que se clasifica el edificio La tabla (tab: 2.1) muestra la estimación de parámetros de personas, flujo de calor proporcionado por máquinas y el ratio de luz. Los valores correspondientes a al ratio se basan en la norma DIN 18599 y SIA 2014 y son los usados en el modelo creado como parámetros fijos. Sabiendo el consumo real del edificio, esos parámetros tienen que calibrarse en el sentido de reducir el cuadrado de la suma de doferencias entre consumo de energía para calefacción y consumo simulado por el modelo. Cuadro 2.1: Estimación de parámetros de cargas internas Nr People (m2/p)Machines (W/m2)Lights(W/m2)Area (m2)Persons Machines Offices 14.0 7.0 15.9 371 26.5 26 Corridor 0.0 0.0 3.0 668.7 0.1 1 Meeting1 3.0 7.0 15.9 41 13.65 3 Meeting2 3.0 6.5 15.9 107.5 35.85 7 Lecture 3.2 4.7 13.0 105.8 33 5 Server 30.0 660.0 7.1 7.58 0.25 50 2.2.1. Métodos elegidos para la estimación de parámetros El modelo de usuario estimado en 2.1 se usa como valores iniciales de nuestros parámetros antes de la estimación. Ésto lleva a incertidumbres en la simulación, sin embargo la comparación es una referencia de la habilidad de simulación para reflejar un edificio real con todas sus influencias. 8 Métodos de optimización 2.2 Edificio de referencia Figura 2.4: Interaction between parameters and demand Los modelos estudiados contienen un gran número de parámetros cuyo comportamiento no se puede predecir. Los métodos de solución numérica en Modelica permiten reutilizar los modelos de simulación para optimización. Además los problemas de optimización dinámica no lineal se pueden tratar de forma eficiente como problemas discretos en el tiempo y resuelto numéricamente por aplicación de métodos no lineales de optimización a gran escala.[Bertsekas, 1999] El objetivo es analizar diferentes métodos y comparar su ajuste para minimizar la diferencia anual demanda de calefacción en el edificio de referencia. El modelo estudiado tiene todos los parámetros fijados y la función objetivo consiste en reducir un escalar calculado como la suma de diferencias al cuadrado. Además, los métodos tienen que considerar el consumo de energía horario sin ecuaciones diferenciales. Así pues, los métodos elegidos para estudiar son los siguientes por las ventajas que se muestran: 1. Función de calibración en Dymola .La función de calibración a está implementada en una librería de Dymola. .Rapidez y facilidad de preparación para la calibración. .Modelica programa el modelo y Dymola lo simula. 2. Herramienta de optimización de Matlab. 9 Métodos de optimización 2.3 Calibración del modelo con Dymola .El método Pattern Search representa una subclase de algoritmo de búsqueda directa en el cual se busca el minimizado de una función continua sin el uso de derivadas. [Audet u. J. E. Dennis, 2000] .Numéricamente robusto abordando los criterios no suaves. [H. Elmqvist, 2005] .Función de optimización implementada en Matlab. .Ajuste de los limites de los parámetros. 3. Método matemático en Matlab (Levenberg-Marquardt method) .Minimiza la función objetivo de forma no diferenciable. .Aproximación del gradiente por diferencias finitas. .Método fiable y rígido para problemas de alto nivel. .Convierte el problema no lineal en uno lineal 2.3. Calibración del modelo con Dymola La calibración del modelo consiste en la estimación de parámetros. Es este proceso los datos medidos de un dispositivo real se usan paraajustar parámetros de tal forma que los resultados de simulación tengan un buen acuerdo con los datos de medida. Los parámetros que ajustamos están referidos en esta herramienta como sintonizadores. Dymola varía los sintonizadores y simula para la búsqueda de soluciones satisfactorias. Matemáticamente, el proceso de ajuste de parámetros es un proceso de optimización para minimizar el error entre resultados de simulación y medidas a través de iteraciones.[AB, 2012] 2.3.1. Configuración de la calibración Antes de sintonizar los parámetros de las medidas, debemos estudiar qué parámetros pueden ser estimados de las medidas disponibles. Cambiando un parámetro a estimar, debe influir en la salida. Sin embargo, dos o más parámetros pueden influir en el resultados de una forma similar (covarianza), los cuales no son posibles estimar individualmente. [H. Elmqvist, 2005] 2.3.2. Ajuste de parámetros El proceso de ajuste de parámetros mediante la herramienta de Dymola se explica detalladamente en la versión de la memoria en inglés (B.3.1) ajustando los parámetros de optimiza- 10 Métodos de optimización 2.3 Calibración del modelo con Dymola ción para un modelo creado a partir del modelo de referencia, en el cual se han reducido las seis zonas a una con las mismas características. Para ajustar los parámetros, primero se realiza una validación de la simulación. Ésto consiste en determinar si el modelo es una representación precisa de un sistema real y el grado en el que se asemeja sin tener que ajustar parámetros. Tras la validación (fig: 2.5), se calcula un error del criterio de 4,7e6 y una diferencia de energía total acumulada de 21.5%. Figura 2.5: Resultados de validación para una semana Una vez validado el modelo hay que seleccionar los parámetros para la calibración, que en este caso son los tres parámetros de las cargas internas estimadas con las normas DIN 18599 and SIA 2024 (tab: 2.1) donde el valor de los tres parámetros es (110, 41, 15.9) y fijamos los límites superior e inferior, se puede calibrar el modelo. Tras 13 iteraciones Dymola devuelve una pantalla con los resultados. En el que se muestra el error del criterio (2,44e5). La diferencia es este caso de energía total acumulada se ha reducido hasta un 7% (fig: 2.6). 11 Métodos de optimización 2.4 Herramienta de optimización de Matlab La sincronización afecta a los resultados disminuyendo el error del criterio en un 94,8%. Figura 2.6: Resultados de calibración para una semana 2.4. Herramienta de optimización de Matlab La herramienta de optimización proporciona un amplio uso de algoritmos para optimización estándar y gran escala. Esos algoritmos resuelven problemas discretos continuos con y sin restricciones. La herramienta incluye funciones para programación lineal, cuadrática, optimización no lineal, mínimos cuadrados no lineales, sistemas de ecuaciones no lineales y optimización multiobjetivo.[The MathWorks, 2003] El método Pattern Search se usa con la herramienta de optimización (fig: 2.7). Es un algoritmo basado en la búsqueda de el valor mínimo para la función objetivo. El algoritmo construye una malla que se evalúa de forma que actualiza el parámetro óptimo cada vez que una búsqueda se realiza con éxito. Poll ocurre cuando el paso buscado no está disponible para obtener un punto en la malla actual que disminuya el valor. Si no disminuye la función objetivo en los puntos de la malla alrededor de la iteración actual, la distancia entre puntos de la 12 Métodos de optimización 2.4 Herramienta de optimización de Matlab malla se reducen a la mitad ∆x=∆x/2, es decir, se refina la malla y el proceso se repite hasta encontrar el punto adecuado. [Audet u. J. E. Dennis, 2000] Figura 2.7: Basis of pattern search method 2.4.1. Conexión entre Dymola y Matlab La conexión entre ambos programas es una parte importante del proyecto. El script de Matlab ,dymolaM, hace que desde Matlab se pueda ejecutar cualquier comando que pueda usarse en Dymola. De esta forma podemos mandar las siguientes ordenes para calibrar el modelo desde Matlab: 1. Transladar el modelo de Modelica a Dymola 2. Ajustar los parámetros estimados por Pattern Search 3. Simulación del modelo 13 Métodos de optimización 2.4 Herramienta de optimización de Matlab Figura 2.8: Conexión entre Dymola - Matlab 2.4.2. Algoritmo y optimización El algoritmo de optimización desarrolado se ha desarrollado para ser más eficiente respec- to a reducción en el tiempo de cálculo, precisión y estabilidad de resultados. Sin embargo, esta herramienta necesita un gran esfuerzo para ser implantada (Ver código de Matlab en 5. Añadir también que la única forma de conseguir un programa eficiente y eficaz es ajustar correctamente las opciones de la herramienta, este trabajo conlleva un trabajo duro de análisis del comportamiento del algoritmo en Matlab. (fig: 2.9). Figura 2.9: Evaluación de las opciones para mejora del algoritmo Por otra parte, una vez implementada la herramienta se puede adaptar de forma sencilla para diferentes casos mediante los siguientes pasos: 1. Modificación del script 14 Métodos de optimización 2.4 Herramienta de optimización de Matlab .Dar los datos de medida real como fichero de matlab .mat (durante el proyecto no fue necesario cambiarlo). .Ajustar las condiciones iniciales del vector de parámetros (x0). .Ajustar los límites superior e interior (lb,ub) así como restricciones si fuera necesario (A,b,Aeq,beq). .Ajustar la tolerancia de los resultados de la función objetivo (Tolmesh : 1e−06). 2. Modificación de la función .Dar el nombre de los parámetros cuyos parametros se cambiarán en Dymola. .Cambiar las características de simulación (tiempo a simular y tiempo entre intervalos, tolerancia y solver). .Cambiar el criterio de la función objetivo a minimizar (|Total Heat −ydata|2). Figura 2.10: Programa desarrollado en Matlab para estimación de parámetros En la función de optimización (fig: 2.10), los parámetros iniciales se cargan y simulan en Dymola y se comparan con los datos medidos para obtener la función objetivo. El programa 15 Métodos de optimización 2.4 Herramienta de optimización de Matlab varía los parámetros y los cambia en Dymola para volver a simular y compara de nuevo con los datos medidos. A partir de ahí, el algoritmo elige el resultado de la función objetivo óptimo y sigue un proceso iterativo hasta que la resta de funciones objetivo es menor que 1e-6. Figura 2.11: Valor de la función y evaluación en función de las iteraciones con Pattern Search 2.4.3. Precisión del algoritmo El grado de optimización de parámetros se verifica mediante un proceso inverso. Es decir, con los parámetros por defecto del modelo se simula el modelo y la energía total de calefacción obtenida va a ser el input de datos medidos. Así pues, se quiere estimar 11 parámetros, correspondientes a las 11 horas que está abierto el edificio y que representan el porcentaje de la presencia de personas que hay en el edificio con un intervalo de una hora entre las 7 a.m. y las 18 a.m.. Se han estudiado dos casos, en el primero no se han fijado límites, es decir, la presencia de personas puede oscilar libremente entre 0% y 100% de la ocupación. La precisión de resultados respecto al modelo por defecto es del 96,8% para una semana, en cambio de aprecia como el algoritmo a estimado las primeras horas de la mañana con mucha mayor presencia que las de mitad de día. Por ello se ha estudiado el segundo caso donde los límites se han fijado a ±20%. Los resultados obtenidos reflejan una muy buena precisión del algoritmo del 99,1%. 16 Desarrollo de optimizaciones 3.2 Estudio de parámetros fijos para el modelo de una zona más tiempo e iteraciones. Los resultados de Dymola son interesantes porque los tres parámetros son inferiores a los obtenidos con los métodos de Matlab, esto es ilógico por una razón, el criterio para los tres métodos es la reducción de la suma de errores al cuadrado. En éste caso el método matemático es también tan preciso como el método de optimización de Matlab. Sin embargo, se a pesar de seguir siendo el más rápido, se necesita más tiempo para la estimación de los parámetros. La siguiente gráfica es el resultado de aplicar los tres parámetros estimados con el programa de optimización durante una semana (fig: 3.3). La primera gráfica representa la potencia instantánea del equipo de calefacción; mientras que la segunda representa el número de personas y máquinas de forma dinámica; y la última gráfica muestra la potencia generada por las personas, máquinas y luces. Figura 3.3: Resultados dinámicos de potencia térmica, valores de los parámetros y potencia de las cargas internas 3.2.3. Once parámetros: perfil de usuario de personas El edificio de oficinas está abierto entre las 07 a.m. y las 6 p.m. y la ocupación a lo largo del día cambia constantemente, por ello el estudio horario de cargas internas debidas a personas 23 Desarrollo de optimizaciones 3.3 Modelo de edificio: J1615 como seis zonas térmicas puede presentar una importante disminución en la función objetivo. Estas cargas las pueden modificar los programas creados con Matlab pero no la función de Dymola por tratarse de datos externos al programa de Dymola (fig: 3.4). Figura 3.4: User profile for one zone Los resultados para este caso se representan en el cuadro 3.3, donde la herramienta de optimización de Matlab gasta cinco veces más tiempo que el método matemático a pesar de que el error se reduce hasta prácticamente el mismo valor. El valor de los parámetros estimado por el método 2 es muy similar al del modelo por defecto debido a su restricción en los límites superior e inferior, mientras que el método 3 por no tener los límites establecidos tiene una estimación de parámetros muy diferida. Cuadro 3.3: Calibration for eleven parameters Hour 7h 8h 9h 10h 11h 12h 13h 14h 15h 16h 17h Default% ocup 0.2 0.4 0.6 0.8 0.8 0.4 0.6 0.8 0.8 0.4 0.2 Toolbox Matlab 0.278 0.398 0.590 0.780 0.780 0.390 0.59 0.827 0.782 0.398 0.205 math method 0.421 0.426 0.410 0.450 0.435 0.418 0.427 0.401 0.423 0.460 0.428 Error2(kW h) Iterations Time (sec) Toolbox Matlab 2.06e10 65 13750 math method 2.07e10 6 2821 3.3. Modelo de edificio: J1615 como seis zonas térmicas El modelo del edificio de referencia con seis zonas es el más complejo que se ha estudiado. Cada zona se comporta de forma diferente y posee su horario para el modelo de usuario. De esta forma, la energía total de calefacción se calcula como la suma de necesidades de cada 24 Desarrollo de optimizaciones 3.3 Modelo de edificio: J1615 como seis zonas térmicas una de las zonas que dependerá a su vez y entre otros parámetros del perfil de usuario de esa zona (fig: 3.6). Para la estimación de parámetros en este modelo se han fijado los porcentajes de cargas internas, de forma que solo es necesario estimar el número total de personas, máquinas y ratio de luz sobre el edificio y como ya se conoce la contribución de energía de cada zona y cada parámetro, los métodos desarrollados reparten el correspondiente del total a cada zona (fig: 3.5). El modelo complejo implica un aumento en el número de parámetros de la estructura hasta 4581 parámetros. Esto significa cinco veces más ecuaciones que el modelo simple de una zona, es decir, aproximadamente el doble de parámetros por cada zona de más que tiene el modelo. (a) Contribución debido a personas (b) Contribución debido a máquinas (c) Contribución debido a luces (d) Contribución de todas las cargas internas Figura 3.5: Reparto de energía total aportada por las cargas internas 25 Desarrollo de optimizaciones 3.4 Estudio de parámetros fijos para el modelo de seis zonas Figura 3.6: Edificio de referencia J1615 con seis zonas 3.4. Estudio de parámetros fijos para el modelo de seis zonas 3.4.1. Tres parámetros: personas, máquinas y ratio de luz La estimación de parámetros en cada unos de los métodos aplicados es muy diferente. Ésto muestra la gran posibilidad de combinación ente parámetros que convergen en una solución mínima. EL método matemático no es tan rápido comparado con el programa que utiliza el algoritmo Pattern Search y mucho más lento que los casos anteriores por la complejidad del modelo actual. La función de Dymola tiene el error más pequeño de los tres pero requiere 4 veces más tiempo y 9 iteraciones más que los dos programas de Matlab. 26 Desarrollo de optimizaciones 3.4 Estudio de parámetros fijos para el modelo de seis zonas Cuadro 3.4: Calibración tres parámetros en el modelo de seis zonas para un año Dymola Toolbox Matlab Math method Nr People 26 222 117 Machines 102 92 268 Light Ratio (W/m2) 25.0 14.5 13.7 Error2(kW h) 1.77e9 1.37e10 1.15e10 Iterations 23 14 14 Time (s) 11402 3420 3520 3.4.2. Seis zonas con flujo de potencia fijo Se trata del mismo caso que el estudiado para el modelo simple de edificio 3.2.1 con la diferencia de las seis zonas y que hay un parámetro de potencia fija para cada una de las zonas del edificio. El principal objetivo del caso es la estabilidad de los métodos, especialmente la convergencia de la función de calibración con Dymola. La calibración con Dymola fue imposible debido al fallo en la estimación de los parámetros (fig: 3.7). Ésto indica que el número máximo de parámetros que Dymola es capaz de estimar para nuestro modelo es 3. Para este caso ambos métodos presentan estabilidad durante la simulación y son capaces de obtener resultados, pero no es necesario estimar los parámetros, ya que el objetivo es la comparación de los tres métodos. 27 Desarrollo de optimizaciones 3.4 Estudio de parámetros fijos para el modelo de seis zonas Figura 3.7: Fallo de calibracion con la función de Dymola 28 4 Resultados 4.1. Análisis de sensibilidad y peso de los parámetros La dependencia entre parámetros y energía de calefacción no es linear, solo los modelos matemáticos consiguen hacer a los modelos lineales. Cambiando el número de personas a estimar en un edificio influye en la salida, sin embargo el número de máquinas y ratio de luz también influye en el resultado de forma que no se pueden analizar individualmente. Cuando un conjunto de parámetros se sincroniza, el modelo se valida primero y los parámetros se ajustan de nuevo para que haya buena relación entre medidas y resultados de simulación. Para unos datos medidos dados es posible conseguir buen ajuste incrementando la complejidad del modelo y el número de parámetros ajustados. Sin embargo, ésto no garantiza que el resultado también sea bueno para otras condiciones de operación. El modelo de una zona tiene unos valores iniciales de los parámetros de 110 personas, 41 máquinas and 15.9 W /m2para el ratio de luz. Tras el proceso que desarrolla el método de optimización (2.10) y teniendo en cuenta los datos reales de calefacción, los valores adecuados son respectivamente 113.5 personas, 33 máquinas y 14.25 W /m2. Por otro lado, mediante una función de Dymola que mide la sensibilidad de los parámetros se ha obtenido la siguiente dependencia entre parámetros: NrPeopleMachines +0,8258·Nr People +8,5919·LightingPower Esta ecuación relaciona los tres parámetros entre sí. De forma que no existe solución única que minimice la función objetivo sino una línea a lo largo del valle con posibles combinaciones entre número de personas y máquinas y fijando el ratio de luz a 14.25 (fig: 4.1). El peso de los parámetros indica como afecta el aumento en una unidad de cada uno de ellos sobre en las cargas internas. El peso de los tres parámetros son los siguientes: 7.9% para personas, 9.6% para máquinas y 82.5% para el ratio de luz. 29 Resultados 4.1 Análisis de sensibilidad y peso de los parámetros Figura 4.1: Solution after sweep two parameter on one zone building En el caso estudiado para este análisis, la estimación de cargas internas para un edificio de oficinas es más efectiva que para un laboratorio, porque hay valores estándar conocidos. En cambio existen factores estocásticos cruciales que son imposibles de prever; los dos más importantes son la ventilación manual y el ajuste del termostato de temperatura (fig: 2.4). El flujo de potencia de calefacción dependiente de los parámetros que varían con el tiempo se representa en la figura (4.2), y se puede observar que la influencia del intercambio de aire y las infiltraciones en el edificio hacen disminuir la temperatura en el sensor de temperatura por debajo de la fijada a 20,5oCdurante las horas de trabajo. Como consecuencia de esta disminución el equipo de calefacción se enciende hasta que la temperatura supera de nuevo los 20,5oC. Mientras que en el edificio no se trabaja, la temperatura de control se disminuye a 16oC, a pesar de que nunca baja por debajo de 19,5oC. 30 Resultados 4.2 Comparación de métodos y simulaciones Figura 4.2: Variable parameters 4.2. Comparación de métodos y simulaciones Gracias a las simulaciones de diferentes casos a lo largo del capítulo (3), se va podido analizar el comportamiento de tres métodos empleados. El análisis de los resultados de simulación pueden ayudar a determinar el mejor método usado para cada caso de estudio en función de las necesidades del usuario de programa. 1. Caso 1 (1 parámetro en una zona) (3.1). .La estimación de parámetros con el método matemático es muy rápida y converge rápido con la solución óptima para estudios de tiempo de una semana; además con una precisión mejor que la obtenida por el método de Dymola. .la precisión no es tan buena cuando el tiempo de estudio es de un año, a pesar de conseguirlo en escasos 5 minutos frente a las más de dos horas que necesita Dymola. 31 Resultados 4.2 Comparación de métodos y simulaciones 2. Caso 2 (3 parámetros en una zona) (3.2). .Para este caso los algoritmos creados en Matlab convergen en una solución donde los parámetros son diferentes, en cambio es un resultado muy lógico teniendo en cuenta que el ratio de luz tiene un peso ocho veces superior al de personas o máquinas. .El resultado de los parámetros obtenido con Dymola no es comparable porque es menor en los tres parámetros que los óptimos con los otros dos métodos. .La reducción en el número de parámetros fijos en el modelo de una zona permite que el método matemático sea el más rápido, exactamente nueves veces más que Dymola. 3. Caso 3 (11 parámetros en una zona) (3.3). .El error se reduce en ambos métodos hasta valores similares pero el método matemático necesita nueve veces menos tiempo y 59 iteraciones menos. .La estimación de parámetros para el método 2 de la herramienta de optimización con Matlab es muy similar a los parámetros por defecto del programa por su tolerancia del ±20%. 4. Caso 4 (3 parámetros en seis zonas) (3.4). .El tiempo necesario para el método matemático aumenta ligeramente por encima del tiempo necesario por el método 2. .El cuadrado del error se reduce con el método matemático en un 16.06%. .Cada método estima los parámetros de forma distinta. 5. Caso 5 (6 parámetros en seis zonas) (3.4.2). .El máximo número de parámetros que Dymola puede resolver es 3. Por eso Dymola devuelve un error con este estudio. La simulación de un modelo simplificado necesita 3 veces menos tiempo y 5 veces menos número de parámetros, a pesar de ésto, la diferencia en los resultados no es muy notable (D.3). Comparando los resultados de calibración del modelo, la línea azul es la simulación de referencia y las otras tres líneas son casos estudiados, donde el ajuste de energía total de la línea roja (caso 1) es del 88.92%. El caso 2 representado por la línea verde tiene una adaptación a la de referencia del 90.43%. Mientras que el caso 4 es el que mejor se ajusta a la curva de referencia con un 91.11%. 32 Conclusión Anexo A: Código del script de optimización outputInterval =3600, tolerance =1, resultFile =" J1615_onezone3Parameters") ’) ; 106 end 107 end 108 end 109 end 110 end 111 112 % Loads data from r e s u l t f i l e .mat in a structure ’d ’ . 113 d=dymload( [onezone , ’/optimierung/J1615_onezone3Parameters .mat ’ ]) ; 114 115 % Extract from structure ’d ’ values of simulation . 116 TotalHeat= dymget(d, ’House. idealHeaterCooler_var_setpoint . heatMeter .q_kwh’ ) ; 117 TotalHeat= TotalHeat ( 1:( length(TotalHeat)−1)) ; 118 119 F=sum(( TotalHeat−ydata ) .^2) 39 Anexo B: Código de la función del algoritmo matemático 1%%%%%%%%FUNCTION FOR MATHEMATICAL ALGORITHM %%%%%%%% 2function [ x_star ] = levenberg_marquardt11 3tic ; 4 5% Path of Dymola m f i l e s . 6DymolaPath = ’C:/ Program Files /Dymola 2013 FD01/Mfiles ’ ; 7 8% Path of workspace . 9MatlabPath = ’D:/mfu−lja /workspaces ’ ; 10 addpath(genpath(DymolaPath) ) ; 11 addpath(MatlabPath) ; 12 addpath ([ MatlabPath , ’/onezone ’ ] ) ; 13 addpath(genpath ([ MatlabPath , ’ /onezone/optimierung ’ ]) ) ; 14 15 % Change directory in Dymola. 16 onezone = [MatlabPath , ’/onezone ’ ] ; 17 dymolaM([ ’cd(" ’ ,onezone , ’/optimierung") ’ ]) ; 18 19 load ([ MatlabPath , ’/dataJ .mat ’ ]) ; % load 1615dinamicTable . csv data . 20 ydata=( dataJ (49:8809 ,5) /1000) −16540; % measurement data for one year . 21 22 % solve nonlinear lea s t squares problem using Levenberg−Marquardt 23 % algorithm . 24 % i n i t i a l guess . 25 x = [0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5]; 26 27 % increment in x for f i n i t e d if fe re nce approximation . 28 delta_x = [0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01]; 29 30 %great lambda means short stepsize . 31 lambda = 1; 32 33 %counter 40 Conclusión Anexo B: Código de la función del algoritmo matemático 34 counter = 0; 35 progress = 1; 36 37 while progress > 1e−6 38 %counter 39 counter = counter+1; 40 41 % Simulation of Total Heat at vector of parameters ’ x ’ . 42 % Call to dymola here with value of x 43 userprofil; % SCRIPT for editting the column corresponding to UserProfilesOffice . txt . 44 45 % Model translation 46 res=dymolaM( ’ translateModel (" Optimization2 . J1615_onezone3Parameters ") ’ ) ; 47 48 % Model simulation 49 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=0, stopTime=31536000, numberOfIntervals=0, outputInterval=3600, tolerance =1, resultFile ="J1615_onezone3Parameters") ’ ) ; 50 51 d=dymload( [onezone , ’/optimierung/J1615_onezone3Parameters .mat ’ ]) ; 52 TotalHeat= dymget(d, ’House. idealHeaterCooler_var_setpoint . heatMeter .q_kwh’ ) ; 53 TH= TotalHeat (1 :( length(TotalHeat)−1)) ; %[8761x1 ] 54 55 % Simulation of Total Heat at x+dx for each parameter . 56 for k = 1: length(x) 57 58 % x plus dx for parameter k . 59 x_plus_dx = x(k) + delta_x (k) ; 60 61 % Simulation of Total Heat at x plus dx . 62 x(k)=x_plus_dx ; 63 userprofil; 64 x(k)= x_plus_dx −delta_x (k) ; 65 res=dymolaM( ’ translateModel (" Optimization2 . J1615_onezone3Parameters ") ’ ) ; 66 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=0, stopTime=31536000, numberOfIntervals=0, outputInterval=3600, tolerance =1, resultFile ="J1615_onezone3Parameters") ’ ) ; 67 68 d=dymload( [onezone , ’/optimierung/J1615_onezone3Parameters .mat ’ ]) ; 69 TotalHeat= dymget(d, ’House. idealHeaterCooler_var_setpoint . heatMeter .q_kwh’ ) ; 70 TH_dx( : , k)= TotalHeat (1 :( length(TotalHeat)−1)) ; % TH_dx [8761x11 ] 41 Conclusión Anexo B: Código de la función del algoritmo matemático 71 72 F=sum(abs(TH_dx( : , k)−ydata ) .^2) %criterion 73 end 74 75 % Right hand side 76 RHS = ydata −TH; %[8761x1 ] 77 78 % gradient 79 for k = 1: length(x) 80 81 % Partial derivative in k−direction . 82 P(: , k) = (TH_dx( : , k) −TH) ./ delta_x (k) ; % P [8761x11 ] each column i s divided by the column delta_x 83 end 84 85 % stepsize 86 delta_x = (P’*P + lambda*diag (diag(P’*P) ) ) \(P’*RHS) ; % [11x1 ] = { ([ 11 x8761 ][8760x11 ] ) + [11x11 ] } \ [11x1 ] 87 88 % update 89 x = x + delta_x ; 90 91 % progress 92 progress = ( delta_x ’*delta_x ) /(x ’*x) ; % [1 x11 ][1 1 x1 ] / [1 x11 ][1 1 x1 ] 93 94 if mod( counter ,10) == 0 95 disp ([ ’Counter : ’ ,num2str( counter ) ]) ; 96 lambda = 0.1*lambda; 97 end 98 end 99 100 % optimal solution 101 x_star = x 102 103 % number of iterations 104 disp ([ ’Number of iterations : ’ ,num2str( counter ) ]) ; 105 106 tiempo = toc; 107 fprintf( ’The process took %d seconds ’ , tiempo) ; 108 109 42 Conclusión Anexo B: Código de la función del algoritmo matemático 110 %%%%%%%%SCRIPT USERPROFIL %%%%%%%% 111 f i l e = [ MatlabPath , ’ /onezone/optimierung/Tables/J1615 ’ ] ; 112 filename = [ file , ’/ UserProfilesOffice . txt ’ ] ; 113 114 Time=[0;3540;3600;7140;7200; . . . ;597600;601140;601200;604740;]; 115 People= [ 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; x (1) ; x (1) ; x(2) ; x (2) ; x (3) ; x (3) ; x (4) ; x (4) ; x (5) ; x (5) ; x(6) ; x(6) ; x (7) ; x (7) ; x(8) ; x(8) ; x (9) ; x(9) ; x(10) ; x(10) ; x(11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; 0; 0 ; 0 ; x (1) ; x (1) ; x (2) ; x (2) ; x (3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x(10) ; x(10) ; x(11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ;0 ; 0 ; 0; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; x (1) ; x (1) ; x (2) ; x (2) ; x(3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x(10) ; x(10) ; x (11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ;0 ; 0 ; 0; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; x (1) ; x (1) ; x (2) ; x(2) ; x(3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x(10) ; x (10) ; x(11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; 0; 0 ; 0 ;0 ; 0 ; x(1) ; x (1) ; x(2) ; x(2) ; x(3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x (10) ; x(10) ; x(11) ; x(11) ; 0; 0 ;0 ;0 ;0 ; . . . ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; ] ; 116 Machines = [0. 1; 0.1 ;0 .1; 0. 1; . . . ; 0 . 1 ; 0 . 1 ; 0 . 1 ; 0 . 1 ; ] ; 117 Light1=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1;1;0.3; ... ;0;0;0;0;]; 118 Light2=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.3;0.3;0.3; ... ;0;0;0;0;]; 119 A=[Time, People , Machines , Light1 , Light2 ] ; 120 121 122 fileID = fopen ([ file , ’/ UserProfilesOffice . txt ’ ] , ’wt ’ ) ; 123 fprintf( fileID , ’#1\n ’ ) ; 124 fprintf( fileID , ’double UserProfilesOffice (336 , 5) \n ’ ) ; 125 126 127 for i =1:length(A) 128 fprintf( fileID , ’ %f \ t %f \ t %f \ t %f \ t %f \n ’ ,A( i , : ) ) ; 129 end 130 131 fclose ( fileID ) 43 Bibliografía [AB 2012] AB, Dassault S.: Dymola: Dynamic Modelling Laboratory User Manual. Volume 2, 2012 [Audet u. J. E. Dennis 2000] AUDET, Charles ; J. E. DENNIS, JR: Pattern Search algorithms for mixed variable programming. (2000), S. 573–594 [Bernard P. Zeigler 2000] BERNARD P. ZEIGLER, Tag Gon K. Herbert Praehofer P. Herbert Praehofer: Theory of Modeling and Simulation: Integrating Discrete Event and Continuous Complex Dynamic System. (2000) [Bertsekas 1999] BERTSEKAS, Dimitri P.: Nonlinear Programming. 1999 [H. Elmqvist 2005] H. ELMQVIST, S.E. Mattsson D. Brueck C. Schweiger D. Joos M. O. H. Ollson O. H. Ollson: Optimization for Design and Parameter Estimation. (2005), S. 255–266 [M. Lauster 2012] M. LAUSTER, M. Fuchs R. Streblow D. M. J. Teichmann T. J. Teichmann: Dynamic building model for city quartier simulation. (2012) [Stefan Finsterle 2010] STEFAN FINSTERLE, Michael B. K.: A truncated Levenberg-Marquardt algorithm for the calibration of highly parameterized nonlinear models. (2010), S. 731–738 [The MathWorks 2003] THE MATHWORKS, Inc.: Optimization Toolbox Users Guide of Matlab, 2003 [Wetter u. Polak 2003] WETTER, Michael ; POLAK, Elijah: A convergence optimization method using pattern search algorithms with adaptative precision simulation. (2003), S. 1393–1400 [Wetter u. Wright 2003] WETTER, Michael ; WRIGHT, Jonathan: Comparison of a generalized pattern search and a genetic algorithm optimization method. (2003), S. 1400–1408 44 English version Diploma Thesis Optimization Procedures for Parameter Estimations in dynamic Simulations of Buildings Aachen, May 2013 Luis Jarque Catalán matriculation number: 317373 Supervisors: Marcus Fuchs Prof. Dr.-Ing. Dirk Müller RWTH Aachen University E.ON Energy Research Center | ERC Institut for Energy Efficient Buildings and Indoor Climate | EBC Mathieustraße 6, D-52074 Aachen 2/2 A Introduction A.1. Motivation Energy efficiency technologies are the key to reducing energy demand. In special, building sector displays a great potential in the reduction of heating and cooling energy. Dynamic simulations are the properly tool for this investigations. However, it is necessary to simplify building for simulation in a suitable level to achieve requirements of convergence, stability, calculation time and effort of parametrization. Usually, a lot of time is spent in creating the input for a simulation model, but once this is done, the user usually does not determine the parameter values that lead to optimal system performance. This can be because there is no time left to do the tedious process of changing input values, running the simulation, interpreting the new results and guessing how to change the input for the next trial, or because the system being analysed is so complex that the user is just not capable of understanding the non linear interactions of the various parameters. However, using mathematical programming, it is possible to do automatic single- or multi-parameter approximation with search techniques that require only little effort. A.2. Objective When a technician create a model of a real building, he tries to represent it as specific as possible. After the model is made some parameter can be optimized in order to reduce the complexity of model. For this action, we need real measured data. Then, comparison of model simulation and measured data can be evaluated. As studied model in this project is already built, all the parameters referring to materials are well know. However, presence of people, number of machines like computers and lights in each room have a stochastic behaviour and it is no possible to determine exactly. This three I A.2 Objective parameters compose inner loads and they contribute decreasing heating consumption. Along the project, we cite the concept Üser modelrefers to the influence and conduct of inner loads on building (A.1). Figura A.1: User model in building That User Models are considered deterministic because heating consumption evolution is determined by inner loads, that means, inputs influence outputs without taking into account of relations of equations. That is because there are no direct equations that relate inner loads and heating energy consumption. In such a way main task consists of tune simulation results and measured data of heat energy. The goal is to develop an optimization program and a mathematical method with Matlab to simulate through Dymola building behaviour along one year. After that, finding with less labour time the independent variables that yield better performance of such systems. This two methods are compare with another calibration tool implemented in Dymola. This tool provides the means for estimate parameters in the model to minimize error between simulation results and measurements through iterations. Also, for parameter optimization of II B.2 Reference building Figura B.4: Interaction between parameters and demand Models studied contains a huge number of parameter whose behaviour can not be predicted. Numerical solution methods in Modelica allows the reuse of simulation models for optimization. Also nonlinear dynamic optimization problems can be treated efficiently as discretetime optimal problem and solved numerically by applying largescale nonlinear optimization methods. The goal is to analyse different methods and compare their performance in minimizing the annual difference of heating energy consumption on reference building. The studied model has all parameters fixed and the objective function consists on a scalar calculated as the sum of squared difference. Furthermore, the methods have to consider energy consumption hourly without differential equations. The method chosen to study are: 1. Dymola calibration function .Calibration function already implemented in a library of Dymola. .Fast method for preparation of calibration. .Modelica for model creation and Dymola for simulations. 2. Matlab optimization toolbox .Pattern search methods represents a subclass of direct search algorithms, in which the minimizer of a continuous function is sought without the use of derivatives. [Audet u. J. E. Dennis, 2000] IX B.3 Model Calibration with Dymola .Numerically robust in tackling with non-smoth criteria. [H. Elmqvist, 2005] .Optimization function implements in Matlab. .Adjust of parameter bounds. 3. Mathematical method in Matlab (Levenberg-Marquardt method) .Minimize objective functions non differentiable. .Gradient approximation by finite difference. .Stable method for high level problems. .Convert a non linear problem to lineal. B.3. Model Calibration with Dymola Model calibration consists in a parameter estimation. In this process measured data from a real device is used to tune parameters such that the simulation results are in good agreement with the measured data. The parameters that we tune are often referred to as tuners. Dymola varies the tuners and simulates when it searches for satisfactory solutions. Mathematically, the tuning procedure is an optimization procedure to minimize the error between the simulation results and the measurements via iterations. B.3.1. Calibration setting Before tuning parameters from measurements, we have to study which parameters can be estimated from the available measurements. Changing a parameter to be estimated must of course influence the output. However, two or several parameters may influence the result in a similar way (co-variance), such that it is not possible to estimate them individually. When a set of parameters have been tuned, the model is validated and the tuned parameters against other measured data to check that there is a good agreement between the simulation result and the new measurements. [AB, 2012] B.3.2. Tune the parameters Previous to calibration of parameters, we have to validate the model. Validation is utilized to determine that a model is an accurate representation of the real system. Validation is usually X B.3 Model Calibration with Dymola achieved through the calibration of the model, an iterative process of comparing the model to actual system behavior and using the discrepancies between the two, and the insights gained, to improve the model. This process is repeated until model accuracy is judged to be acceptable. By way of example for tuning parameters, a simplified model with three parameter (C.1) has been chosen for validation and calibration of model. Let us first check how the model with nominal parameters compares with measured data.Validation is set up very similar to calibration. A basic difference is that no tunable parameters need to be specified for the validation. First, we specify the model to calibrate. Then, the model is translated in order to gather information needed to build browsers and selectors to support the remaining setting up of the calibration task (fig: B.5). Figura B.5: Model specification The next task is to specify the measurements and how they are stored (fig: B.6). This measurement file has to be the same as we will use with for algoriths (B.4 and B.5). XI B.3 Model Calibration with Dymola Figura B.6: Case of studio If measured data are given in some unit different than that used in the model (kW h and kW ), the scale column allows scaling of the measurements: var i able =data ∗scale +o f f set. In case the deviations of several variables shall be used to specify the criterion, the weight column allows the user to give them different weights (fig: B.7). Figura B.7: Details of measurements Integrator element allows the specification of a global simulation interval, this is important on load file when time is different to 0 as this case. Here it is studied the building during one year (fig: B.8). XII B.3 Model Calibration with Dymola Figura B.8: Time interval for simulation After validation (fig: B.9), an error of criterion of 4,7e6 has been detected and total energy difference for one week respect reference is 21.5%. Now, we change the task in Cases >Calibrate and Select parameters in Tuner parameters (fig: B.10). Figura B.9: Results of Validation for a week Parameter selected for the calibration are the three inner loads This three parameters corres- XIII B.3 Model Calibration with Dymola pond to estimated with norms DIN 18599 and SIA 2024 (tab: B.1) for this reason initial values are (110, 41, 15.9). Also, lower and upper bounds are from 0 until 250 for the two parameters people and machines, but bounds for ratio of lights is from 0 until 25. Bounds should permit to find optimal solution. Figura B.10: Tuner parameters After iterations Dymola returns a screen with results. Number of iterations (105), best result of parameter (180, 52.7, 0.2)and result of criterion (2,44e5). Which entails a total energy difference of 7% respected measure.(fig: B.11). Mainly Tuning affects the results decreasing the error by 94,8% as we can see comparing figures B.9 and B.11. XIV B.4 Matlab Optimization toolbox Figura B.11: Results of calibration for a week B.4. Matlab Optimization toolbox Optimization Toolbox provides widely used algorithms for standard and large-scale optimization. These algorithms solve constrained and unconstrained continuous and discrete problems. The toolbox includes functions for linear programming, quadratic programming, binary integer programming, nonlinear optimization, nonlinear least squares, systems of nonlinear equations, and multiobjective optimization [The MathWorks, 2003]. They can be used to find optimal solutions, perform tradeoff analyses, balance multiple design alternatives, and incorporate optimization methods into algorithms and models. Pattern search method is used with optimization toolbox (fig: B.12).The method is based on a search of the minimum value for objective function. Flow chart shows the method used by Pattern Search. It constructs a mesh which is then explored. Poll occurs when the search step was unable to obtain a point on the current mesh that decreased the incumbent value. If no decrease in cost is obtained on mesh points around the current iterate, then the distance between the mesh points is reduced (refine mesh), and the process is repeated. [Audet u. XV B.4 Matlab Optimization toolbox J. E. Dennis, 2000] Figura B.12: Basis of pattern search method B.4.1. Dymola to Matlab connection Connection between both programs is an important part of the project. The script "dymolaMçan execute every command that can be used in Dymola. Thus, the achieved actions for model calibration are: 1. Translate model into Dymola. 2. Set estimated parameters with Pattern Search. 3. Simulate the model. Also, dymload and dymget scripts load a result file from Dymola to Matlab and get the values of the searched parameter respectively. XVI B.4 Matlab Optimization toolbox Figura B.13: Dymola - Matlab conecction B.4.2. Algorithm and solver The optimization algorithm is developed to be more efficiency with regard time effort, accuracy and stability of results. However, this tool needs great effort to be implemented (See Matlab code in 8). Also, the only way to get efficiency on program is to set optimization toolbox options in Matlab with simulation tests (fig: B.14). Figura B.14: Options evaluation for algorithm in a simple one week study Furthermore, all specific cases studied in project can be adapted in an easy way with next steps: 1. Script modification .Give your measured data as a .mat file (for the project was not necessary to change). XVII B.4 Matlab Optimization toolbox .Set initial condition vector (x0) of parameters. .Setlowerand upper bounds(lb,ub) as wellas constrainsinsuch case (A,b,Aeq,beq). .Adjust tolerance of objective function results (Tolmesh : 1e−06). 2. Function modification .Give the parameter names to change their values in Dymola . .Change simulation characteristics (time, number of intervals, tolerance or solver). .Change criterion if it was necessary (|Total Heat −ydata|2). In the optimization function (fig: B.15), at each iteration either a comparison shows improvement (simple decrease) in the function value and the improving point becomes the new iterate or no decrease in value is found at any point in the pattern (fig: B.16). Figura B.15: Algorithm running with Optimization toolbox XVIII C.2 Study of fixed parameters in one thermal zone model C.2. Study of fixed parameters in one thermal zone model On the basis of simplified model (C.1), next cases will provide a first study to the methods with such the error of calibration is bigger than in complex cases. This is because number of parameters is lower than six zones model. C.2.1. One parameter: fixed heat flow The most simplified model correspond to substitute all the user model for a heat power load (fig: C.2). One zone model has 827 parameters versus 4581 parameters that reference model with six zones has. Figura C.2: J1615 people heat contribution XXV C.2 Study of fixed parameters in one thermal zone model After calibration at easy model, fixed heat flow and error calculated by optimization toolbox and mathematical method have a good agreement between them. Also, mathematical method is the fastest and needs less iterations to converge the solution (tab C.1). Cuadro C.1: Results of one zone calibration Time for simul. 03.01-10.01 31.01-07.02 28.02-07.03 05.12-12.12 1 year FixedFlow (W) Dymola 6299 9422 6169 6469 4346 Toolbox Matlab 11720 7624 5576 5240 6392 Math method 11698 7650 5553 5244 6396 Error 2 Dymola 5.79e5 9.29e5 5.92e5 1.45e6 6.69e9 Toolbox Matlab 4.62e5 5.91e5 4.62e5 1.09e6 2.18e10 Math method 4.69e5 5.91e5 4.62e5 1.09e6 2.18e10 Iterations Dymola 25 29 17 13 25 Toolbox Matlab 13 14 11 11 15 Math method 5 8 4 4 5 Time (s) Dymola 186 188 136 100 8312 Toolbox Matlab 76 63 65 81 560 Math method 16 30 18 14 336 C.2.2. Three parameters: people, machines and light ratio In next case, the influence of three parameter on heating consumption is studied. In fact, these parameters are multiply by the percentage of building presence at each hour giving in UserProfileOffice.txt (fig: C.4). The profile is repeated every day from Monday until Friday and null at weekend. Normally, more presence of people on building is related with more machines and light contribution. Nevertheless, we consider heating demand as fixed input and then increasing number of person reduces power due to machines and or lights (tab: C.2). In spite of the lower error of Dymola calibration toolbox, it requires enough time and iterations. Dymola results are interesant because all three parameters are lower than obtained with Matlab, this is illogical for one reason, criterion for all method is the sum of squared difference. Now, mathematical method is also as accurate as optimization method. However, it is necessary more time for parameter estimation than before. The next graphic is the result of apply the three estimate parameters with the toolbox of Matlab along a weeek (fig: C.3). First graphic represents a instantaneous heating power, whi- XXVI C.2 Study of fixed parameters in one thermal zone model Cuadro C.2: Calibration of one zone with three parameters for one year Dymola Toolbox Matlab Math method Nr People 104 113 109 Machines 18 33 40 Light Ratio (W/m2) 8.0 14.2 15.8 Error 2(kW h) 6.57e9 2.07e10 2.07e10 Iterations 31 14 5 Time (s) 6520 1498 776 le second one represents dynamic number of person and machines estimate. Last graphic shows the power generates by people, machines and lights. Figura C.3: Dynamic results of heating power, parameter values and inner loads heat power C.2.3. Eleven parameters: user profile people Office building is opened between 07 a.m. until 6 p.m. and occupation along the day change constantly, that is the reason why hourly presence of persons could indicate a significant decrease in objective function. This kind of parameters can be modified by both developed XXVII C.3 Building model: J1615 as six thermal zones programs of Matlab. Nevertheless, calibration function of Dymola is not allow to estimate them (fig: C.4). Figura C.4: User profile for one zone Optimization toolbox of Matlab has an expensive waste of time calibrating all parameter. However, mathematical method fit elf parameters with a similar error but much faster. Difference between setting of parameters depending of used method is for bound estimation. In case of toolbox Matlab, the function of optimization allows adjustment by upper and lower bounds. Nevertheless it is not the same for mathematical method where the step size of parameters depend on an equation that determine step and direction without bounds. Cuadro C.3: Calibration for eleven parameters Hour 7h 8h 9h 10h 11h 12h 13h 14h 15h 16h 17h Default% ocup 0.2 0.4 0.6 0.8 0.8 0.4 0.6 0.8 0.8 0.4 0.2 Toolbox Matlab 0.278 0.398 0.590 0.780 0.780 0.390 0.59 0.827 0.782 0.398 0.205 math method 0.421 0.426 0.410 0.450 0.435 0.418 0.427 0.401 0.423 0.460 0.428 Error 2(kW h) Iterations Time (sec) Toolbox Matlab 2.06e10 65 13750 Math method 2.07e10 6 2821 C.3. Building model: J1615 as six thermal zones In this model are diferentiated six zones. Each zone behaves in many different ways. However, total heating energy is calculate as addition of individual demands. For creation of this model all contribution of inner loads were held (fig: C.5). Complex model implies an increase in number of parameters making up the model; exactly to 4581 parameters. That means, a model with almost five times more equations due it has five rooms mores than model of one zone. XXVIII C.3 Building model: J1615 as six thermal zones (a) People contribution (b) Machines contribution (c) Light contribution (d) Inner loads together Figura C.5: Total inner loads energy Figura C.6: J1615 model as six zones XXIX C.4 Study of fixed parameters in six thermal zone model C.4. Study of fixed parameters in six thermal zone model C.4.1. Three parameters: people, machines and Light ratio After applying the methods on this case, a varied parameter estimation is set for each method. Thus shows the huge possibility of parameter combinations that converge in a minimum solution. Mathematical method is not so fast than before cases despite squared difference is reduced significantly respect to optimization algorithm of Matlab. However, Dymola has the lower error obtained with almost four times more requirements of time and 9 iterations more. Cuadro C.4: Calibration of six zone with three parameters for one year Dymola Toolbox Matlab Math method Nr People 26 222 117 Machines 102 92 268 Light Ratio (W/m2) 25.0 14.5 13.7 Error 2(kW h) 1.77e9 1.37e10 1.15e10 Iterations 23 14 14 Time (s) 11402 3420 3520 C.4.2. Six zones with fixed heat flow Similar case is carried out as studied for one parameter in one zone. The difference is only six spaces instead one. Main propose of case is the stability of methods, specially the consistence of calibration function with Dymola. Calibration with Dymola was impossible due to a fail in parameter estimation (fig: C.7). That means, the maximal number of parameter that Dymola is capable to perform in our model is 3. XXX C.4 Study of fixed parameters in six thermal zone model Figura C.7: Failure in calibration with Dymola for six parameters Matlab algorithms submit stability in simulation. Thus, no parameter estimation is necessary because the goal is to compare the three methods. XXXI D Results D.1. Sensibility analysis and weight of parameters Parameters on model have not a linear dependence between them, only mathematical model obtains linearise. Changing number of people in building to be estimated influences the output. However, number of machines and light ratio also influence the result in a similar way such that every parameter is not possible to be estimated individually. When a set of parameters have been tuned, the model is validated and the tuned parameters against other measured data to check that there is a good agreement between the simulation result and the new measurements. For a specific series of measured data it is possible to get good fits by increasing the model complexity and the number of tuned parameters. However, this does not guarantee that the result is that good for other operating conditions. The one zone model has a default values of 110 persons, 41 machines and 15.9 W /m2. After the procedure in (B.4) having into account the real heating energy data, fitting values with optimization algorithm are respectively 113.5 persons, 33 machines and 14.25 W /m2. The weight of parameters has been calculated with Dymola, the linear parameter combination should be fixed in such a way to avoid one parameter take hold of the others. The next relation has been estimated for one week with a Dymola function that is able to calculate dependence between parameters: NrPeopleMachines +0,8258·Nr People +8,5919·LightingPower This equation relates three parameters with each other. Then, no unique solution exists and all punts along the valley are possible solutions. (fig: D.1). The weight of parameters refer is 7.9% number of people, 9.6% number of machines and 82.5% ratio light. However, it is no possible to apply it for eleven parameter because this function is specific of Dymola that does not permit to estimate parameter of external file. XXXII D.1 Sensibility analysis and weight of parameters Figura D.1: Solution after sweep two parameter on one zone building In the studied case, the estimation of internal loads for an office building is more practicable than for a laboratory, as standard values are available for offices but not for laboratories. The differences between the simulation and measurement of the heat consumption are clearly correlated to the accuracy of the internal loads estimation. The manual ventilation and heater set temperature are also crucial factors (fig: B.4). After one week simulation, the main heat power flow depending of the parameters that variate with the time are represented in (fig: D.2), as we can see the influence of the air exchange and infiltration makes the sensor temperature on room drops bellow the temperature 20,5oC during office hours. As consequence of this decreasing heating equipment turns on until temperature goes beyond 20,5oC. During no office hours control temperature between 18 p.m. and 07 a.m. of next day is fixed to 16oC, but never goes beyond 19,5oC. XXXIII D.2 Comparison of methods and simulations Figura D.2: Variable parameters D.2. Comparison of methods and simulations Through the simulation of different cases along the project we can estimate behaviour of three used methods depending of results in section C. The analysis of simulation results help to determine the best method in every case depending of the interest of program user. 1. Case 1 (1 parameter in onezone) (C.2.1). .Parameter estimation for both Matlab methods are almost the same for each time of simulation. .The major assessment about case 1 is referred to time, twelve times faster than Dymola method. 2. Case 2 (3 parameters in onezone) (C.2.2). XXXIV Annex A: Code for optimization algorithm script 75 end 76 end 77 end 78 end 79 end 80 81 % dymolaM−command can execute every command that can be used in Dymola. 82 % Translate model J1615_onezone3Parameters . 83 res=dymolaM( ’ translateModel (" Optimization2 . J1615_onezone3Parameters ") ’ ) ; 84 85 % Parameters to s et for optimization . 86 TotalNrPeople=x(1) 87 NrPeopleMachines=x(2) 88 RatioLights=x (3) 89 90 % Perform of parameter values in Dymola . 91 dymolaM([ ’House.zoneParam. NrPeople= ’ ,num2str( TotalNrPeople ) ] ) ; 92 dymolaM( [ ’House. zoneParam. NrPeopleMachines= ’ ,num2str(NrPeopleMachines) ] ) ; 93 dymolaM([ ’House.zoneParam. LightingPower= ’ ,num2str( RatioLights ) ]) ; 94 95 % Simulation of building . 96 if simulate_week==1 97 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=0, stopTime =604800, numberOfIntervals=0, outputInterval =3600, tolerance =1, resultFile =" J1615_onezone3Parameters") ’) ; %Simulation in Dymola with the value . k =x ( 1) 98 else i f simulate_week==2 99 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime =2419200, stopTime=3024000, numberOfIntervals=0, outputInterval=3600, resultFile ="J1615_onezone3Parameters") ’ ) ; 100 else i f simulate_week==3 101 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime =4838400, stopTime=5443200, numberOfIntervals=0, outputInterval=3600, tolerance =1 , resultFile ="J1615_onezone3Parameters ") ’ ) ; 102 else i f simulate_week==4 103 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=29030400, stopTime=29635200, numberOfIntervals=0, outputInterval =3600, tolerance =1, resultFile =" J1615_onezone3Parameters") ’) ; 104 else i f simulate_week==52 105 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=0, stopTime=31536000, numberOfIntervals=0, XLI Annex A: Code for optimization algorithm script outputInterval =3600, tolerance =1, resultFile =" J1615_onezone3Parameters") ’) ; 106 end 107 end 108 end 109 end 110 end 111 112 % Loads data from r e s u l t f i l e .mat in a structure ’d ’ . 113 d=dymload( [onezone , ’/optimierung/J1615_onezone3Parameters .mat ’ ]) ; 114 115 % Extract from structure ’d ’ values of simulation . 116 TotalHeat= dymget(d, ’House. idealHeaterCooler_var_setpoint . heatMeter .q_kwh’ ) ; 117 TotalHeat= TotalHeat ( 1:( length(TotalHeat)−1)) ; 118 119 F=sum(( TotalHeat−ydata ) .^2) XLII Annex B: Code for mathematical algorithm function 1%%%%%%%%FUNCTION FOR MATHEMATICAL ALGORITHM %%%%%%%% 2function [ x_star ] = levenberg_marquardt11 3tic ; 4 5% Path of Dymola m f i l e s . 6DymolaPath = ’C:/ Program Files /Dymola 2013 FD01/Mfiles ’ ; 7 8% Path of workspace . 9MatlabPath = ’D:/mfu−lja /workspaces ’ ; 10 addpath(genpath(DymolaPath) ) ; 11 addpath(MatlabPath) ; 12 addpath ([ MatlabPath , ’/onezone ’ ] ) ; 13 addpath(genpath ([ MatlabPath , ’ /onezone/optimierung ’ ]) ) ; 14 15 % Change directory in Dymola. 16 onezone = [MatlabPath , ’/onezone ’ ] ; 17 dymolaM([ ’cd(" ’ ,onezone , ’/optimierung") ’ ]) ; 18 19 load ([ MatlabPath , ’/dataJ .mat ’ ]) ; % load 1615dinamicTable . csv data . 20 ydata=( dataJ (49:8809 ,5) /1000) −16540; % measurement data for one year . 21 22 % solve nonlinear lea s t squares problem using Levenberg−Marquardt 23 % algorithm . 24 % i n i t i a l guess . 25 x = [0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5]; 26 27 % increment in x for f i n i t e d if fe re nce approximation . 28 delta_x = [0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01;0.01]; 29 30 %great lambda means short stepsize . 31 lambda = 1; 32 33 %counter 34 counter = 0; XLIII Annex B: Code for mathematical algorithm function 35 progress = 1; 36 37 while progress > 1e−6 38 %counter 39 counter = counter+1; 40 41 % Simulation of Total Heat at vector of parameters ’ x ’ . 42 % Call to dymola here with value of x 43 userprofil; % SCRIPT for editting the column corresponding to UserProfilesOffice . txt . 44 45 % Model translation 46 res=dymolaM( ’ translateModel (" Optimization2 . J1615_onezone3Parameters ") ’ ) ; 47 48 % Model simulation 49 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=0, stopTime=31536000, numberOfIntervals=0, outputInterval=3600, tolerance =1, resultFile ="J1615_onezone3Parameters") ’ ) ; 50 51 d=dymload( [onezone , ’/optimierung/J1615_onezone3Parameters .mat ’ ]) ; 52 TotalHeat= dymget(d, ’House. idealHeaterCooler_var_setpoint . heatMeter .q_kwh’ ) ; 53 TH= TotalHeat (1 :( length(TotalHeat)−1)) ; %[8761x1 ] 54 55 % Simulation of Total Heat at x+dx for each parameter . 56 for k = 1: length(x) 57 58 % x plus dx for parameter k . 59 x_plus_dx = x(k) + delta_x (k) ; 60 61 % Simulation of Total Heat at x plus dx . 62 x(k)=x_plus_dx ; 63 userprofil; 64 x(k)= x_plus_dx −delta_x (k) ; 65 res=dymolaM( ’ translateModel (" Optimization2 . J1615_onezone3Parameters ") ’ ) ; 66 dymolaM( ’simulateModel (" Optimization2 . J1615_onezone3Parameters " , startTime=0, stopTime=31536000, numberOfIntervals=0, outputInterval=3600, tolerance =1, resultFile ="J1615_onezone3Parameters") ’ ) ; 67 68 d=dymload( [onezone , ’/optimierung/J1615_onezone3Parameters .mat ’ ]) ; 69 TotalHeat= dymget(d, ’House. idealHeaterCooler_var_setpoint . heatMeter .q_kwh’ ) ; 70 TH_dx( : , k)= TotalHeat (1 :( length(TotalHeat)−1)) ; % TH_dx [8761x11 ] 71 XLIV Annex B: Code for mathematical algorithm function 72 F=sum(abs(TH_dx( : , k)−ydata ) .^2) %criterion 73 end 74 75 % Right hand side 76 RHS = ydata −TH; %[8761x1 ] 77 78 % gradient 79 for k = 1: length(x) 80 81 % Partial derivative in k−direction . 82 P(: , k) = (TH_dx( : , k) −TH) ./ delta_x (k) ; % P [8761x11 ] each column i s divided by the column delta_x 83 end 84 85 % stepsize 86 delta_x = (P’*P + lambda*diag (diag(P’*P) ) ) \(P’*RHS) ; % [11x1 ] = { ([ 11 x8761 ][8760x11 ] ) + [11x11 ] } \ [11x1 ] 87 88 % update 89 x = x + delta_x ; 90 91 % progress 92 progress = ( delta_x ’*delta_x ) /(x ’*x) ; % [1 x11 ][1 1 x1 ] / [1 x11 ][1 1 x1 ] 93 94 if mod( counter ,10) == 0 95 disp ([ ’Counter : ’ ,num2str( counter ) ]) ; 96 lambda = 0.1*lambda; 97 end 98 end 99 100 % optimal solution 101 x_star = x 102 103 % number of iterations 104 disp ([ ’Number of iterations : ’ ,num2str( counter ) ]) ; 105 106 tiempo = toc; 107 fprintf( ’The process took %d seconds ’ , tiempo) ; 108 109 110 %%%%%%%%SCRIPT USERPROFIL %%%%%%%% XLV Annex B: Code for mathematical algorithm function 111 f i l e = [ MatlabPath , ’ /onezone/optimierung/Tables/J1615 ’ ] ; 112 filename = [ file , ’/ UserProfilesOffice . txt ’ ] ; 113 114 Time=[0;3540;3600;7140;7200; . . . ;597600;601140;601200;604740;]; 115 People= [ 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; x (1) ; x (1) ; x(2) ; x (2) ; x (3) ; x (3) ; x (4) ; x (4) ; x (5) ; x (5) ; x(6) ; x(6) ; x (7) ; x (7) ; x(8) ; x(8) ; x (9) ; x(9) ; x(10) ; x(10) ; x(11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; 0; 0 ; 0 ; x (1) ; x (1) ; x (2) ; x (2) ; x (3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x(10) ; x(10) ; x(11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ;0 ; 0 ; 0; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; x (1) ; x (1) ; x (2) ; x (2) ; x(3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x(10) ; x(10) ; x (11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ;0 ; 0 ; 0; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; x (1) ; x (1) ; x (2) ; x(2) ; x(3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x(10) ; x (10) ; x(11) ; x(11) ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0; 0 ; 0 ;0 ; 0 ;0 ; 0 ; 0; 0 ; 0 ;0 ; 0 ; x(1) ; x (1) ; x(2) ; x(2) ; x(3) ; x(3) ; x(4) ; x(4) ; x(5) ; x(5) ; x(6) ; x(6) ; x(7) ; x(7) ; x(8) ; x(8) ; x(9) ; x(9) ; x (10) ; x(10) ; x(11) ; x(11) ; 0; 0 ;0 ;0 ;0 ; . . . ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; ] ; 116 Machines = [0. 1; 0.1 ;0 .1; 0. 1; . . . ; 0 . 1 ; 0 . 1 ; 0 . 1 ; 0 . 1 ; ] ; 117 Light1=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1;1;0.3; ... ;0;0;0;0;]; 118 Light2=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.3;0.3;0.3; ... ;0;0;0;0;]; 119 A=[Time, People , Machines , Light1 , Light2 ] ; 120 121 122 fileID = fopen ([ file , ’/ UserProfilesOffice . txt ’ ] , ’wt ’ ) ; 123 fprintf( fileID , ’#1\n ’ ) ; 124 fprintf( fileID , ’double UserProfilesOffice (336 , 5) \n ’ ) ; 125 126 127 for i =1:length(A) 128 fprintf( fileID , ’ %f \ t %f \ t %f \ t %f \ t %f \n ’ ,A( i , : ) ) ; 129 end 130 131 fclose ( fileID ) XLVI Bibliografía [AB 2012] AB, Dassault S.: Dymola: Dynamic Modelling Laboratory User Manual. Volume 2, 2012 [Audet u. J. E. Dennis 2000] AUDET, Charles ; J. E. DENNIS, JR: Pattern Search algorithms for mixed variable programming. (2000), S. 573–594 [Bernard P. Zeigler 2000] BERNARD P. ZEIGLER, Tag Gon K. Herbert Praehofer P. Herbert Praehofer: Theory of Modeling and Simulation: Integrating Discrete Event and Continuous Complex Dynamic System. (2000) [Bertsekas 1999] BERTSEKAS, Dimitri P.: Nonlinear Programming. 1999 [H. Elmqvist 2005] H. ELMQVIST, S.E. Mattsson D. Brueck C. Schweiger D. Joos M. O. H. Ollson O. H. Ollson: Optimization for Design and Parameter Estimation. (2005), S. 255–266 [M. Lauster 2012] M. LAUSTER, M. Fuchs R. Streblow D. M. J. Teichmann T. J. Teichmann: Dynamic building model for city quartier simulation. (2012) [Stefan Finsterle 2010] STEFAN FINSTERLE, Michael B. K.: A truncated Levenberg-Marquardt algorithm for the calibration of highly parameterized nonlinear models. (2010), S. 731–738 [The MathWorks 2003] THE MATHWORKS, Inc.: Optimization Toolbox Users Guide of Matlab, 2003 [Wetter u. Polak 2003] WETTER, Michael ; POLAK, Elijah: A convergence optimization method using pattern search algorithms with adaptative precision simulation. (2003), S. 1393–1400 [Wetter u. Wright 2003] WETTER, Michael ; WRIGHT, Jonathan: Comparison of a generalized pattern search and a genetic algorithm optimization method. (2003), S. 1400–1408 XLVII