Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1)

Equation Chapter 1 Section 1 Trabajo Fin de Grado Ingeniería Electrónica, Robótica y Mecatrónica Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) Autor: Sofía Tarín Roldán Tutor: Manuel Gil Ortega Linares Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 iii Trabajo Fin de Grado Ingeniería Electrónica, Robótica y Mecatrónica Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) Autor: Sofía Tarín Roldán Tutor: Manuel Gil Ortega Linares Catedrático de Universidad Dpto. Ingeniería Automática y Sistemas Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Trabajo Fin de Grado: Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) Autor: Sofía Tarín Roldán Tutor: Manuel Gil Ortega Linares El tribunal nombrado para juzgar el Trabajo arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A todos los que nunca dejan de creer en sus sueños y luchan por ellos cada día latido a latido ix Agradecimientos Gracias a Manuel Gil, mi tutor durante este trabajo. Desde la primera clase consiguió transmitirme algo que creo que es imprescindible en el aprendizaje: la curiosidad. Y también desde la primera tutoría en la que le hablé de este trabajo para proponérselo me brindó todo su apoyo, un apoyo que se ha mantenido a lo largo del mismo hasta acabarlo. Sin su gran ayuda este trabajo no sería posible. Gracias a mi familia, a mis padres, a mis tías y a mis abuelos, por todo lo que me han dado en la vida, y por todo lo que me han enseñado. Gracias por todo el esfuerzo que cada uno de ellos han hecho, de una manera u otra, para que a día de hoy yo pueda estar escribiendo estas líneas y acabando esta etapa de mi vida. Gracias por el apoyo y el cariño incondicional. Gracias a Rafa, el mejor regalo que me ha podido dar esta aventura en la universidad. Gracias por estar siempre a mi lado para sacarme una sonrisa cuando más lo necesito, por querer compartir conmigo esta montaña rusa que es la vida y por no soltarme nunca de la mano. Estoy muy orgullosa de lo que hemos construido y seguimos construyendo día tras día juntos. Gracias a todos mis amigos, especialmente a los que me ha dado el grado, tanto dentro de clase como fuera de ella. Me habéis acompañado y apoyado durante estos 4 años siendo mi segunda familia aquí en Sevilla, y espero que así sea toda la vida. En especial a Pablo, que tantos y tantos momentos hemos vivido juntos desde el comienzo de esta etapa, algunos mejores y otros peores, pero sabemos que siempre estamos ahí cerca para apoyarnos el uno al otro. Gracias a los profesores y maestros que he tenido a lo largo de mi vida por todo lo que me han podido enseñar. Cada uno de ellos y ellas han puesto su granito de arena para que hoy sea la persona que soy. En especial, gracias a Francisco Javier León, un verdadero apasionado por el mundo de la tecnología que me transmitió su sana locura, especialmente por los robots, y con ello que hace cuatro años yo eligiera este grado y decidiera dedicarme a esto. Por último, gracias a todas las personas que, de una forma u otra, me han apoyado y han hecho posible este trabajo y lo que ello ha conllevado. Sofía Tarín Roldán Sevilla, 2021 ÍNDICE DE TABLAS Tabla 2-1. Parámetros del Modelo de Bergman 9 Tabla 2-2. Parámetros de la ecuación de perturbaciones del Modelo de Bergman 9 Tabla 2-3. Protocolo de perturbaciones de comida 12 Tabla 2-4. Protocolo de perturbaciones de ejercicio físico 12 Tabla 5-1. Comparación cuantitativa de los controladores 44 xvii ÍNDICE DE FIGURAS Figura 1-1. Poster informativo sobre la diabetes 1 Figura 1-2. Medición de la glucemia con un glucómetro 3 Figura 1-3. Medición de la glucemia con un dispositio de monitorización continua 4 Figura 1-4. Plumas de distintos tipos de insulina y aguja desechable 5 Figura 1-5. Ejemplo de bomba de insulina 6 Figura 2-1. Esuqma del Modelo de Bergman 8 Figura 2-2. Modelo de Bergman en Simulink 11 Figura 2-3. Interior del bloque "Modelo_Bergman" 11 Figura 3-1. Bode Gm(s) 16 Figura 3-2. Bode Gm(s)*C(s) sin ajuste de Kp 17 Figura 3-3. Bode de Gm(s)*C(s) final 18 Figura 3-4. Modelo con controlador PID continuo en Simulink 19 Figura 3-5. Comportamiento del modelo de Bergman sin perturbaciones 19 Figura 3-6. Comportamiento del modelo de Bergman con perturbaciones 20 Figura 3-7. Control PID continuo sin perturbaciones 21 Figura 3-8. Contro PID continuo con perturbaciones 21 Figura 3-9. Modelo con controlador PID discreto en Simulink 23 Figura 3-10. Bloque del controlador PID en tiempo discreto en Simulink 23 Figura 3-11. PID discreto sin perturbaciones 24 Figura 3-12. PID discreto con perturbaciones 24 Figura 3-13. PID con algoritmo anti-windup 25 Figura 3-14. PID con tiempo de muestreo de 3 minutos 27 Figura 3-15. PID con tiempo de muestreo de 5 minutos 27 Figura 3-16. PID con tiempo de muestreo de 10 minutos 28 Figura 3-17. Bloques del control PID con inyección de insulina por lotes 29 Figura 3-18. PID con inyección de insulina por lotes 29 Figura 3-19. Zoom del PID con inyección de insulina por lotes 30 Figura 3-20. Tiempo de inyección de insulina 30 Figura 3-21. Modelo del controlador PID Feedforward 31 Figura 3-22. Esquema para la obtención de Cff 32 Figura 3-23. Bloque del control PID Feedforward 33 Figura 3-24. Controlador PID Feedforward con inyección de insulina por lotes 33 xix Figura 4-1. Modelo del controlador MPC 37 Figura 4-2. Bloque del control MPC 37 Figura 4-3.Estructura del controlador MPC 38 Figura 4-4.Barra de herramientas para el diseño del controlador MPC 38 Figura 4-5. Controlador MPC HP: 300 s; Hc: 10 s; wu: 0.1; wg: 1 38 Figura 4-6. Controlador MPC HP: 300 s; Hc: 100 s; wu: 0.1; wg: 1 39 Figura 4-7. Controlador MPC HP: 300 s; Hc: 100 s; wu: 0.1; wg: 10 39 Figura 4-8. Controlador MPC HP: 500 s; Hc: 50 s; wu: 0.1; wg: 50 40 Figura 4-9. Controlador MPC HP: 300 s; Hc: 50 s; wu: 0.1; wg: 50 40 1 1 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) 1 INTRODUCCIÓN a diabetes mellitus (DM) es una enfermedad crónica del sistema metabólico que se produce cuando el páncreas deja de producir insulina, o bien la insulina que segrega tiene una calidad inferior y el organismo no es capaz de usarla de una manera eficaz. La insulina es una hormona cuya función es regular la glucosa (azúcar) en sangre. El resultado principal de la diabetes son las hiperglucemias (aumento de la glucosa en la sangre por encima de los valores normales), que a largo plazo pueden causar daños muy importantes en diferentes órganos, así como provocar el desarrollo de otras enfermedades. Según la OMS (Organización Mundial de la Salud), en 2014 el número total de personas en el mundo que sufrían diabetes era de 422 millones, con un ascenso de 314 millones desde el año 1980. La prevalencia en adultos se sitúa entorno al 8.5%, y alrededor de 1.6 millones de muertes en 2016 fueron causadas a la diabetes. Asimismo, la diabetes es una de las principales causas de ceguera, infartos y amputación de los miembros inferiores. [1] L “Cuando el corazón y la mente están unidos, todo es posible”. - Diego Pablo Simeone - Figura 1-1. Poster informativo sobre la diabetes Introducción 2 2 Una persona no diabética mantiene siempre sus niveles de glucosa en un rango de 70 a 110 mg/dL, aunque lo que haya comido tenga una alta cantidad de azúcares o hidratos de carbono. Dichos valores de glucosa forman los llamados “valores normales”. Este mantenimiento constante dentro del rango es gracias a la secreción rápida de insulina. Cuando una persona ingiere alimentos, los azúcares presentes en los mismos se absorben y pasan a la sangre. Los niveles de glucosa se elevan, pero inmediatamente las células que se encargan de producir insulina detectan dicho aumento de glucosa y segregan de forma rápida esta hormona, haciendo que la glucosa entre en los tejidos y músculos donde posteriormente se metabolizará para producir energía, y de esta forma se disminuye los niveles de glucosa en sangre. Sin embargo en una persona con diabetes, como ya hemos dicho antes, la producción de insulina es muy baja o incluso nula, por lo que cuando la glucosa aumenta en sangre no se produce la secreción instantánea de insulina y los niveles se mantienen altos, por encima de los valores altos, imposibilitando la conversión de dicha glucosa en energía para las células. [2] 1.1 Tipos de diabetes Se puede diferenciar 3 tipos de diabetes: - Diabetes tipo 1: también denominada diabetes insulinodependiente o juvenil. Su principal característica es la producción muy deficiente de insulina, debida a la destrucción de tipo autoinmune de las células beta del páncreas. No se conoce a día de hoy la causa que conlleva a dicha destrucción autoinmune y, por tanto, a la diabetes tipo 1 (DM1). Se suele desarrollar a una edad temprana, sobre todo en niños y jóvenes, y sus síntomas se conocen como las 3 P: polidipsia (excesiva sed), poliuria (excesivas ganas de orinar) y polifagia (hambre constante). Además, también conlleva una importante pérdida de peso, cansancio extremo y trastornos visuales, y en caso más extremos donde la diabetes no se diagnostica a tiempo, puede llegar a causar cetoacidosis (ocurre cuando el cuerpo descompone la grasa demasiado rápido y el hígado la convierte en cetonas, volviendo a la sangre ácida), coma diabético o incluso la muerte. Una vez se diagnostica, se pone en tratamiento con administración externa de insulina para controlar la glucosa e intentar que los niveles sean lo más cercanos posible a los valores normales. Dicho tratamiento es de por vida ya que no existe aún cura para la diabetes tipo 1, de ahí el nombre de insulinodependiente. - Diabetes tipo 2: también conocida como no insulinodependiente o de edad adulta. Esta diabetes se caracteriza por el uso poco eficaz de la insulina por parte del organismo, debido principalmente a que su calidad es peor. Las causas de la diabetes tipo 2 (DM2) son el sobrepeso, el sedentarismo (muy poca actividad física) y una dieta desequilibrada y rica en azúcares y grasas; de ahí a que se suela dar con mayor prevalencia en personas con una edad avanzada, aunque cada vez más se da en personas más jóvenes por malos hábitos de salud. Los síntomas característicos de la DM2 son parecidos a los de la de tipo 1, si bien normalmente no son tan intensos, lo cual hace que sea más difícil un detección precoz y se pueden pasar años con la enfermedad sin que esta haya sido diagnosticada. Por este motivo, cuando se detecta, suele haber algún tipo de complicación más. El tratamiento dependerá de lo avanzada que esté la diabetes. Si está en una fase inicial y los niveles de glucosa son ligeramente mayores a los normales, bastará con una dieta saludable, baja en azúcares y grasas, y actividad física. También se recomendará determinados hábitos saludables como no fumar o no tomar alcohol, para evitar complicaciones. Podrá tratarse con medicamentos orales para controlar la glucemia. Si, por el contrario, la diabetes está en una fase más avanzada, donde los niveles de glucosa son demasiado altos, se tratará con inyección externa de insulina (al igual que en la DM1), aunque si los niveles se van corrigiendo con el tiempo se podrá dejar este tratamiento y aplicar el mencionado anteriormente para fases más leves. 3 3 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) - Diabetes Gestacional: se caracteriza principalmente por hiperglucemias que se dan en algunas mujeres embarazadas. Aunque no se llegan a tener niveles de glucosa tan altos como para diagnosticar diabetes, sí que estos superan ligeramente los valores normales y recomendados. Sin embargo, este tipo de diabetes es igualmente peligrosa, pues hace que tanto la madre como el hijo puedan sufrir más complicaciones durante el embarazo y el parto. El diagnóstico de este tipo de diabetes se realiza mediante pruebas durante el embarazo, y para evitarla se recomienda un dieta equilibrada y ejercicio físico. [1] 1.2 Monitorización de la glucemia Actualmente existen dos tipos de dispositivos con los cuales podemos realizar una monitorización de la glucosa en sangre o glucemia: - Glucómetros: siendo el más usado en las últimas décadas, este dispositivo es un instrumento de medida de la glucemia capilar de forma inmediata. Su funcionamiento consiste en insertar una tira reactiva en el glucómetro, pinchar la yema del dedo con una aguja, y poner una pequeña gota de sangre en la tira reactiva. Pasados unos 5 segundos, obtendremos en la pantalla del glucómetro el valor de la glucemia capilar. Si bien es relativamente rápido, la principal y obvia desventaja de este método “tradicional” es lo invasivo que resulta, ya que puesto que una persona diabética necesita controlar variar veces al día sus niveles de glucosa en sangre, esto se traduce en una gran cantidad de pinchazos. Por otra parte, nos da datos de medidas puntuales y discontinuas, sin poder saber cómo ha variado la glucemia entre una medida y otra. - Sistema de monitorización continua de glucosa: son los nuevos dispositivos que, cada vez con mayor frecuencia, se usan a día de hoy. Se diferencian de los glucómetros porque, como su nombre indica, miden de forma continua la glucosa en sangre, cada minuto aproximadamente, lo cual es mucho más beneficioso. Estos sistemas se componen de un sensor con un filamento flexible el cual se inserta debajo de la piel, un monitor para observar los valores medidos (este monitor puede ser el propio teléfono móvil del paciente con una aplicación específica para esta tarea), y un transmisor que permite comunicarse al sensor con el monitor. Estos dispositivos miden el nivel de glucosa en el líquido intersticial (el espacio entre las células), que si bien difiere ligeramente del valor capilar, al final resulta ser prácticamente idéntico. Además, los sistemas de monitorización continua de glucosa suelen incluir alarmas que avisan cuando detectan que la glucosa supera, bien por arriba o por abajo, el rango establecido como deseado, disminuyendo así la probabilidad de tener hiperglucemias e hipoglucemias severas. [3] Figura 1-2. Medición de la glucemia con un glucómetro Introducción 4 4 1.3 Tratamiento con insulina en la Diabetes tipo 1 Puesto que para este trabajo me centraré en el control específicamente de la diabetes tipo 1, realizaré a continuación un breve repaso de los distintos tipos de tratamientos o terapias que existen para este tipo de diabetes. Como ya he mencionado antes, la diabetes tipo 1 se caracteriza por la nula secreción de insulina en el páncreas y, por consecuencia, es imprescindible para el paciente la administración exógena de insulina. Este procedimiento se podrá llevar a cabo de dos formas: mediante inyecciones de jeringuillas o plumas de insulina, o bien mediante una bomba de insulina. Para entender mejor estos dos tipos de administración, primero explicaré los tipos de insulina que existen. Principalmente, podemos distinguir dos tipos de esta hormona: la basal, que es la que genera el páncreas en pequeñas cantidades de forma continuada a lo largo del día, y los llamados bolos, que se corresponde con la insulina que se genera en mayor cantidad cuando la persona come y hay que contrarrestar la cantidad de glucosa ingerida. 1.3.1 Formas de administración de la insulina En las personas diabéticas se usan dos tipos de insulina exógena que imitan el comportamiento de las dos mencionadas anteriormente: la insulina de acción prolongada y la de acción rápida o ultrarrápida. - Insulina de acción prolongada: este tipo de insulinas tienen una duración aproximada de 24 horas. El comienzo de su efecto es alrededor de 2 horas después de la inyección y su curva de acción o actividad es prácticamente plana, es decir, que hace el mismo efecto desde esas 2 horas después de inyectarse hasta que se cumplen las 24 horas, sin tener picos de acción. Se utiliza como sustitución de la insulina basal. - Insulina de acción rápida o ultrarrápida: la duración de estas insulinas son de aproximadamente 3 o 4 horas, empezando a hacer efecto a los 15 minutos después de la inyección. Su curva de actividad Figura 1-3. Medición de la glucemia con un dispositivo de monitorización continua 5 5 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) presenta un pico en la hora posterior a la administración, es decir, cuando hace su máximo efecto. Esta insulina sustituye a los llamados bolos, la insulina que genera el páncreas después de la ingesta de alimentos. Es por eso que la insulina de acción ultrarrápida se administra antes de cada comida (desayuno, almuerzo y cena, y meriendas si es que se hacen), así como para corregir la glucosa cuando esta se encuentra a valores demasiado altos después de comer, si no se ha calculado bien la cantidad de insulina administrada antes de comer. [3] 1.3.2 Formas de administración de la insulina Actualmente existen dos formas principales de administrar la insulina, tanto basal como prandial (bolos): - Plumas: también conocidos como “bolígrafos” de insulina por su forma, estas vienen cargadas con un cartucho de un tipo de insulina y sirven para varias inyecciones, que serán más o menos dependiendo de la cantidad de insulina que necesita el paciente. El pinchazo es fácil y preciso, y se realiza colocando una pequeña aguja desechable tras la administración de unos 5 mm en la punta. Una vez se carga la cantidad deseada, solo hay que pulsar un botón y la insulina saldrá de forma inmediata, en apenas unos segundos. Las zonas más frecuentes donde se administra la insulina con las plumas son la zona superior de los brazos, los muslos, la parte inferior del abdomen y los glúteos. - Bombas de insulina: son unos pequeños dispositivos programables que contienen insulina de acción rápida o ultrarrápida que se va administrando de forma continua durante todo el día al paciente gracias a un catéter o cánula que se inserta normalmente en la zona inferior del abdomen. Las bombas son portátiles, lo que permite que el paciente vaya de una parte a otra con ella sin que esta le moleste en la mayoría de sus tareas cotidianas. Normalmente las bombas funcionan en bucle abierto, es decir, no reciben información de los valores de glucemia, y por tanto es el propio paciente quien tiene que programarla en función de la comida que vaya ingerir o de los valores de glucemia que tenga para que se administre más o menos insulina. Figura 1-4. Plumas de distintos tipos de insulina y aguja desechable Modelo del sistema glucoregulatorio 12 12 Perturbaciones de Comida Hora de la ingesta Gramos de C6H12O6 8:00 55 14:00 90 21:00 45 Tabla 2-3. Protocolo de perturbaciones de comida Perturbaciones de Ejercicio Hora Tipo de ejercicio 0:00-18:00 Reposo/Ejercicio suave 18:00-19:00 Ejercicio intenso 19:00-0:00 Reposo/Ejercicio suave Tabla 2-4. Protocolo de perturbaciones de ejercicio físico Dentro de las funciones de Matlab descritas antes se define este protocolo dependiendo del tiempo de simulación. La simulación durará un total de 86400 segundos, los equivalentes de 1 día completo. Se ha elegido este protocolo por ser estándar, con horas habituales de desayuno, almuerzo y cena en un paciente medio que además realiza ejercicio intenso una hora al día. La referencia de glucosa, como se verá en el capítulo siguiente dedicado al control, se pondrá a 90 mg/dL, un valor adecuado y entorno al cual suelen encontrarse las personas no diabéticas. 2.3 Linealización del modelo Como he dicho en apartados anteriores, el modelo de Bergman es un modelo no lineal. Por lo tanto, es necesario linealizarlo previamente a diseñar el control, y para ello se debe obtener las ecuaciones de estado que definen el modelo, decidir los puntos de equilibrio y obtener las matrices correspondientes a la linealización para poder sacar la función de transferencia del modelo G(s). Se quiere linealizar el modelo entorno a los siguientes puntos de equilibrio: Geq=0 mmol/L, Ieq=0 mU/L, Xeq=0 min-1, Deq=0, Ejeq=0, y Ueq=16.667 mU/min. Se obtienen las siguientes matrices: A = [−𝑃10 −𝐺𝑏 0 −𝑛 0 0 𝑃3−𝑃2] B = [0 1 −𝐺𝑏 1 𝑉𝑖0 0 0 0 0 ] C = [1 0 0] D = [0 0 0] (7) 13 13 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) Por lo tanto, y sabiendo que: 𝑥󰇗=𝐴𝑥+𝐵𝑢 𝑦=𝐶𝑥+𝐷𝑢 Siendo 𝑥, 𝑢 y 𝑦 los valores de los incrementos que se le suman a los puntos de equilibrio, tendríamos que las ecuaciones del modelo serían las siguientes: [𝐺󰇗𝐼󰇗 𝑋󰇗] = [−𝑃10 −𝐺𝑏 0 −𝑛 0 0 𝑃3−𝑃2][𝐺𝐼𝑋]+ [0 1 −𝐺𝑏 1 𝑉𝑖0 0 0 0 0 ][𝑈−𝑈𝑒𝑞 𝐷 𝐸𝑗 ] (8) 𝑦=[1 0 0][𝐺𝐼𝑋]+ [0 0 0][𝑈−𝑈𝑒𝑞 𝐷 𝐸𝑗 ] (9) Para calcular la función de transferencia del modelo, se tiene la siguiente ecuación: 𝐺(𝑠)=𝐶(𝑠𝐼−𝐴)−1𝐵+𝐷 (10) Dicha ecuación se calcula en el script de Matlab llamado “función_transferencia.m”, y se obtienen 3 funciones de transferencia: una para la señal de control u(t) y otras dos para cada una de las perturbaciones, D y Ej. En este caso, me quedo con la de la señal de control, que es la siguiente: 𝐺𝑚(𝑠)=−31.59 (15000𝑠+7)(2400𝑠+1)(648𝑠+1) (11) Ahora que ya está calculada la función de transferencia del modelo se puede empezar a diseñar el control en el siguiente capítulo. Modelo del sistema glucoregulatorio 14 14 15 15 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) 3 CONTROL PID n este capítulo se desarrollará un control PID básico, primero en tiempo continuo y después en tiempo discreto, al que se le irán añadiendo distintas modificaciones para mejorarlo. Se realizarán diversos experimentos para observar si su comportamiento es el deseado. 3.1 Controlador PID en tiempo continuo Primero se va a diseñar un controlador PID en tiempo continuo, es decir, con un tiempo de muestreo igual a 0. Este controlador no se puede usar en la realidad, ya que el algoritmo de control se implementa y calcula dentro de una computadora que tiene un cierto tiempo de muestreo, por pequeño que sea, y en donde por lo tanto tendremos un controlador en tiempo discreto. Sin embargo, este primer diseño en tiempo continuo servirá para definir especificaciones, así como la frecuencia de control y distintos parámetros del controlador. También se podrá observar al final si los experimentos con el controlador salen tal y como se esperaban y si es válido para el control de la glucosa. 3.1.1 Especificaciones definidas para el controlador Se va a diseñar el controlador usando el dominio de las frecuencias. La frecuencia de control (a la que trabajará el controlador) se puede elegir según dos criterios, del cual se escogerá el valor más restrictivo. El primer criterio será el de la frecuencia de Nyquist. Siguiendo dicho criterio se marcará que la frecuencia de control ωc sea al menos 10 veces superior que la frecuencia de Nyquist. Para calcular la frecuencia de Nyquist se tiene que: 𝜔𝑁=𝜋 𝑇𝑚 (12) Siendo Tm el tiempo de muestreo. Si bien en tiempo continuo el tiempo de muestreo como he dicho antes es 0, se va definir un tiempo de muestreo inicial de 1 minuto (60 segundos) en vistas al controlador en tiempo discreto. E “If you want to make the world a better place, take a look at yourself and then make a change”. - Michael Jackson - Control PID 16 16 Por lo tanto, obtenemos una frecuencia de Nyquist: 𝜔𝑁=𝜋 60= 0.05 𝑟𝑎𝑑/𝑠 (13) Sabiendo que la frecuencia de control debe ser al menos 10 veces mayor que la frecuencia de Nyquist: 𝜔𝑐=10∗0.05= 0.5 𝑟𝑎𝑑/𝑠 (14) El segundo criterio que se va a usar es el de limitación de fase, ya que no es recomendable que el controlador trabaje a frecuencias que se encuentren por encima de la frecuencia que marca los -180º de fase en el bode la función de transferencia del modelo. Para ver cuál es esta frecuencia se dibuja el bode de Gm(s) (de la ecuación (11)), que podemos ver en la figura 3-1. Podemos observar que la frecuencia en la que el sistema tiene una fase de -180º es ωc = -0.0012 rad/s. Puesto que se debe escoger la frecuencia de control más restrictiva de los dos criterios usados, me quedaré en este caso la frecuencia de control determinada por el criterio de limitación de fase, ya que es la más pequeña de las dos. Por lo tanto, mi frecuencia de control a partir de este punto es: 𝜔𝑐= 0.0012 𝑟𝑎𝑑/𝑠 (15) 3.1.2 Diseño del controlador PID en tiempo continuo Se va a diseñar el controlador, en el dominio de la frecuencia tal y como he mencionado antes, optando para el diseño del mismo que sea sin cancelación de dinámica. Esto es debido a que nos enfrentamos a un problema de regulación principalmente, y no de seguimiento. Se va a definir un valor constante de referencia de glucosa y el controlador se encargará de regular la respuesta del sistema de forma que ante las distintas perturbaciones que sufra bien por comida o bien por ejercicio esta se mantenga en el valor de referencia. Por lo tanto, el principal problema serán las perturbaciones y no los cambios a lo largo del tiempo de valores de referencia que el controlador deba ir siguiendo (problema de seguimiento). Puesto que los controladores diseñados con Figura 3-1. Bode de Gm(s) 17 17 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) cancelación de dinámica (cancelando los polos de la función de transferencia del modelo) trabajan peor ante las perturbaciones, se elige el diseño sin cancelación de dinámica, para obtener mejores resultados. El margen de fase actual que tenemos se calcula como: 𝑀𝑓𝑎𝑐𝑡=180+(−180)=0º (16) El margen de fase deseado se fija en 75º, para que posteriormente la fase del Gm(s)*C(s) sea aproximadamente -105º a la frecuencia de control. De esta forma, la fase que debe aportar el PID se calcula de la siguiente forma: 𝜑=𝑀𝑓𝑑𝑒𝑠−𝑀𝑓𝑎𝑐𝑡=75−0=75º (17) Se calcula la función de transferencia del controlador tipo PID: 𝐶(𝑠)=𝐾𝑝(𝜏1𝑠+1)2 𝜏1𝑠 (18) Para obtener el valor de 𝜏1se tiene la siguiente ecuación: 𝜑=2∗𝑎𝑟𝑐𝑡𝑔(𝜏1∗𝜔𝑐)−90 (19) De donde se despeja y se obtiene que: 𝜏1=6028.4 𝑠 (20) Sustituyendo en la ecuación (18) con el valor obtenido tenemos que la función de transferencia del controlador es: 𝐶(𝑠)=𝐾𝑝(6028.4𝑠+1)2 6028.4𝑠 (21) Se añade a esta función de transferencia un polo de alta frecuencia, que hace que la función del controlador sea bipropia y además ayuda a que la respuesta sea menos ruidosa. Para calcular este polo de alta frecuencia se tiene que: 1 𝜏𝑎𝑓≅10∗𝜔𝑐=0.012 𝑟𝑎𝑑/𝑠 (22) 𝜏𝑎𝑓=83.33 𝑠 (23) Con el polo de alta frecuencia añadido C(s) queda de la forma: 𝐶(𝑠)=𝐾𝑝(6028.4𝑠+1)2 6028.4𝑠 1 83.33𝑠+1 (24) Por tanto, solo faltaría hallar el valor de Kp para tener definida la función. Dicho valor se obtiene ajustando el Figura 3-2. Bode de Gm(s)*C(s) sin ajuste de Kp Control PID 18 18 bode para que la ganancia En la figura 3-2 se representa el bode resultante de multiplicar Gm(s)*C(s), es decir, la función de transferencia del modelo por la del controlador: Se observa que la fase a la frecuencia de control es -110º, 5 grados por debajo de lo que se esperaba según los cálculos de la ecuación (17). Estos 5 grados de pérdida de fase se deben al polo de alta frecuencia añadido, pero no afectan a la estabilidad del sistema en ningún caso y considero aceptable el resultado. Mirando el diagrama de ganancia, se ajusta la Kp de C(s) de forma que la ganancia sea 0 a la frecuencia de control, obteniendo un valor de Kp=-0.3385. Por lo tanto y finalmente se tiene que: 𝐶(𝑠)=−0.3385(6028.4𝑠+1)2 6028.4𝑠 1 83.33𝑠+1 (25) 3.1.3 Implementación del controlador PID en tiempo continuo en Simulink Una vez calculada C(s), se implementa el controlador en Simulink. Para ello se sustituye el bloque constante para U(t) que vemos en la figura 2-2 por un bloque de función de transferencia con la función de C(s) calculada. A este bloque le llega el error que resulta de restar la referencia de glucosa con la glucosa real, que es la salida del sistema. Puesto que al hablar de función de transferencia hablamos necesariamente de incrementos, la salida del controlador es un incremento al que habrá que sumarle el punto de equilibrio, que en este caso es Ueq definido en el capítulo 2.3, donde se linealizaba el modelo. Ya que en Simulink estoy trabajando en segundos, tenemos que Ueq es igual a 0.277 mU/s. Se aplica a la señal de control un bloque de saturación, tanto por debajo ya que no se pueden inyectar valores negativos de insulina (sería el equivalente a sacar insulina del cuerpo del paciente, y obviamente esto no se puede hacer en la realidad), como por arriba, para no inyectar demasiada insulina. Para el valor límite superior usamos el escogido también por F. Barragán et al. [2] y P. Dua et al. [5], que es 100 mU/min. El rango de valores de 0 a 100 mU/min se considera adecuado para la administración de insulina usando bombas, ya que hay que tener en cuenta también el diámetro del catéter por el que se administra, el cual es de un tamaño reducido para ser lo menos invasivo posible. La referencia de glucosa será 90 mg/dL, y se mantendrá para los siguientes experimentos. En la figura 3-4 se puede observar este modelo de control PID continuo en Simulink. Figura 3-3. Bode de Gm(s)*C(s) final 19 19 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) 3.1.4 Pruebas con el controlador PID en tiempo continuo Antes de mostrar las gráficas obtenidas de la simulación del control PID en tiempo continuo, enseñaré primero las gráficas correspondientes a la simulación del modelo de la figura 2-2, donde no hay ningún tipo de control aplicado, solo el propio comportamiento del sistema en sí. Se realizan dos pruebas: sin perturbaciones y con las perturbaciones de comida y ejercicio. Esta primera prueba servirá para comparar con las posteriores en donde se aplica el control y comprobar que realmente este control mejora la respuesta del sistema según las especificaciones. Aunque el tiempo de simulación en Simulink sea en segundos, en todas las gráficas de pruebas, tanto en este capítulo como en los siguientes, el tiempo en el eje X se muestra en minutos, que resulta una escala mucho más adecuada a la hora de observar los resultados. Figura 3-4. Modelo con controlador PID continuo en Simulink Figura 3-5. Comportamiento del modelo de Bergman sin perturbaciones Control PID 20 20 En la figura 3-5 se observa el comportamiento del modelo de Bergman sin aplicar perturbaciones, es decir, con D(t) = 0 y Ej = 0. La señal de control U(t) se aplica constante y es igual al valor de equilibrio, es decir, 0.277 mU/s. Si bien las simulaciones son pocos útiles ya que no son realistas (en la vida real el paciente tendrá distintos tipos de perturbaciones a lo largo del día, y es muy raro que no tenga alguna en algún momento), se realizan para comprobar posteriormente con el control que la salida alcanza rápidamente el valor de referencia y se mantiene constante alrededor de él, tal y como se espera. En este caso, el sistema definido por las ecuaciones del modelo de Bergman parte de un valor de glucosa cercano a 120 mg/dL y pasados unos 200 minutos alcanza un valor permanente de 80 md/dL, que es un valor aceptable de glucosa. Esta curva estaría muy bien y no necesitaría control si en la vida real no hubiera perturbaciones, pero como he dicho antes, en la realidad el paciente tendrá constantemente perturbaciones de comida, ejercicio y demás. Por su parte en la figura 3-6 vemos el comportamiento del modelo ahora sí con las perturbaciones, por lo tanto es mucho más realista. El tiempo de la perturbaciones, como en todas las gráficas siguientes, viene dado por el protocolo de pruebas definido en el capítulo 2.2.2. Se observa como al principio, cuando todavía no hay perturbación alguna, la glucosa tiene a irse a 80 mg/dL como en la gráfica anterior. En cuanto empieza a haber perturbaciones de comida, la glucosa tiende a subir mucho, alcanzando valores máximos de casi 250 mg/dL, que están muy por encima de los valores deseados. Por ejemplo, en el almuerzo, la glucosa no vuelve al valor deseado (entorno a 100 mg/dL) hasta pasados un poco más de 200 minutos, es decir, más de 3 horas. El paciente estaría todo ese periodo de tiempo a niveles de glucosa por encima de lo deseado, lo cual puede tener consecuencias negativas en su salud si esto se repite día tras día. Esta gráfica correspondería a una bomba de insulina trabajando en bucle abierto, con un valor constante de insulina administrándose durante todo el tiempo. Es más que evidente que no resulta muy efectiva, ya que si bien a lo largo de las horas hace que la glucosa acabe bajando, ese tiempo es demasiado grande y además los valores máximos de glucosa que se alcanzan son demasiado alejados de los valores deseados. A continuación se realiza la simulación correspondiente al modelo de la figura 3-4, el del control PID continuo. En las figuras se mostrarán dos gráficas fundamentales en las pruebas. Por una parte, arriba está la correspondiente a la glucosa del paciente, donde se muestra la glucosa real (la salida G(t) del sistema), la glucosa de referencia, y los valores límite aceptables en una persona con diabetes tipo 1, que serán 60 mg/dL por abajo (en amarillo, marcando la hipoglucemia) y 180 mg/dL por arriba (en morado, marcando la hiperglucemia), que ayudarán a ver con facilidad si la salida del sistema está dentro de ese rango. Por otra parte, abajo se encuentra la gráfica de la señal de control, es decir, la insulina que se administra y que ha sido calculada previamente por Figura 3-6. Comportamiento del modelo de Bergman con perturbaciones 21 21 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) el algoritmo de control. Para la obtención de las mismas se usa el script “graficas.m”. En la figura 3-7 se puede observar cómo es la respuesta del sistema sin perturbaciones, y se comprueba que ésta es capaz de seguir la referencia de 90 mg/dL sin ningún problema, y para ello se administra un valor constante de insulina de alrededor 10 mU/min una vez alcanzado el valor de referencia de glucosa para mantenerla. En la figura 3-8 se ve cómo el sistema responde cuando sí hay perturbaciones, y podemos observar a simple vista una mejora importante con respecto a la figura 3-6. Para empezar, cuando se aplica el control, la glucosa del paciente mantiene unos valores aceptables durante toda la simulación, siendo el valor máximo que se alcanza 160 mg/dL después del almuerzo (que es la comida en donde más hidratos de carbono ingiere). Esto sin control no ocurría, y se alcanzaban casi los 250 mg/dL, es decir, se ha conseguido bajar hasta 100 mg/dL la glucosa del Figura 3-7. Control PID continuo sin perturbaciones Figura 3-8. Control PID continuo con perturbaciones Control PID 28 28 3.3.3 Control PID con inyección de insulina por lotes Hasta ahora se había hablado de un controlador PID en tiempo discreto puesto que la señal de control no se calculaba de forma continua, es decir, con un tiempo de muestreo igual a 0, sino que esta se calculaba cada cierto tiempo cuando el sensor tomaba una muestra nueva y se actualizaba el error en el algoritmo de control implementado en la computadora. Sin embargo, aunque el valor de la insulina administrada se calcule cada 5 minutos como sucede en nuestro caso, sí que estamos suponiendo hasta este momento (tal y como se puede apreciar en todas las figuras correspondientes a las pruebas realizadas) que la forma de administrar esta insulina por parte de la bomba es continua. Es decir, que una vez se calcula el valor a inyectar, la bomba está administrando esa cantidad de insulina de forma continua durante los 5 minutos siguientes. Sin embargo, en la realidad esto no es lo que sucede. Lo normal es que la insulina se administre de forma “discreta”, esto es, por lotes. El procedimiento a seguir será el siguiente: cuando se calcule el nuevo valor de insulina, cada tiempo de muestreo, se calculará a su vez el valor total equivalente que se administraría de forma continua durante los 5 minutos. Esta cantidad total será la que se administre de una sola vez al comienzo del periodo de muestreo, y tardará en administrarse el tiempo que venga impuesto por el caudal de la bomba (la máxima cantidad de insulina que sea capaz de administrar en un instante de tiempo). De esta forma, la insulina no se inyectará de forma continua, sino en pequeños lotes cada 5 minutos, haciendo que el control sea más realista. Para esta mejora del control PID se crean dos nuevos archivos nuevos en Matlab, uno para el control llamado “controllotes.m”, con estructura muy similar al usado anteriormente, y otro para mandar la señal de control según el tiempo de inyección calculado, llamado “inyectalotes.m”. En el primer archivo se calcula la señal de control tal y como se hacía en el archivo “control.m”, con un tiempo de muestreo de 5 minutos y anti-windup. La novedad viene en las últimas líneas, donde se calcula el tiempo durante el cual se debe inyectar insulina al caudal máximo para inyectar toda la cantidad total necesaria durante ese periodo de muestreo. El caudal máximo será Q=100 mU/min, el mismo definido anteriormente como valor máximo de la señal de control a la hora de saturar. Esta función, que se llama desde Simulink con un bloque de función igual que en los casos anteriores, devuelve el tiempo de inyección calculado que servirá de entrada, aparte del tiempo de simulación, a la segunda función. En “inyectalotes.m” primero se calcula el tiempo que ha pasado desde que comenzó el último periodo de muestreo, que irá entre 0 y 5 minutos, y que calcula como el resto del tiempo de simulación entre el tiempo de muestreo. A continuación, se tiene la condición siguiente: si el tiempo que se lleva de periodo es menor que el Figura 3-16. PID con tiempo de muestreo de 10 minutos 29 29 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) tiempo calculado de inyección, entonces el valor de Uk (la señal de control y salida de la función) será 100 mU/min. En caso de que el tiempo sea mayor, entonces querrá decir que ya se ha administrado toda la insulina necesaria en ese periodo y no hay que inyectar más hasta el siguiente. Hay que decir que bloque de esta función en Simulink sí que tiene un tiempo de muestreo igual a 0, es decir, se ejecuta de forma continua, ya que es el propio algoritmo quien se encarga que la administración sea “discreta”, como ya se ha explicado. El modelo de Simulink se observa en la siguiente figura: Se realizan pruebas y se obtienen los resultados que se muestran a continuación. Figura 3-17. Bloques del control PID con inyección de insulina por lotes Figura 3-18. PID con inyección de insulina por lotes Control PID 30 30 Tal y como podemos ver en la figura 3-18, la respuesta de la glucosa en sangre en el paciente con esta nueva versión más realista del control es prácticamente idéntica a la obtenida en el caso anterior, cuando la insulina se administraba de forma continua. Esto es excelente, porque podemos implementar el controlador PID con inyección de insulina por lotes con total seguridad ya que sabemos que no va a afectar a la glucosa del paciente, con la mejora de que se diseñando un control más realista. Puesto que en la figura no se puede apreciar del todo bien como es la señal de control y como varía esta cuando se produce una perturbación, se muestra en la figura 3-19 un zoom de la misma realizado en el intervalo de tiempo correspondiente a la primera perturbación de Figura 3-19. Zoom del PID con inyección de insulina por lotes Figura 3-20. Tiempo de inyección de insulina 31 31 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) comida (el desayuno). Se observa como justo antes de la perturbación, cuando el valor de glucosa es muy cercano a la referencia, la cantidad de insulina que se debe inyectar es baja y por lo tanto, el tiempo de inyección, y la señal de control tiene forma de pulso prácticamente. En cambio, cuando se produce la perturbación, se puede ver que se necesita una mayor cantidad de insulina y de esta manera el tiempo de inyección aumenta, teniendo ahora la señal de control forma de escalón. Conforme se va corrigiendo la perturbación, más estrecho es el escalón de la insulina administrada puesto que menos insulina se va necesitando. En cualquier caso, en la figura 3-20 se muestra el tiempo de inyección calculado en la primera función, que tal y como se esperaba nunca sobrepasa el valor del tiempo de muestreo, y que aumenta cuando se producen las perturbaciones y viceversa. Por lo tanto, y tras estudiar el resultado de las pruebas, se decide que la implementación de este control PID con inyección de insulina por lotes es totalmente válido para nuestro caso. 3.3.4 Controlador PID Feedforward El control feed forward, también llamado control anticipativo de prealimentación, es un tipo de control que mide las perturbaciones que le entran al sistema a controlar de forma que puede anticiparse y compensarlas, mejorando la respuesta del sistema ante el problema de regulación. Este control PID feedforward sumaría la señal de control propia del PID que ya se ha obtenido en los apartados anteriores, con una señal de control calculada a partir de la perturbación medida (en nuestro caso dos señales de control, puesto que hay dos perturbaciones), obteniéndose así una señal de control total que se aplica al sistema. Hay que destacar un punto importante en nuestro sistema en particular, y es el hecho de medir las perturbaciones. Esto se podría traducir a la realidad como que el paciente, justo antes de comer o de hacer ejercicio (es decir, cuando se producen las perturbaciones), debería introducir en la bomba o el dispositivo que le administre la insulina la cantidad más o menos exacta de los hidratos de carbono que va a ingerir o el tiempo que va a hacer ejercicio, respectivamente. De esta forma, el algoritmo de control tendría la información necesaria para medir las perturbaciones que va a tener el sistema, ya que no se pueden medir directamente de otra forma. En la simulación que haré, sin embargo, y puesto que en Simulink sí es posible medir directamente las perturbaciones, lo que le llegue al bloque de control, además de los parámetros que ya conocemos, serán las salidas directas de cada uno de los bloques de perturbaciones, como podemos ver en la figura 3-21. Para obtener las respectivas señales de control de cada perturbación se calcula una función de transferencia de forma que la entrada sea la perturbación medida y la salida sea la propia señal de control (de forma similar a como se hacía en el caso del controlador en tiempo continuo, aunque ahí la entrada era el error de glucosa). Para calcular Cffp1 y Cffp2, las funciones de transferencia para el control feedforward para la perturbación 1 (comida) y 2 (ejercicio físico), partimos de la siguiente ecuación observando el siguiente esquema: Figura 3-21. Modelo del controlador PID Feedforward Control PID 32 32 0=𝑦=𝑢∗𝐺𝑚+𝑝∗𝐺𝑝 (31) Siendo u la señal de control aplicada al sistema, Gm la función de transferencia del modelo, p la perturbación medida, Gp la función de transferencia de la perturbación, e y la salida. De esta forma en la ecuación se muestra el caso que queremos, es decir, que el incremento de la salida sea 0 y esté en equilibrio, y para esto la señal de control que le entra al sistema debe compensar la perturbación, tal que: 𝑢=−𝑝∗𝐺𝑝 𝐺𝑚 (32) Por tanto la función de transferencia, a la que en este caso general llamaré Cff, será igual a: 𝐶𝑓𝑓(𝑠)=−𝐺𝑝 𝐺𝑚 (33) En nuestro caso particular, se tienen las siguientes funciones de transferencia para las perturbaciones, calculadas en el script “funcion_transferencia.m” ya mencionado en el capítulo 2: 𝐺𝑑(𝑠)=15000 (15000𝑠+7) (34) 𝐺𝑒𝑗(𝑠)=−67500 (15000𝑠+7) (35) En el script “feed_forward.m” se calculan las dos funciones de transferencia para el control feed forward tal como se ha visto en la ecuación (33), y se obtiene: 𝐶𝑓𝑓𝑑(𝑠)=3.49∗1014𝑠3+8.49∗1011𝑠2+5.45∗108𝑠+105000 (473850𝑠+221.1) (36) 𝐶𝑓𝑓𝑒𝑗(𝑠)=−(1.575∗1015𝑠3+3.821∗1012𝑠2+2.453∗109𝑠+472500) (473850𝑠+221.1) (37) Una vez se tienen las funciones se puede implementar el controlador en Simulink, el cual podemos ver en la siguiente figura 3-23. Hay que decir que dentro de los bloques de función de transferencia para Cffd y Cffej , puesto que estas no son bipropias (mismo número de polos que de ceros) y Simulink pide que lo sean, se añaden polos de alta frecuencia hasta que el denominador tiene tantos polos como ceros el numerador. Los polos de alta frecuencia que se añaden son, en este caso, los mismos que se añadieron en el controlador PID de tiempo continuo y que se pueden ver en la ecuación (23), y no afectan al comportamiento del sistema. De esta forma, y como ya se ha explicado antes, se tiene la señal de control propia del PID, que se obtiene en el bloque “ControlPID”, a la que se suman las señales de control de las perturbaciones calculadas. Puesto que en el apartado anterior se había estudiado que la inyección de insulina por lotes era válida e incluso favorable al hacer el control más realista, se opta por seguir manteniéndola en este apartado también (así como en el siguiente capítulo). Para esto, la señal de control total se pasa por el bloque “Control”, donde dentro se llama a la función “lotes.m” que primero satura la señal y después calcula el tiempo de inyección. Por último, se tiene de nuevo el bloque “Lotes” que a partir del tiempo de inyección calcula el valor de insulina a administrar. Figura 3-22. Esquema para la obtención de Cff 33 33 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) Cuando se simula se obtienen los siguientes resultados: A primera vista, y comparando con la figura 3-18, observamos una clara diferencia en cuanto al tiempo en corregir las perturbaciones, especialmente la comida. En este caso el controlador es capaz de corregir las perturbaciones en menos tiempo de forma notable, reduciendo el tiempo por lo tanto en el que la glucemia se encuentra a niveles por encima del de referencia. En el caso del desayuno, por ejemplo, observamos que con el feed forward el pico se corrige y la glucosa en sangre vuelve al valor de referencia hacia el minuto 600, mientras que sin el feed forward este no se corregía hasta el minuto 700. Esto 100 minutos menos (casi 2 horas) son vitales y de gran importancia para la salud del paciente, ya que significan maximizar el tiempo con una glucosa a niveles de una persona no diabética, reduciendo las hiperglucemias y los problemas que ellas conllevan. En el caso del almuerzo y cena también se reduce de forma igualmente notable este tiempo. También podemos destacar que los picos (valores máximos y mínimos alcanzados) son similares a los obtenidos con el controlador sin feed forward, o en cualquier caso ligeramente menos pronunciados, lo cual también es bueno para nuestro sistema. Por estos motivos, y sin olvidar que el feedforward conllevaría que el paciente proporcione a la bomba los Figura 3-23. Bloque del control PID Feedforward Figura 3-24. Controlador PID Feedforward con inyección de insulina por lotes Control PID 34 34 valores que tomarán las perturbaciones como he mencionado anteriormente (sin que esto conlleve un problema en sí, ya que es factible), se determina que el controlador feedforward con inyección de insulina por lotes es una excelente opción para este sistema en concreto, ya que aporta unas mejoras muy importantes a la hora de mantener la glucemia del paciente en los valores correctos. 35 35 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) 4 CONTROL PREDICTIVO BASADO EN EL MODELO (MPC) l Control Predictivo basado en el Modelo o MPC (del nombre inglés, Model based Predictive Control) es un tipo de control avanzado ampliamente usado, especialmente en procesos industriales, pero también en ámbitos científicos, que se basa en el modelo del sistema a controlar para predecir las salidas o respuestas de este para de esta forma optimizar la señal de control a aplicar en el futuro. En este capítulo se aplicará este control al sistema formado por el Modelo de Bergman usando para ello la herramienta que ofrece Matlab y Simulink llamada “MPC Toolbox”. 4.1 Introducción al Control Predictivo Si bien en el anterior párrafo he dicho que el Control Predictivo basado en el Modelo (a partir de ahora se llamará por sus siglas MPC para resultar más corto) es un tipo de control, realmente esto no es así ya que el MPC se puede definir como un conjunto de distintos tipos de control, como el óptimo o el de procesos multivariables. Tal y como se ha mencionado antes, el MPC se basa en dos principios básicos: el uso del modelo para poder predecir las salidas en un futuro del sistema, y la optimización de la señal de control. S. Diego Fernando [6] apunta en su artículo que el objetivo principal actualmente en el mundo del control es tener la posibilidad de actuar sobre las variables manipuladas o señales de control, para que estas puedan adaptarse según los requisitos del proceso, que cada vez son más exigentes y complejos (ya sea por calidad del producto, causas económicas, ambientales…). El MPC es una gran opción a la hora de enfrentarse a este objetivo, puesto que aporta flexibilidad en cuanto al modelo y a las restricciones, haciendo que sea probablemente la estrategia de control basada en modelo más escogida a la hora de controlar este tipo de procesos [7]. Algunas de las ventajas que se pueden destacar del MPC son las siguientes: es multivariable, funciona tanto con modelos lineales como no lineales, al ser predictivo sirve para solucionar el problema de los tiempos muertos, ayuda a compensar perturbaciones medibles, se le pueden añadir restricciones rígidas… Como contrapunto, el MPC conlleva una implementación más compleja, sobre todo comparado con los controladores tipo PID, así como un mayor coste E “Si quelqu’un aime une fleur qui n’existe qu’à un exemplaire dans les millions et les millions d’étoiles, ça suffit pour qu’il soit heureux quand il les regarde”. - “El Principito”, Antoine de Saint-Exupéry - Control Predictivo basado en el Modelo (MPC) 36 36 computacional. Aunque como todo, será el proceso en sí a controlar y los requisitos que este conlleve los que harán que se termine optando por un tipo de control u otro [6]. En el control predictivo tenemos variables manipulables (equivalente a las señales de control) y variables controladas (las salidas del sistema), así como las perturbaciones. Puesto que es un control en tiempo discreto, se tiene también un tiempo de muestreo, que será cada cuanto tiempo se realice el algoritmo de control. Se define un horizonte de predicción, que se corresponde con el rango de tiempo (o el número de instantes en discreto) que el controlador es capaz de predecir la salida del sistema desde el instante actual. Es decir, si nos encontramos en el instante 200 segundos y el horizonte de predicción es de 300 segundos, esto quiere decir que el controlador MPC puede predecir la salida del sistema hasta el instante 500 segundos. Se debe encontrar un valor adecuado para este parámetro, pues si es demasiado bajo, el controlador no podrá controlar de forma adecuada si se encuentra una perturbación o un cambio de referencia para los cuales el sistema necesite más tiempo de adaptación o corrección que el propio horizonte de control; igualmente, si es muy alto, se pueden estar haciendo predicciones en vano si se produce una perturbación de repente y ya no sirven las predicciones hechas. En cada instante de muestreo, el controlador primero mide la salida real del sistema y la compara con el valor que había predicho en el instante anterior que ésta tendría (podrá o no ser igual, dependiendo de las perturbaciones que haya y que no se puedan medir). Una vez hecho esto, predice una serie de escenarios del comportamiento del sistema aplicando en cada uno de ellos distintas señales de control. Se busca que el error entre la referencia y la salida real del sistema sea el mínimo posible, pero también que se cumplan, si las hay, una serie de restricciones. De esto se encarga el optimizador, ya que mediante unas funciones de coste en donde se asignan unos pesos dependiendo de si queremos penalizar acciones en la señal de control o en la salida (por ejemplo, si no queremos que la señal de control realice cambios muy bruscos, o si preferimos que sí haya cambios bruscos a cambio de darle un mayor peso a que la salida se mantenga al valor de referencia el máximo tiempo posible) es capaz de elegir la secuencia de control más adecuada en cada caso, correspondiéndose esta a la función con un coste menor. El tiempo o el número de instantes de muestreo que dura la secuencia de control se llama horizonte de control, e igual que el horizonte de predicción se debe definir al principio. Cuanto mayor sea el horizonte de control, mayor será el coste computacional, pero tampoco se puede definir demasiado bajo ya que entonces el control podría ser menos efectivo. Una vez se tiene la secuencia a aplicar elegida por el optimizador, se aplica al sistema solo el primer valor de la señal de control calculada (no la secuencia entera). En el siguiente instante de muestreo se volverán a repetir todos los pasos, empezando de nuevo con la medición de la salida real. El horizonte de predicción se va deslizando, moviéndose cada vez a un instante de muestreo más [7]. 4.2 Implementación del controlador MPC en Matlab y Simulink Una vez hecha una breve introducción de este tipo de control, pasaré a explicar como se ha implementado usando para ello el Toolbox que ofrece Matlab y Simulink. Siguiendo los modelos anteriores de Simulink para los otros controladores (y más en concreto, para el PID feedforward), ahora se cambia el interior del bloque de control. A este bloque le sigue llegando el valor de glucosa real, la glucosa de referencia y las perturbaciones, estas últimas como una señal multiplexada, como se puede ver en la figura 4-1, y saca la insulina que se administra siguiendo con la inyección por lotes. Por su parte, dentro el bloque de control, el bloque de la función de Matlab que llamaba a la función que implementaba el PID se cambia por un bloque de controlador MPC, perteneciente al Toolbox mencionado. A este bloque le llega como entradas el valor de salida del sistema, la referencia y las perturbaciones, y saca la señal de control que posteriormente se pasará por las funciones correspondientes a la inyección por lotes al igual que con el PID feedforward, como vemos en la figura 4-2. 37 37 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) El bloque MPC permite diseñar de forma intuitiva un controlador predictivo. Para ello, primero se introducen el número de variable manipulables, variables controladas y perturbaciones. En este caso, tal y como se muestra en la figura 4-3, se tienen una variable manipulable (la insulina), una variable controlada medida (la glucosa) y dos perturbaciones medidas (la comida y el ejercicio). Se escogen las señales del modelo que corresponden a cada una de las señales mencionadas anteriormente, y Simulink linealiza automáticamente el sistema. Una vez hecho esto, se abre una ventana donde se pueden realizar distintas acciones para modificar y diseñar el controlador, y dentro de la pestaña “Tunning” se puede ir modificando los horizontes de control y predicción, así como los pesos de la función de coste del optimizador y las restricciones en las variables. Esta barra de herramientas se puede observar en la figura 4-4. Se pone como restricción que la variable manipulable solo pueda valer de 0 a 100 mU/min (tal y como estaba en el saturador del PID), y se define el tiempo de muestreo como 300 segundos. El horizonte de predicción se define como 300 segundos también al principio (para que prediga al menos hasta el siguiente tiempo de muestreo), mientras que el de control se deja el que viene por defecto, que es 10 segundos. Los pesos igualmente se dejan los que vienen por defecto en esta primera simulación. Se realizan distintas simulaciones, cambiando para cada una de ellas el valor de los parámetros anteriormente mencionados, para ir observando cómo afectan estos cambios y cuáles son los valores más adecuados. Se obtienen los siguientes resultados: Figura 4-1. Modelo del controlador MPC Figura 4-2. Bloque de control del MPC Comparaciones y conclusiones finales 44 44 que es la mayor. Tabla 5-1. Comparación cuantitativa de los controladores Como se puede ver, en la tabla también se añade los valores correspondiente a la respuesta del modelo en bucle abierto, sin ningún tipo de control, para poder observar el cambio favorable y la mejora en la respuesta cuando se implementa un control en bucle cerrado en este sistema, que también era uno de los objetivos de este trabajo. Por una parte, en los 3 controladores se reduce el valor de pico máximo, casi 100 mg/dL respecto al bucle abierto, alcanzándose en los 3 un valor muy similar y que resulta adecuado pues no se llega ni si quiera a una hiperglucemia (que sería por encima de 180 mg/dL). Lo mismo se puede decir de los valores de pico mínimo, alcanzados cuando se realiza el ejercicio, también prácticamente iguales en los 3 controladores y con un valor cercano a la hipoglucemia, pero que tampoco la alcanza (ésta se marca a los 60 mg/dL). En cuanto al tiempo de establecimiento antes de las perturbaciones, aquí el MPC es el que destaca con claridad, tardando unos 100 minutos menos que los PID. Sin embargo, esto no quiere decir que estos tengan una respuesta muy alejada de la referencia, ya que como podemos ver en la figura 3-24 el valor de la glucosa durante este tiempo es cercano a la referencia, unos 5 o 10 mg/dL por debajo, siendo esto poco significativo ya que esos valores de glucosa son igualmente muy aceptables. En el MPC sí que llega al valor de referencia y se mantiene ahí hasta la perturbación. Por último para finalizar esta comparación cuantitativa, se estudia el tiempo de corrección de la segunda perturbación, un parámetro más importante si cabe que el anterior, ya que cuanto menor sea este tiempo, mejor será a la larga la salud del paciente. Es aquí donde se nota también que el control tiene información de las perturbaciones para poderse anticipar y compensarlas con una señal de control más alta. Mientras que en el PID sin feedforward se tarda 210 minutos (3 horas y media), en el PID feedforward y el MPC se tarda una hora menos, lo cual es bastante significativo y resulta muy beneficioso para la salud del paciente como acabo de decir. Desde luego, esta mejora hace que estos dos controladores sean mucho más efectivos que el PID sin feedforward, sin olvidarse que para ello el paciente debe proporcionar cuánto valen las perturbaciones. Una vez realizada la comparación, toca sopesar distintos aspectos y elegir uno de los controladores. Antes de nada hay que destacar un aspecto importante, y es el hecho de poder mejorar los controladores diseñados. Como la mayoría de trabajos de ingeniería, este trabajo queda abierto a futuras mejoras con técnicas de control más complejas o con posibles innovaciones que surjan en este campo del control de la diabetes. De hecho, los controladores que se han diseñado a lo largo del trabajo podrían realizar un mejor control del sistema si no fuese por una restricción que se puso desde el principio y que se ha ido comentando a lo largo del mismo: la limitación de la señal de control. Desde el diseño del primer controlador se impuso que la insulina máxima que se podía administrar era 100 mU/min. Un valor más alto en esta limitación (o incluso no tener limitación) supondría un control mucho más exacto, sobre todo en el PID feedforward o el MPC, que reducirían al mínimo los picos que se provocan por las perturbaciones consiguiendo que la glucosa se mantuviera prácticamente todo el tiempo al nivel de referencia. Pero ahora bien, esto implicaría una bomba con un caudal más grande, y por lo tanto un tubo Comparación cuantitativa de los controladores Tipo de control Valor máximo de pico Valor mínimo de pico Tiempo establecimiento antes de las perturbaciones Tiempo de corrección después del almuerzo Sin control (bucle abierto) 241.9 mg/dL 72.32 mg/dL - 274 min PID sin feedforward 166.1 mg/dL 69.14 mg/dL 346.5 min 210 min PID feedforward 166.3 mg/dL 70.74 mg/dL 333 min 163 min MPC 168.4 mg/dL 69.45 mg/dL 230 min 146 min 45 45 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) o catéter más grande e invasivo para el paciente, lo cual no sería para nada cómodo y, quizás, tampoco sería viable si tuviera que ser demasiado ancho por ejemplo. Puesto que lo principal en este trabajo es la salud y el bienestar del paciente con diabetes tipo 1 (para ello se realiza un control automático de la glucosa, para que su salud pueda ser mejor), y puesto que igual de importante es el buen control como que el paciente se sienta bien y pueda realizar una vida normal sin tener que estar cargando con un dispositivo incómodo, invasivo o que le limite en sus actividades diarias, esta limitación en la señal de control es obligada, sabiendo, eso sí, que el control también tendrá un punto a partir del cual ya no se podrá mejorar debido a esto. Por la comparativa se puede deducir que las mejores opciones son el PID feedforward y el MPC, puesto que son los dos que corrigen más rápido las perturbaciones. Entre estos dos controladores, sabiendo que el MPC supone en la realidad una implementación más compleja y que conlleva un mayor coste computacional, yo optaría por quedarme con el PID feedforward para nuestro sistema. Este ofrece unos resultados muy similares a los del MPC y supondría un ahorro económico y más facilidad de implementación. Esto último es otro aspecto importante a destacar, puesto que si se quisiera llevar a la realidad este control automático e implementarlo en un dispositivo de inyección de insulina (ya sea bomba o dispositivos que se desarrollen a lo largo de estos años siguientes), es necesario que este conlleve un coste no muy alto para que se facilite esta implementación, y por consiguiente que el precio del producto no fuera demasiado alto y así poder llegar al máximo de pacientes con diabetes tipo 1 posible y que estos se pudieran beneficiar de este dispositivo, suponiendo una gran ayuda en su día a día. Puesto que esto último resulta mucho más importante que los 5 o 10 mg/dL de glucosa que hay de diferencia respecto a la referencia al comienzo de la simulación en el PID feedforward (y que hemos visto que no hay en el MPC), se añade otra razón de peso más para elegir este controlador. Como conclusión, después de los estudios y comparaciones realizados a lo largo del trabajo, se termina optando como mejor opción de los controladores diseñados para el control automático de la glucosa en pacientes con diabetes tipo 1 el PID feedforward, tanto por sus buenos resultados a la hora de controlar, haciendo que la glucosa se mantenga en el rango adecuado y corrigiendo las perturbaciones de una forma eficaz, como por su facilidad de implementación y menor coste, posibilitando que así un mayor número de personas pudieran tener acceso a él si se llevara a cabo en la realidad. Para que este controlador funcionase de la forma más efectiva posible (la que hemos visto en este trabajo) haría falta que el paciente proporcionara la cantidad aproximada de hidratos de carbono que va a ingerir antes de comer y el ejercicio físico, sin que esto supusiera un esfuerzo extraordinario pues sería algo rápido y sencillo de llevar a cabo; si bien tampoco ocurría nada grave si alguna vez por algún motivo no proporcionara esta información, ya que simplemente el controlador actuaría en este caso como el PID sin feedforward (que hemos visto que realiza un buen control también, solo tardaría algo más en corregir la perturbación correspondiente). Comparaciones y conclusiones finales 46 46 47 47 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) REFERENCIAS [1] Organización Mundial de la Salud. Diabetes. (13 de abril de 2021). [En línea]. Disponible: https://www.who.int/es/news-room/fact-sheets/detail/diabetes [Consulta: 16 de abril de 2021]. [2] F. Barragán, M.G. Ortega Linares, y F. Rubio. "Técnicas de control robusto para la regulación de glucosa en sangre en pacientes con Diabetes Mellitus tipo 1". [3] A. Campos Carceller. "Diseño de un sistema de control multivariable del nivel de glucosa en sangre para pacientes con Diabetes Mellitus tipo 1". Trabajo de fin de grado en Ingeniería en Tecnologías Industriales. Escuela Superior de Ingenieros Industriales de Valencia, 2014. [4] A. Nath et al. "Physiological Models and Control for Type 1 Diabetes Mellitus: A Brief Review". IFAC PapersOnLine 51-1, 2018, pp. 289-294. [5] P. Dua, F. J. Doyle, E. N. Pistikopoulos. "Multi-objective Parametric Control of Blood Glucose Concentration for Type 1 Diabetes", en la 44ª Conferencia de IEEE sobre Decisión y Control, y la Conferencia Europea de Control, 2005. [6] S. Diego Fernando. "¿Qué es el Control Predictivo y hacia dónde se proyecta?". Revista Especializada en Ingeniería. Universidad Abierta y a Distancia. Publicaciones e Investigación ISSN: 1900-6608 Volumen 7, 2013. [7] J. M. Vargas Lara. "Control predictivo multivariable: evolución histórica y conceptos". TECSUP. Referencias 48 48 49 49 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) ANEXOS A continuación se adjuntan como anexos los códigos de Matlab que se han realizado: “Modelo_Bergan.m” function [salida] = Modelo_Bergman(in) % Variables de entrada en la funcion D = in(1); U = in(2); G = in(3); I = in(4); X = in(5); Ej = in(6); %Parametros del modelo P1=0.028/60; %s^-1 P2=0.025/60; %s^-1 Gb=4.5; %mmol/L n=(5/54)/60; %s^-1 Ib=15; %mU/L Vi=12; %L P3=0.000013/3600; %L/(mU s^2) % Ecuaciones diferenciales % Gd=(-P1*G) - X*(G+Gb) + D; Gd=(-P1*G) - (X+Ej)*(G+Gb) + D; Id=(-n*(I+Ib)) + (U/Vi); Xd=(-P2*X) + (P3*I); salida=[Gd;Id;Xd]; “comida.m” function [D] = comida(in) % Variables de entrada en la funcion t = in(1); %Parametros del modelo Ag=0.8; %adimensional Tmaxi=55*60; %s Vg=13.79; %L Tmaxg=40*60; %s t1=0; Dg=0; f=1000/180; %factor de conversion de g a mmol %Comidas (perturbaciones) a lo largo de un día Anexos 50 50 if t>=480*60 && t<840*60 %desayuno: 55 g de HC a las 8:00 t1=480*60; Dg=55*f; end if t>=840*60 && t<1260*60 %almuerzo: 90 g de HC a las 14:00 t1=840*60; Dg=90*f; end if t>=1260*60 %cena 45 g de HC a las 21:00 t1=1260*60; Dg=45*f; end % Ecuacion D= (Dg*Ag*(t-t1)*exp(-(t-t1)/Tmaxi))/(Vg*Tmaxg^2); “ejercicio.m” function [Ej] = ejercicio(in) % Variables de entrada en la funcion t = in(1); %Parametros del modelo cons = 0.005/60; %Ejercicio (perturbaciones) a lo largo de un día if t>=1080*60 && t<1140*60 %ejercicio de 18:00 a 19:00 Ej=cons+ 0.005/60; else Ej=cons; end “función_transferencia.m” %Modelo linealizado y funcion de transferencia format short %Parametros del modelo (en segundos) P1=0.028/60; P2=0.025/60; Gb=4.5; n=(5/54)/60; Ib=15; Vi=12; P3=0.000013/3600; syms s; %Matrices del modelo linealizado A=([-P1,0,-Gb;0,-n,0;0,P3,-P2]); B=[0,1,-Gb;1/Vi,0,0;0,0,0]; C=[1,0,0]; D=[0;0;0]; I=eye(3); %Funcion de transferencia entera G=C*(inv((s*I)-A))*B; G=vpa(G,5); 51 51 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) %Gm=tf(-2.6325,conv(conv([250 7],[40 1]),[54 5])) en rad/min Gm=tf(-31.59,conv(conv([15000 7],[2400 1]),[648 1])) %en rad/s Gd= tf(15000,[15000 7]) Gej=tf(-67500, [15000 7]) “control.m” function [U] = control(in) % Variable de entrada en la funcion: Gref Gref = in(1); G =in(2); tiempo = in(3); % Tiempo de muestreo (revisar el utilizado en Simulink) persistent Int_err; Tm=300; %Cambiar en simulink también persistent Err_G_ant; if(tiempo < 1e-2) %al principio el error integral es 0 Int_err= 0; Err_G_ant= 0; end %Implementacion PID Kp=-0.3385; Td=3014.2; Ti=12057; Err_G= Gref-G; Kd=Kp*Td; Ki=Kp/Ti; u=Kp.*Err_G + Ki.*Int_err + Kd.*(Err_G-Err_G_ant)/Tm; Err_G_ant=Err_G; U=u + 0.277; Int_err=Int_err + (Err_G*Tm); %Anti windup if (U<=0) || (U>=100/60) Int_err=Int_err - (Err_G*Tm); end %Saturacion if U<0 U=0; end if U>100/60 %mU/s U=100/60; end “controllotes.m” function [tinyec] = controllotes(in) % Variables de entrada en la funcion: Gref = in(1); G =in(2); tiempo = in(3); % Tiempo de muestreo (revisar el utilizado en Simulink) persistent Int_err; Tm=300; %Cambiar en simulink también Anexos 52 52 persistent Err_G_ant; if(tiempo < 1e-2) Int_err= 0; Err_G_ant= 0; end %Implementacion PID Kp=-0.3385; Td=3014.2; Ti=12057; Err_G= Gref-G; Kd=Kp*Td; Ki=Kp/Ti; u=Kp.*Err_G + Ki.*Int_err + Kd.*(Err_G-Err_G_ant)/Tm; Err_G_ant=Err_G; U=u + 0.277; Int_err=Int_err + (Err_G*Tm); %Anti windup if (U<=0) || (U>=100/60) Int_err=Int_err - (Err_G*Tm); end %Saturacion if U<0 U=0; end if U>100/60 %mU/s U=100/60; end % %Inyeccion de la insulina "por lotes" Uq=Tm*U; %cantidad de insulina administrada durante 5 min (tm) Q= 100/60; %el caudal maximo de inyeccion es de 100 mU/min tinyec=Uq/Q; %tiempo en segundos que debe estar inyectando insulina para inyectar toda la necesaria en ese minuto “inyectalotes.m” %funcion a la que le entra el tiempo real y el tiempo de control calculado %durante el que se tiene que inyectar insulina function [Uk] = inyectalotes(in) tinyec = in(1); tiempo =in(2); Tm=300; t=mod(tiempo,Tm); %tiempo que lleva en el periodo de muestreo if t<=tinyec Uk=100/60; else Uk=0; end “feed_forward.m” %Calculo de Cffd y Cffej para el control feed forward Gm=tf(-31.59,conv(conv([15000 7],[2400 1]),[648 1])); %en rad/s Gd= tf(15000,[15000 7]); Gej=tf(-67500, [15000 7]); 53 53 Modelo y control automático de la glucosa en sangre en pacientes con Diabetes Mellitus tipo 1 (DM1) Cffd=-Gd/Gm; Cffej=-Gej/Gm; Cffd=tf([3.49e14 8.49e11 5.45e8 105000],[473850 221.1]) Cffej=tf([-1.575e15 -3.821e12 -2.453e9 -472500],[473850 221.1]) “lotes.m” function [tinyec] = lotes(in) % Variables de entrada en la funcion: U = in(1); % Tiempo de muestreo (revisar el utilizado en Simulink) Tm=300; %Cambiar en simulink también %Saturacion if U<0 U=0; end if U>100/60 %mU/s U=100/60; end % %Inyeccion de la insulina "por lotes" Uq=Tm*U; %cantidad de insulina administrada durante 5 min (tm) Q= 100/60; %el caudal maximo de inyeccion es de 100 mU/min tinyec=Uq/Q; %tiempo en segundos que debe estar inyectando insulina para inyectar toda la necesaria en ese minuto “graficas.m” %Para dibujar las gráficas obtenidas t=out.t; G=out.G; I=out.I; X=out.X; U=out.U; Ref=out.Ref*18; tmin=t/60; Umin=U*60; Gmin=60; Gmax=180; % plot(tmin,(out.tiempoiny/60));grid;xlabel('t(min)');ylabel('tiempo inyección(min)');title('Tiempo de inyección de insulina');axis([0 1440 0 10]); subplot(2,1,1);plot(tmin,G,tmin,Ref,[0,1440],[Gmin,Gmin],[0,1440],[Gmax,Gmax]);grid;xl abel('t(min)');ylabel('Glucosa(mg/dL)');axis([0 1440 40 200]);legend('Glucemia real', 'Glucemia referencia');title('Glucemia del paciente'); subplot(2,1,2);plot(tmin,Umin);grid;xlabel('t(min)');ylabel('Insulina inyectada(mU/min)');axis([0 1440 0 120]);title('Insulina administrada');