Tesis Modelado e Identificacin de Bioprocesos
Tesis Modelado e Identificacin de Bioprocesos
Tesis Modelado e Identificacin de Bioprocesos
MODELADO E IDENTIFICACIÓN DE
BIOPROCESOS
Vigo, 2006
Autorización
CERTIFICA:
Son muchas las personas e instituciones que durante estos años han participado
en este trabajo y a todas ellas quiero expresar mi gratitud por el apoyo prestado.
En primer lugar deseo expresar mi más profundo agradecimiento al Dr. Julio
Rodrı́guez Banga por su labor de dirección en este trabajo y por su constante apoyo
y asesoramiento en todos los aspectos de la investigación y elaboración de esta Tesis.
A él debo agradecer la confianza depositada en mi al brindarme la oportunidad de
formar parte de su grupo ası́ como todo su tiempo y permanente disponibilidad.
Quisiera agradecer a todo el personal del Instituto de Investigaciones Marinas del
CSIC, en especial al grupo de Ingenierı́a de Procesos por todos los momentos com-
partidos, ¡os echaré de menos!. Me gustarı́a destacar aquı́ a Eva por la colaboración
prestada en todo momento.
Asimismo, agradezco al Ministerio de Educación y Ciencia por la financiación de
los proyectos AGL2001-2610-C02-02 y AGL2004-05206-C02-01/ALI y a la Xunta de
Galicia por el proyecto PGIDIT02PXIC40211PN y por concederme una beca que me
permitió realizar una estancia en la UCSB (University of California Santa Barbara).
Quiero mandar desde aquı́ mi más sincero agradecimiento al Prof. Francis J.
Doyle III, a su grupo de investigación de la UCSB y al equipo del Sansum Diabetes
Research Institute por la amabilidad y el afecto con el que me acogieron y por lo
mucho que he aprendido de ellos.
A la Comisión Europea, a través del CTS (Control Training Site), por darme la
oportunidad de realizar una estancia en Supèlec y a la Prof. Françoise Lamnabhi-
Lagarrigue y al Prof. Eric Walter por la supervisión del trabajo allı́ realizado.
A mi familia y amigos debo el cariño y los ánimos necesarios para que este trabajo
llegase a buen puerto, especialmente a mis padres y a mis hermanos, gracias es lo
menos y a la vez lo más que puedo deciros.
No puedo acabar los agradecimientos sin recordar a la persona que más ha sufrido
este trabajo. Gracias Juan, por estar siempre ahı́.
A todos los mencionados y a todos los que quedan en el tintero, muchas gracias
por el apoyo prestado, el cariño recibido y los ánimos proporcionados.
A mis padres y a mis hermanos
A Juan
All models are wrong... but some are useful.
George E. P. Box
Resumen
xi
xii Resumen
Debido a la naturaleza no lineal del modelo dinámico, estos dos problemas son fre-
cuentemente multimodales (no convexos) y, por lo tanto, si se resuelven con métodos
tradicionales de optimización local es muy probable que converjan a óptimos locales.
Además, en el caso de un mal ajuste de los parámetros, no hay modo de saber si
éste se debe a una mala formulación del modelo o si es debido a la convergencia a
una solución de naturaleza local. Ésta es una clara motivación para la utilización
de métodos que proporcionen más garantı́as de convergencia al óptimo global tanto
para resolver el problema de calibración como para resolver el problema de diseño
óptimo de experimentos.
La creciente demanda de los consumidores con respecto a la calidad de los ali-
mentos y el endurecimiento de las normas de seguridad, han motivado el desarrollo
de métodos de computación basados en modelos para la simulación, la optimización
y el control de técnicas de procesamiento de alimentos. Las aproximaciones basadas
en modelos son también un tema central en la biologı́a de sistemas ya que proporcio-
nan nuevos modos de analizar los datos procedentes de la genómica y la proteómica,
proporcionando un gran entendimiento sobre el lenguaje de las células y los orga-
nismos. Además, estas aproximaciones proporcionan estrategias sistemáticas para
cuestiones clave de la medicina y la industria farmacéutica y biotecnológica como,
por ejemplo, el desarrollo de fármacos teniendo en cuenta los efectos de posibles
nuevos medicamentos en rutas bioquı́micas y en la fisiologı́a.
En este trabajo se estudia el modelado y la identificación de una serie de bioproce-
sos. Relativos a la industria alimentaria, se han considerado procesos de conservación
basados en técnicas de secado y procesamiento térmico de alimentos. En relación a
la biologı́a de sistemas, se han considerado modelos de rutas bioquı́micas de gran
interés ası́ como la modelización de la cinética de la glucosa en pacientes diabéticos
que es un paso clave en el desarrollo del deseado “páncreas artificial”.
Para llevar a cabo estas tareas, se presentan nuevas metodologı́as de optimización
global que aumentan significativamente la eficiencia de las hasta ahora utilizadas
garantizando su robustez. Se hace también una revisión de los métodos de análisis de
sensibilidades ası́ como de los tipos de funciones de sensibilidad y de su aplicabilidad,
especialmente para cuantificar la importancia de los parámetros estableciendo un
ranking de los mismos. Además, se analizan las técnicas existentes para el estudio de
la identificabilidad y se presenta un programa desarrollado en Matlab° R
que, como se
explicará detalladamente a lo largo de este trabajo, permite automatizar las tareas de
análisis de identificabilidad, ranking de parámetros, calibración del modelo, cálculo
de intervalos de confianza y diseño óptimo de experimentos dinámicos.
Objetivos
xiii
xiv Objetivos
I Introducción 1
1. Modelos matemáticos 3
1.1. Desarrollo de modelos matemáticos . . . . . . . . . . . . . . . . . . . 4
1.2. Tipos de modelos matemáticos . . . . . . . . . . . . . . . . . . . . . . 5
II Metodologı́a 9
2. Estimación de parámetros 11
2.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.1. Caracterización del modelo . . . . . . . . . . . . . . . . . . . . 11
2.1.2. Datos experimentales . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.3. Funciones de coste . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2. Métodos de estimación . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.1. Métodos de valor inicial . . . . . . . . . . . . . . . . . . . . . 18
2.2.2. Método multiple shooting . . . . . . . . . . . . . . . . . . . . . 19
3. Análisis de sensibilidad 21
3.1. Métodos numéricos para el cálculo de sensibilidades locales . . . . . . 22
3.1.1. Aproximación por diferencias finitas . . . . . . . . . . . . . . . 22
3.1.2. Métodos directos . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.3. Método de la función de Green . . . . . . . . . . . . . . . . . 24
3.2. Tipos de funciones de sensibilidad . . . . . . . . . . . . . . . . . . . . 24
3.2.1. Función de sensibilidad absoluta . . . . . . . . . . . . . . . . . 25
3.2.2. Función de sensibilidad relativa . . . . . . . . . . . . . . . . . 25
3.2.3. Función de sensibilidad semirelativa . . . . . . . . . . . . . . . 26
3.3. Ranking de parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4. Análisis de identificabilidad 29
4.1. Identificabilidad estructural . . . . . . . . . . . . . . . . . . . . . . . 30
xv
xvi Índice general
5. Intervalos de confianza 37
5.1. Regiones de confianza exactas . . . . . . . . . . . . . . . . . . . . . . 37
5.2. Método basado en la FIM . . . . . . . . . . . . . . . . . . . . . . . . 38
5.3. Método basado en la matriz Hessiana . . . . . . . . . . . . . . . . . . 39
5.4. Métodos de Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.4.1. Jackknife . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.4.2. Bootstrap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7. Métodos de optimización 51
7.1. Métodos locales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
7.1.1. Métodos para problemas sin restricciones . . . . . . . . . . . . 53
7.1.2. Métodos para problemas con restricciones . . . . . . . . . . . 54
7.1.3. Métodos locales empleados . . . . . . . . . . . . . . . . . . . . 55
7.2. Métodos globales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
7.2.1. Métodos deterministas . . . . . . . . . . . . . . . . . . . . . . 56
7.2.2. Métodos estocásticos . . . . . . . . . . . . . . . . . . . . . . . 57
7.2.3. Métodos hı́bridos . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.3. Desarrollo de un método hı́brido secuencial . . . . . . . . . . . . . . . 62
7.3.1. Ajuste del método hı́brido secuencial . . . . . . . . . . . . . . 63
7.4. Método hı́brido paralelo sincrónico . . . . . . . . . . . . . . . . . . . 64
III Aplicaciones 77
9. Secado de alimentos 79
9.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
9.2. Modelo matemático . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
9.2.1. Transferencia de masa . . . . . . . . . . . . . . . . . . . . . . 81
9.2.2. Transferencia de energı́a . . . . . . . . . . . . . . . . . . . . . 82
9.3. Análisis de identificabilidad estructural . . . . . . . . . . . . . . . . . 83
9.4. Ranking de parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . 86
9.5. Estimación de parámetros . . . . . . . . . . . . . . . . . . . . . . . . 87
9.5.1. Caso 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
9.5.2. Caso 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
9.6. Identificabilidad a posteriori . . . . . . . . . . . . . . . . . . . . . . . 91
9.7. Intervalos de confianza . . . . . . . . . . . . . . . . . . . . . . . . . . 91
9.8. Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
IV Conclusiones 173
V Apéndices 181
VI Bibliografı́a 189
xxi
xxii Índice de tablas
xxiii
xxiv Índice de figuras
12.1. Esquema de reacción para la inhibición irreversible de la proteasa del HIV 119
12.2. Parámetros ordenados por orden decreciente de δ msqr . . . . . . . . . . . 121
12.3. Frecuencia de las soluciones de un SQP en modo multi-start . . . . . . . . 123
12.4. Curvas de convergencia de SRES, DE y SSm . . . . . . . . . . . . . . . . . 124
12.5. Datos experimentales versus valores predichos por el modelo . . . . . . . 125
12.6. Valores de los residuos en función del tiempo . . . . . . . . . . . . . . . 125
12.7. Matriz de correlación a posteriori . . . . . . . . . . . . . . . . . . . . . 126
Tipografı́a
Itálica Escalar
Negrita minúscula Vector
Negrita mayúscula Matriz
Abreviaturas
xxvii
xxviii Notación
Sı́mbolos
² Errores de observación (ruido)
C Matrix de covarianza
F Jacobiano paramétrico de un conjunto de ODEs
F1 , F2 y F3 Condiciones frontera de primer, segundo y tercer orden
J Función objetivo
Jmp Función de máxima probabilidad
Jmc Función de mı́nimos cuadrados
J Jacobiano de un conjunto de ODEs
M Estructura de un modelo
N Número total de datos experimentales
Np Número de parámetros a estimar
Nu Número de variables de control
Nx Número de variables de estado distribuidas
Ny Número de variables de estado concentradas
Nz Número de variables medidas
R Matriz de correlación
σ Desviación estándar (ruido de las medidas)
p Vector de los Np parámetros del modelo
p∗ Vector de parámetros verdaderos del proceso
p̂ Estimador de p asociado a óptimo local
p̌ Estimador de p asociado a óptimo global
t Variable temporal
u Vector de las Nu variables de control
x Vector de las Nx variables de estado distribuidas
xξ y xξξ Vectores de la primera y segunda derivada espacial de x
xt Derivada temporal de x
y Vector de las Ny variables de estado concentradas
ẏ Derivada temporal de y
z Vector de las Nz variables medidas en cada experimento
z̃ Vector de las medidas experimentales
ξ Vector de coordenadas espaciales
Parte I
Introducción
Capı́tulo 1
Modelos matemáticos
3
4 Capı́tulo 1. Modelos matemáticos
Cualquiera que sea el objetivo del modelo, éste debe ser formulado explı́citamente
ya que influenciará en gran medida el proceso de modelado. Además, el modelo
obtenido debe ser juzgado en base a la satisfacción de esos propósitos.
Objetivos
Definición
Diseño de experimentos
del marco
Conocimientos
experimentales
Análisis
a priori
Datos
Caracterización
de la estructura
Estimación de
parámetros
NO
Validación
SÍ
Modelo
Definición del marco en el que se establecen los lı́mites del sistema y las
variables de entrada y salida.
Caracterización de la estructura en donde se determina el tipo de modelo
a considerar (lineal-no lineal, continuo-discreto,...), su nivel de compleji-
dad y las relaciones funcionales entre las variables.
Estimación de parámetros que proporciona valores numéricos para las
constantes de las relaciones funcionales.
modelo y la información disponible (no tiene sentido concebir un modelo muy com-
plejo con muchos parámetros si los datos experimentales disponibles son escasos e
imprecisos).
Entre las distintas clasificaciones que se pueden realizar de los modelos en función
de su estructura matemática cabe destacar (Jeppsson, 1996):
Por otra parte, se dice que la estructura de un modelo es lineal en sus paráme-
tros si la salida satisface el principio de superposición con respecto a sus
parámetros, es decir, si:
1988). De este modo, por definición, cualquier sistema de ODEs tiene ı́ndice
cero. Ası́ mismo, los sistemas de DAEs de ı́ndice uno se comportan de modo
bastante parecido a los sistemas de ODEs pudiendo ser resueltos utilizando
métodos similares. Sin embargo, el comportamiento de los sistemas de ı́ndi-
ce superior es cualitativamente diferente y deben ser tratados con métodos
especı́ficos.
Muchos procesos biológicos están distribuidos no sólo en el tiempo sino tam-
bién en el espacio. Matemáticamente, las variables distribuidas en el espacio
pueden describirse mediante ecuaciones diferenciales parciales (PDEs) y los
modelos resultantes son los llamados modelos distribuidos. Muchas veces estas
ecuaciones están combinadas con ecuaciones diferenciales ordinarias y ecua-
ciones algebraicas dando lugar a sistemas de PDAEs.
Metodologı́a
Capı́tulo 2
Estimación de parámetros
11
12 Capı́tulo 2. Estimación de parámetros
Condiciones frontera:
- mixtas o de Robin:
En caso de que se desconozcan las condiciones iniciales del problema, éstas tam-
bién pueden ser estimadas considerándose, a efectos de la calibración, como paráme-
tros adicionales.
Para la resolución de los sistemas distribuidos que aparecen en este trabajo, se
ha empleado el método numérico de las lı́neas (Numerical Method of Lines, NMOL)
(Schiesser, 1991) que transforma el problema original de dimensión infinita en uno
de dimensión finita, es decir, en un conjunto de ODEs o de DAEs. Por este motivo,
de aquı́ en adelante se prescindirá del vector de variables x y de sus derivadas.
Estimador Bayesiano
La estimación Bayesiana considera las medidas y los parámetros del modelo
como variables aleatorias. Si se conoce la densidad de probabilidad a priori
πp (p) para la ocurrencia del vector de parámetros p y la densidad de probabi-
lidad condicional πz (z̃|p) del modelo para medir los valores z̃ para unos valores
dados de los parámetros p, la densidad de probabilidad de los parámetros para
valores dados de las medidas se puede escribir como:
πz (z̃|p)πp (p)
πp (p|z̃) = (2.9)
πz (z̃)
J(p) (2.20)
y que verifica:
- Función objetivo con forma de valle angosto donde es difı́cil progresar hacia
la solución.
18 Capı́tulo 2. Estimación de parámetros
- Malos valores iniciales para los parámetros, lo que hace necesario un gran
número de iteraciones en la optimización.
- Funciones del modelo mal escaladas, en particular, los valores de las medidas.
Debido a estos escollos, se debe prestar especial atención al método elegido para
llevar a cabo la identificación que, a grandes rasgos, pueden clasificarse en dos grupos:
los métodos de valor inicial (o single shooting) y el método de disparo múltiple (o
multiple shooting).
Integración de las
Definición estructura ecuaciones del modelo
del modelo (IVP)
¿Mínimo NO
Nuevos valores
de la función
para los parámetros
objetivo?
SÍ
Mejor estimador
para los parámetros
Valores iniciales
z (t)
8ª iteración
z (t)
Análisis de sensibilidad
21
22 Capı́tulo 3. Análisis de sensibilidad
Hay varios métodos numéricos para el cálculo de sensibilidades locales pero los
valores calculados deben ser idénticos dentro de la precisión del método empleado.
y Kramer, 1985) :
µ ¶ µ ¶
∂f (t) ∂f (t)
Ṡ(t) = S(t) + ; S(0) = S0 (3.4)
∂y p ∂p y
o, de forma matricial:
En esta ecuación, K se conoce como función de Green del núcleo y el método numéri-
co basado en la solución de la ecuación (3.8) se llama método de la función de Green
(Green Function Method, GFM).
El método Magnus analı́ticamente integrado (Analytically Integrated Magnus,
GFM/AIM) es una modificación del método de la función de Green más desarrolla-
do. En esta versión, la matriz K es aproximada por una matriz exponencial dismi-
nuyendo significativamente el esfuerzo de cálculo:
·Z t+∆t ¸
K(t + ∆t, t) = exp J(s)ds (3.9)
t
El método GFM/AIM es varias veces más rápido que otras versiones del método de
la función de Green.
Aplicando el método directo, el esfuerzo numérico aumenta linealmente con el
número de parámetros. En el caso de los métodos de la función de Green, el esfuerzo
numérico es proporcional al número de variables. Sin embargo, en la práctica, el
método GFM sólo es más rápido que el DDM cuando la relación entre el número
de parámetros y el de variables es muy elevada y el error numérico es más difı́cil de
controlar en este caso que utilizando el método DDM, que es mucho más simple.
una breve descripción de cada una de ellas y sus aplicaciones principales destacando
sus ventajas e inconvenientes.
Mientras que las funciones de sensibilidad semirelativas con respecto a las va-
riables de estado tienen la misma forma que las funciones de sensibilidad relativas
(y por lo tanto los mismos problemas de división por cero y sobrepesado), las fun-
ciones de sensibilidad semirelativas con respecto a los parámetros tienen la misma
forma que las funciones de sensibilidad absolutas (están únicamente multiplicadas
por los valores constantes de los parámetros) pero este reescalado permite hacer
comparaciones de los efectos de los distintos parámetros.
Cuando se utilizan funciones de sensibilidad, tanto relativas como semirelativas
con respecto a las variables, se puede definir un valor umbral ymin en el factor
premultiplicador de las ecuaciones (3.11) y (3.12) cuando éste es menor que el valor
de ymin . De este modo, se evitan los errores de sobrepesado cuando la trayectoria de
salida tiende a cero (Versyck, 2000) .
Ny N
1 1 XX
δjmabs = |S ij (tk )| (3.15)
Ny N i=1 k=1
Ny N
1 1 XX
δjmean = S ij (tk ) (3.16)
Ny N i=1 k=1
Análisis de identificabilidad
29
30 Capı́tulo 4. Análisis de identificabilidad
1 1
zi (p, O+ ) = a0i (p) + a1i (p)t + a2i (p)t2 + ... + ani (p)tn (4.4)
2 n!
siendo:
dj zi
aji (p) = ∀i = 1, ..., Nz ; j = 0, ...n (4.5)
dtj
Ya que el vector de medidas es único, todas sus derivadas son también únicas.
Por lo tanto una condición suficiente para que el modelo sea s.g.i. será:
Considerando los valores de los parámetros del conjunto p̂ como “valores ver-
daderos”, las matrices de sensibilidad de los estados medidos, Sz , de dimensión Nz
por Np , se calculan para un número suficientemente grande de puntos de tiempo N
donde: µ ¶
∂zi
Szij = (4.9)
∂pj z=z(t,p̂),p=p̂
La matriz G se construye entonces almacenando las matrices de sensibilidades
para estos puntos:
Sz (t1 )
Sz (t2 )
G= .. (4.10)
.
Sz (tN )
Finalmente la matriz de correlación de los parámetros (Mc ), de dimensión Np
por Np , se calcula como:
Mc = correlación(G) (4.11)
donde z̃i y zi (p) son vectores de N valores medidos y predicciones del modelo a los
tiempos ti (i = 1 a N ), respectivamente, y Qi es una matriz cuadrada proporcionada
por el usuario de coeficientes de peso.
4.3. Identificabilidad práctica o a posteriori 35
donde Vi representa la matriz de covarianza del error de las medidas (Qi se elige
tı́picamente como Vi −1 ).
Una consecuencia importante de la ecuación (4.13) es que para optimizar la
identificabilidad práctica (maximizar la diferencia entre J(p + δp) y J(p)) se tiene
que maximizar el término entre corchetes [·]. Este término se conoce como matriz
de información de Fisher (FIM) y expresa la cantidad de información de los datos
experimentales (Ljung, 1999):
N µ
X ¶T µ ¶
∂z ∂z
FIM = (ti ) Qi (ti ) (4.14)
i=1
∂p ∂p
Los términos ∂z/∂p son las funciones de sensibilidad que son de gran importancia
para la evaluación de la identificabilidad práctica ya que son el componente principal
de la matriz de información de Fisher.
Como se explicará en detalle en el capı́tulo 5, la matriz de información de Fisher
es también una aproximación de la inversa de la matriz de covarianza del error del
mejor estimador lineal no sesgado (Best Linear Unbiased Estimator, BLUE):
" N µ ¶T µ ¶#−1
X ∂z ∂z
C = FIM−1 = (ti ) Qi (ti ) (4.15)
i=1
∂p ∂p
Rij = 1, i = j, (4.17)
parámetros sean identificables de forma única incluso si la salida del modelo es muy
sensible a los cambios en los parámetros individuales.
El problema de analizar la identificabilidad práctica es similar al análisis de la
identificabilidad a priori local pero ahora los puntos de evaluación de las funciones
están limitados a los datos experimentales y éstos tienen error. Si las funciones de
sensibilidad presentan dependencia lineal en los puntos de los datos experimentales,
la matriz de covarianza se vuelve singular y el modelo no es identificable. Una FIM
singular indica la presencia de parámetros no identificables y correlaciones entre
parámetros superiores de 0.99 pueden dar lugar a una FIM singular.
En el presente estudio, se ha utilizado el estimador del número de condición de
una matriz de Matlab, rcond, para determinar si la FIM es singular:
1
rcond(FIM) = (4.18)
norm(FIM, 1)norm(FIM−1 , 1)
Si rcond(FIM) < 10ε, donde ε es la máxima precisión de Matlab en punto
flotante (2.2 10−16 ), la FIM se considera singular.
A pesar de que el análisis de sensibilidad basado en la FIM es una técnica muy
extendida, ésta implica una linealización de primer orden del modelo con respecto a
los parámetros lo que en algunos casos dará lugar a conclusiones erróneas (Petersen,
2000). En el caso de problemas muy no lineales, en esta aproximación se pierde
mucha información lo que podrı́a dar lugar a que los parámetros sean identificables
en la práctica aún cuando la FIM sea singular.
Por este motivo, el análisis basado en las regiones de confianza se presenta en este
trabajo como una una alternativa más robusta aunque mucho más costosa desde el
punto de vista computacional.
Intervalos de confianza
37
38 Capı́tulo 5. Intervalos de confianza
J(p̂) £ ¤−1
CJ (p̂) = J(p̂)T V−1 J(p̂) (5.5)
N − Np
donde el término J(p̂)/N − Np es una aproximación objetiva de la varianza residual
σ 2 y J es la matrix Jacobiana del modelo que puede escribirse en columna, tal que:
£ ¤
J = J1 |J2 |...|Jj |...|JNp (5.6)
siendo ¯
∂z ¯¯
Jj (p̂) = = Sj , j = 1, ..., Np (5.7)
∂pj ¯p̂
5.3. Método basado en la matriz Hessiana 39
Las columnas de la matrix Jacobiana J son las salidas de las funciones de sensibilidad
Sj = (∂z/∂pj ) con respecto a los parámetros.
Asumiendo que el ruido de las medidas no está correlacionado y que éste pre-
senta una distribución normal con media cero y varianza constante, CJ , dada por
la ecuación (5.5), es también la inversa de la FIM, definida como:
N µ
X ¶T µ ¶
∂z(ti ) ∂z(ti )
FIM = Vi−2 (5.8)
i=1
∂p ∂p
es decir,
CJ (p̂) = FIM−1 (5.9)
En este caso, la inversa de la FIM representa la matriz de covarianza del error del
estimador objetivo de varianza mı́nima de acuerdo con el teorema de Cràmer-Rao
(Ljung, 1999). Sustituyendo C de la ecuación (5.5) en la ecuación (5.3), se obtienen
los elipsoides de confianza aproximados:
n o
T −1 1−α
p : (p − p̂) CJ (p − p̂) ≤ Np FNp ,N −Np (5.10)
donde:
∂ 2 J(p̂)
H(p̂) = (5.13)
∂p∂pT
con los elipsoides de confianza dadas por:
n o
p : (p − p̂)T CH −1 (p − p̂)T ≤ Np FN1−α
p ,N −Np
(5.14)
40 Capı́tulo 5. Intervalos de confianza
1−(α/2)
p
σi = ±tN −Np Cii (5.15)
5.4.1. Jackknife
La mayor ventaja de la técnica jackknife estriba en su simplicidad. Sea N = Gh,
el vector de datos experimentales z̃ = (z̃1 , z̃2 , ..., z̃N ) se divide de la forma z̃ =
(z̃01 , z̃02 , ..., z̃0G ), siendo cada z̃0g (g = 1, 2, ..., G) un vector h × 1. Sea p̂ el vector de
parámetros estimados a partir de todos los datos z̃ y p̂(g) el mejor estimador de
parámetros para el modelo omitiendo el grupo z̃0g calculado iterativamente emplean-
do p̂ como punto inicial. Se definen entonces G pseudo-estimadores de la forma:
con media:
G
1 X
p̄Jk = p̃g (5.17)
G g=1
y matriz de varianza-covarianza:
G
£ ¤ 1 X
CJk = (CJkij ) = (p̃g − p̄Jk ) (p̃g − p̄Jk )T (5.18)
G − 1 g=1
Las propiedades de (5.19) y (5.20) fueron estudiadas por Dunkan (1978) para
N = 24 y distintos valores de h. Aunque tener un h lo más grande posible es
conveniente desde el punto de vista computacional, su estudio concluye claramente
que las omisiones de uno en uno (h = 1) son más recomendables, especialmente para
muestras pequeñas. A pesar de su fácil implementación, este método resulta menos
flexible y fiable que el método bootstrap (Walter y Pronzato, 1997).
5.4.2. Bootstrap
El método bootstrap (ver, por ejemplo, DiCiccio y Romano, 1988; Hinkley, 1988)
utiliza solamente los valores de los datos experimentales z̃ y el modelo M (p̂) asu-
miendo que los errores son variables aleatorias con idéntica distribución pero sin
especificar. Se asume por ejemplo que
- ¿Cómo manipular? Se trata de definir cuáles son las posibles variables ma-
nipuladas o controles y el tipo de manipulaciones, es decir, la variación de
los controles a lo largo del experimento. En algunos casos también se podrán
modificar las condiciones iniciales de ciertas variables.
43
44 Capı́tulo 6. Diseño óptimo de experimentos
Criterio A Criterio E
p2
Criterio D
p1
¡ ¢
Criterio A: min traza FIM−1
Este criterio se centra en la minimización de la traza y por lo tanto de la
suma de los autovalores de la matriz de covarianza, es decir, el cuadrado de
la longitud de los ejes de los elipsoides de confianza. Esto es equivalente a
minimizar la media aritmética de los errores de los parámetros. Nótese que
este criterio se basa en la inversión de la FIM por lo que producirá errores
numéricos en el caso de que la FIM sea singular o esté mal condicionada.
h (y, z, u, v, p, t) = 0 (6.4)
g (y, z, u, v, p, t) ≤ 0 (6.5)
48 Capı́tulo 6. Diseño óptimo de experimentos
Métodos de optimización
51
52 Capı́tulo 7. Métodos de optimización
dos globales lo que se pretende es encontrar el valor de p̌ tal que J(p̌) < J(p) para
todos los posibles valores de p. Nótese que, en el caso de los problemas convexos, el
óptimo local coincide con el global.
Por otra parte, los métodos estocásticos calculan las direcciones de búsqueda
empleando secuencias pseudo-aleatorias con un importante componente heurı́stico y
valores previos de la función objetivo, sin hacer uso de información sobre la estruc-
tura del problema. En general, son de carácter global aunque debido a su naturaleza
aleatoria no se puede garantizar la convergencia al óptimo global con total certe-
za. Sin embargo, muchos de estos métodos disponen de pruebas de convergencia
asintóticas y a menudo son estrategias muy eficientes localizando las proximidades
del óptimo global.
A su vez, la mayor parte de los métodos hı́bridos combinan estrategias estocásti-
cas y deterministas para superar las limitaciones inherentes a cada una de ellas
mientras que potencian sus puntos fuertes. No obstante, un método hı́brido es todo
aquel que resulta de la combinación de dos o más algoritmos por lo que también
se pueden encontrar estrategias que combinan varios métodos estocásticos o varios
métodos deterministas. La clave para la obtención de métodos hı́bridos robustos y
eficaces radica en la elección de los métodos a combinar y en el modo de estructurar
dicha combinación.
primer orden y aquellos que requieren las derivadas de primer y segundo orden
se llaman métodos de segundo orden. El primer método de este grupo, propues-
to por Cauchy en el siglo XIX, es el método del gradiente clásico que utiliza
el negativo del gradiente como dirección de búsqueda para la minimización.
En este método la dirección de búsqueda puede interpretarse como ortogo-
nal a una aproximación lineal tangente a la función objetivo en ese punto. Los
métodos de segundo orden tienen también en cuenta la curvatura de la función
para lo que deberá considerarse la matriz Hessiana de la misma. Entre estos
métodos destacan el método de Newton y los populares métodos quasi-Newton
que tratan de aproximar la dirección de Newton asintóticamente.
Estas técnicas han sido originalmente desarrolladas para resolver problemas sin
restricciones pero se adaptan fácilmente a la resolución de problemas con restriccio-
nes mediante el uso de funciones de penalización (Balsa-Canto, 2001).
clsSolve: este algoritmo forma parte del entorno de optimización Tomlab (Holms-
tröm, 2004) y resuelve problemas de estimación de parámetros no lineales,
dispersos o densos, manejando explı́citamente igualdades y desigualdades li-
neales y lı́mites simples en las variables.
npsol: desarrollado por Stanford Systems Optimization Laboratory (ver Gill et al.,
1998), se considera como el estado actual para la resolución de problemas de
optimización no lineal densos.
snopt: desarrollado por Stanford Systems Optimization Laboratory (ver Gill et al.,
2002), se considera el método más puntero para la resolución de problemas
grandes y dispersos de optimización no lineal.
En la actualidad, los tipos de métodos estocásticos más populares son los algorit-
mos genéticos (GAs) y los de templado simulado (SA). Sin embargo, según nuestra
propia experiencia y como muchos autores han señalado en los últimos años (Banga
et al., 2003; Moles et al., 2003b), los GAs y SA, diseñados en primera instancia
para problemas combinatorios (con variables enteras), no son normalmente los al-
goritmos más eficientes y robustos para la optimización global en variables reales.
La elección de un método para un cierto tipo de problemas de GO es complicada
y la bibliografı́a al respecto está muy fragmentada, por lo que en muchos casos la
elección está basada más en la predilección del autor por un tipo de técnicas que en
criterios racionales (Preux y Talbi, 1999). Aunque éste continúa siendo un tema de
gran debate, sin embargo, Wolpert y MacReady (1997) han demostrado a través del
teorema “no hay comida gratis” (No Free Lunch, NFL) que no existe ningún méto-
do que pueda ser considerado mejor que los demás para la resolución de problemas
generales de optimización global de estructura desconocida.
No obstante, para el caso de optimización global de NLOs y para el caso par-
ticular de estimación de parámetros, diferentes trabajos recientes (Balsa-Canto et
al., 1998; Moles et al., 2003a; Moles et al., 2003b) indican que ciertos métodos
estocásticos simples, en concreto el método Evolución Diferencial (Differential Evo-
lution, DE) (Storn y Price, 1997) y ciertas estrategias evolutivas como la Estrategia
Evolutiva con Ranking Estocástico (Stochastic Ranking Evolution Strategy, SRES)
desarrollada por Runarsson y Yao (2000, 2005) presentan un comportamiento mejor
en términos de eficiencia y robustez. Además estas estrategias escalan bien con el
tamaño del problema y permiten ser paralelizados muy fácilmente, lo que signifi-
ca que los problemas de mediana o gran escala podrán ser resueltos en un tiempo
computacional razonable.
Por estos motivos, estos dos serán los métodos estocásticos considerados en el
presente estudio y servirán también de base para la creación de métodos hı́bridos
que mejoren su eficiencia manteniendo su robustez:
Hibridación secuencial
Hibridación paralela
Dependiendo del tamaño del problema a tratar, será conveniente considerar im-
plementaciones paralelas de los métodos a utilizar. Sin embargo, se debe distinguir
entre la mera paralelización de un algoritmo secuencial y una implementación de
un algoritmo paralelo. Una implementación en paralelo de un algoritmo trata de
mantener la esencia de la búsqueda secuencial mientras que los algoritmos hı́bridos
paralelos son una sub-clase de algoritmos estocásticos que se comportan de modo
diferente a los secuenciales. Dentro de estos métodos se debe distinguir entre la
hibridación paralela sincrónica y asincrónica.
La idea fundamental de la hibridación paralela sincrónica es la utilización de un
algoritmo como un operador de otro. Un ejemplo de este tipo de hibridación serı́a
emplear un algoritmo evolutivo como estrategia de búsqueda principal en donde, los
tradicionales operadores matemáticos utilizados en las etapas de recombinación y/o
mutación, son sustituidos por otro algoritmo que puede ser desde un método local
tradicional hasta un método tabú o incluso un algoritmo de templado simulado SA.
Este tipo de hibridación se denomina hibridación paralela sincrónica porque en ella
los diferentes algoritmos están sincronizados con precisión.
62 Capı́tulo 7. Métodos de optimización
SÍ SÍ
¿Llamada ¿Pasa
método filtros?
local?
NO
NO
Regeneración Actualización Cálculo
RefSet RefSet resultados
método local
SÍ
¿Elementos no
combinados?
NO
NO ¿Satisfacción
criterio de
parada?
SÍ
Resultados
para mejorar las soluciones, tanto del conjunto de referencia como las combi-
nadas antes de estudiar su inclusión en el conjunto de referencia.
Cuando ya se han hecho todas las combinaciones entre las soluciones del Ref-
Set, el algoritmo puede parar o continuar mediante la reconstrucción parcial
del conjunto de soluciones de elite. Se ha implementado una nueva estrategia
para reconstruir este conjunto, basada en direcciones de búsqueda ortogonales.
El usuario puede elegir entre un amplio número de métodos locales SQP como
fmincon (The MathWorks Inc.), solnp (Ye, 1987), npsol (Gill et al, 1998),
snopt (Gill et al, 2002), métodos directos como NOMADm (Abramson, 2002)
para casos con datos con mucho ruido, y otros especı́ficamente diseñados para
problemas de estimación de parámetros como n2fb/dn2fb de Dennis et al.
(1981).
Capı́tulo 8
69
70 Capı́tulo 8. GOSBio: entorno para modelado e identificación
Diseño Cálculo
experimentos sensibilidades
NO
SÍ NO
NO SÍ
SÍ
Cómputo
FIM
Paso 9: Validación del modelo: a pesar de que no existe ningún método que pueda
garantizar la validez de un modelo con total certeza, es necesario investigar el
comportamiento del mismo de modo que, si supera todas las pruebas de invali-
dación a las que es sometido, se pueda considerar satisfactorio. Este programa
permite realizar distintas pruebas de invalidación:
opt solver.name:
opt solver.u guess/u min/u max: valores iniciales y lı́mites inferiores y su-
periores para los controles en OED.
opt solver.t f/tf min/tf max: valor inicial y lı́mite inferior y superior para
el tiempo final en OED.
8.3.2. Figuras
Además se generan una serie de figuras que se detallan a continuación. Éstas se
muestran en pantalla y son posteriormente almacenadas en la carpeta de resultados
como figuras de Matlab (.fig) y en formato PostScript encapsulado (.eps).
fit plot: representación de los valores predichos por el modelo para cada
uno de los estados medidos (linea continua) y de los datos experimentales
correspondientes (marcador) a lo largo del tiempo.
Aplicaciones
Capı́tulo 9
Secado de alimentos
9.1. Introducción
La creciente demanda de los consumidores con respecto a la calidad de los ali-
mentos y el endurecimiento de las normas de seguridad, han motivado el desarrollo
de métodos de computación basados en modelos para la simulación, la optimización
y el control de técnicas para su procesamiento (Datta, 1998; Banga et al., 2003).
El modelado de secado por aire de alimentos ha recibido una gran atención
durante las últimas décadas ya que es uno de los métodos de preservación más
importantes. Los primeros modelos matemáticos simples de este proceso aparecieron
a finales de los años 70. El desarrollo de nuevas técnicas numéricas y el incremento
de las capacidades computacionales han permitido el aumento de la complejidad
de modelos posteriores. Muchos autores han revisado los distintos avances, ver p.ej.
Bruin y Luyben (1980), Jayaraman y Das Gupta (1992) o Waananen et al. (1993) o,
más recientemente, Ruiz-Lopez et al. (2004) en el contexto de modelado de procesos
y Banga y Singh (1994) y los trabajos allı́ citados, en el contexto de optimización.
La combinación de las leyes fı́sicas de transferencia de masa y de energı́a con las
propiedades fı́sicas del producto alimentario, permiten la predicción de la variación
a lo largo del tiempo de las variables de estado relevantes (contenido de humedad y
temperatura), sujetas a diferentes condiciones de secado. Aunque los modelos más
simples asumen que la contracción del alimento durante el proceso es despreciable y
que las propiedades de transporte son constantes, se ha ilustrado experimentalmente
que estas suposiciones no son realistas (Balaban, 1989; Park, 1998; Simal et al., 1998)
y que, por lo tanto, cualquier modelo riguroso deberı́a considerar estos efectos (ver
una revisón más extensa en Mayor y Sereno (2004)).
En la bibliografı́a se han propuesto muchos modelos (la mayorı́a empı́ricos) para
diferentes tipos de alimentos. La mayor parte de ellos son muy no lineales, frecuen-
79
80 Capı́tulo 9. Secado de alimentos
corriente de aire
Tdb
Ts , m
L3 = 0.2 cm
L2 = 3.5 cm
L1 = 3.4 cm
dm
= ∇ (D∇m) (9.1)
dt
Debido al pequeño grosor de la lámina en comparación con las otras dimensiones,
ésta puede ser considerada como un sistema semi-infinito en donde el contenido de
humedad depende sólo de la posición con respecto a la dimensión menor. Además,
con objeto de tener en cuenta el efecto de la contracción, se asume que la difusividad
D es una función no lineal de ambos, el contenido de humedad y la temperatura,
por lo que la ecuación (9.1) resulta:
µ ¶ µ ¶2
dm ∂ 2m ∂D ∂m
=D + (9.2)
dt ∂x2 ∂m ∂x
82 Capı́tulo 9. Secado de alimentos
ρ0 L1 L2 L3
ms = (9.7)
mavg,0 + 1
λw = α1 − α2 Ts (9.9)
A0 = 2 (L1 L2 + L1 L3 + L2 L3 ) (9.11)
ṁ1 = 0 (9.12)
µ ¶ µ ¶µ ¶
mi−1 − 2mi + mi+1 Di+1 − Di−1 mi+1 − mi−1
ṁi = Di + (9.13)
2dx 2dx 2dx
µ ¶
mnx−1 − mnx
ṁnx = Dnx (9.14)
dx
nx
1X
ṁavg = ṁi (9.15)
L i=1
µ ¶
1
Ṫs = (A0 (p1 mavg + p2 ) (Tdb − Ts ) (9.16)
ms Cps + ms Cpw mavg
+ms (α1 − α2 Ts ) ṁavg )
a01 = mavg,0
a02 = Ts,0
1 1
a11 = − dxDm,0 m0 + dx2 (Dm,0 − Dm1 ,0 ) m0
2 4
λw,0 a11 A0 (a01 p1 + p2 ) (Tdb,0 − a02 )
a12 = +
Cps + Cpw a01 (Cps + Cpw a01 ) ms
2
µ µ
dx m0 a12 (ED,0 − b4 Dm1 ,0 ) Dm,0 dx m0 (b2 + b3 log (Dref,0 ))
a21 = 2
− a11 1 −
4Ra02 2 1 + b3 m0
¶ ¶
m0 (−b5 + b6 ED,0 ) (log (Dref,0 ) − log (Dm,0 )) m0 ED,0 a12
+ +
ED,0 (1 + b6 m0 ) Ra202
³ ´
A0 (a01 p1 + p2 ) Ṫdb,0 − a12 + A0 p1 a11 (Tdb,0 − a02 ) + ms a21 λw,0
a22 =
(Cps + Cpw a01 ) ms
¡ ¢
a11 −a12 (α1 + Cpw ) − 12 dxDm,0 λw,0
+
Cps + Cpw a01
donde el sub-ı́ndice “,0” se refiere al valor de la magnitud en t = 0, tal que Dm,0 =
D (m0 ) y Dm1 ,0 = D (m1 (t = 0)).
La complejidad del sistema hace imposible el estudio de la identificabilidad es-
tructural global de todo el conjunto de parámetros p = [b1 b2 b3 b4 b5 b6 p1 p2 ]T . Sin
embargo, el primer coeficiente de Taylor permite extraer algunas conclusiones sobre
la identificabilidad estructural local para algunos subconjuntos:
i. Fijando todos los bi y utilizando a12 y a21 se puede obtener una solución
única para los parámetros relacionados con la transferencia de energı́a, p1 y
p2 , siempre que Tdb,0 sea diferente de Ts,0 .
Por lo tanto pk = [p1 p2 ]T son s.l.i.
18
6.755
16
6.75
14
6.745
12
b6
b3
10 6.74
8 6.735
6 6.73
4
6.725
6.4 6.5 6.6 6.7 6.8 6.9 7 7.1 8.71 8.715 8.72 8.725 8.73 8.735 8.74 8.745 8.75
b3 p1 −4
x 10
Figura 9.2: Lı́neas de contorno para Figura 9.3: Lı́neas de contorno para
los parámetros b3 y b6 los parámetros p1 y b3
iii. Las condiciones suficientes (ecuación 4.6) aplicadas a a22 revelan que algunos
conjuntos de cuatro parámetros son s.l.i. A modo ilustrativo, aquı́ se consi-
deró el caso pk = [b1 b4 p1 p2 ]T .
³h i´
De la ecuación a22 ([b1 b4 p1 p2 ]) = a22 b̂1 b̂4 p̂1 p̂2 se concluye que:
p1 = p̂1 (9.17)
a01 p1 + p2 = a01 p̂1 + p̂2 (9.18)
µ ¶ µ ¶
1 1 1 1
−b1 + b4 − + = −b̂1 + b̂4 − + (9.19)
RTs,0 RTref RTs,0 RTref
µ ¶
b1 b4 1 1 b̂1
− + + =− (9.20)
1 + b3 m0 1 + b6 m0 RTs,0 RTref 1 + b̂3 m0
µ ¶
b̂4 1 1
+ − +
1 + b6 m0 RTs,0 RTref
h i
[b1 b4 p1 p2 ] = b̂1 b̂4 p̂1 p̂2
100.25
0.01874
100.2
0.01873
100.15
0.01872 100.1
0.01871 100.05
b4
p2
0.01870 100
0.01869 99.95
0.01868 99.9
0.01867 99.85
99.8
0.01866
99.75
8.71 8.715 8.72 8.725 8.73 8.735 8.74 8.745 8.75 34.12 34.14 34.16 34.18 34.2 34.22 34.24 34.26 34.28
p1 x 10
−4 b1
Figura 9.4: Lı́neas de contorno para Figura 9.5: Lı́neas de contorno para
los parámetros p1 y p2 los parámetros b1 y b4
4
dmsqr
d
3 mabs
dmean
dmax
2 dmin
Valor del criterio
−1
−2
−3
−4
b3 b2 b1 b6 p2 b4 b5 p1
Parámetros
Estos resultados reflejan grandes diferencias en los valores de δ msqr para los
distintos parámetros lo que indica que la salida del modelo es considerablemente
sensible a algunos de ellos y muy poco sensible a otros. Ası́ se ve, por ejemplo, que
el parámetro ante el cual el modelo presenta una mayor sensibilidad es b3 y el de
menor influencia es p1 . Las lı́neas de contorno representadas en la Figura 9.3 también
apuntan en esta dirección mostrando cómo un cambio en el valor de b3 produce una
gran variación en el valor de la función objetivo mientras que un cambio en p1
produce un efecto mucho menor.
Las pequeñas diferencias entre δ msqr y δ mabs , indican que no existe una gran va-
riabilidad en las sensibilidades (Sj ) ni valores alejados (outliers). Una comparación
de δ max y δ min indica que todos los parámetros presentan sensibilidades tanto posi-
tivas como negativas siendo el efecto global negativo para tres de los parámetros y
positivo para los demás, como puede verse por el signo de δ mean .
donde wij corresponde a los diferentes pesos considerados con objeto de normalizar
la contribución de cada término:
9.5.1. Caso 1
9.5.2. Caso 2
Aquı́ se consideró el problema de estimación del conjunto pk = [b1 b4 p1 p2 ]T
que, como se ha demostrado, es estructuralmente identificable y multimodal. Para
acercarse más a las condiciones experimentales reales, se consideró un error relativo
gaussiano de un 5 % en las medidas pseudo-experimentales.
Con objeto de enfatizar la necesidad de métodos globales de optimización, en
primer lugar se intentó resolver el problema utilizando un método de optimización
local de tipo SQP (solnp) en modo multi-start. El histograma de la Figura 9.7
representa la frecuencia de las soluciones mostrando que la mayorı́a de ellas están
lejos del óptimo global (soluciones locales).
40
35
30
25
Frecuencia
20
15
10
0
0 0.5 1 1.5 2 2.5 3
4
Función Objetivo x 10
Por otra parte, el uso de métodos globales (SRES, DE y SSm) dio lugar solucio-
nes óptimas globales siendo el valor del vector de parámetros encontrado próximo
al vector nominal. El algoritmo SSm convergió casi dos órdenes de magnitud más
rápido que los demás como se ilustra en la Figura 9.8 con las curvas de convergencia
(evolución del valor de la función objetivo con respecto al tiempo de computación).
4
10
DE
SRES
SSm
Función objetivo
3
10
2
10
Las Figuras 9.9 y 9.10 muestran una comparación entre los valores predichos a
partir del mejor vector de decisión y los datos pseudo-experimentales correspondien-
tes a la temperatura de la lámina (Ts ) y al contenido medio de humedad (mavg ).
Los parámetros estimados mediante los métodos de optimización global permiten
reproducir adecuadamente los datos pseudo-experimentales.
100 1
Tdb=100
90 0.9
T =85
db
0.8
80
T =75 0.7
db
70
Ts (ºC)
mavg
0.6
Tdb=65
60
0.5
T =55
50 db
0.4 Tdb=55
T =65
db
40 0.3 T =75
db
T =85
db
T =100
30 0.2 db
0 1000 2000 3000 4000 5000 6000 7000 0 1000 2000 3000 4000 5000 6000 7000
Tiempo (s) Tiempo (s)
Figura 9.9: Valores predichos versus Figura 9.10: Valores predichos versus
datos experimentales para Ts datos experimentales para mavg
9.6. Identificabilidad a posteriori 91
0.8
b4
0.6
0.4
b1 0.2
−0.2
p2
−0.4
−0.6
p1 −0.8
−1
p1 p2 b1 b4
resultados del análisis de correlación (R3,4 = 0.82). Sin embargo, la baja correlación
entre los parámetros p2 y b4 (R2,4 = 0.61) se traduce en una región de confianza más
esférica.
120 120
115 115
110 110
b4
105
b4
105
100 100
95 95
90 90
0.0125 0.013 0.0135 0.014 0.0145 0.015 0.0155 0.016 0.0165 33.8 34 34.2 34.4 34.6 34.8
p2 b
1
9.8. Conclusiones
En este capı́tulo se presentó un análisis detallado del problema de estimación
de parámetros relativo al secado por aire de alimentos que revela la necesidad de
llevar a cabo un procedimiento en dos pasos para garantizar que la solución alcan-
zada sea capaz de reproducir el comportamiento del sistema en un amplio rango de
condiciones de operación.
En una primera fase, se estudió la identificabilidad estructural del modelo por
medio de la aproximación de Taylor lo cual permitió concluir que sólo un subconjunto
de parámetros relacionados con la contracción y la transferencia de energı́a pueden
9.8. Conclusiones 93
Procesamiento térmico de
alimentos
10.1. Introducción
El procesamiento térmico es una de las operaciones más importantes para la
conservación de alimentos. La idea subyacente es procesar el producto alimentario a
una temperatura elevada durante un cierto periodo de tiempo con objeto de reducir
los microorganismos dañinos, garantizando la seguridad alimentaria y aumentando
su tiempo de conservación. Sin embargo, los cambios que experimentan los alimentos
durante el procesamiento térmico dan lugar, normalmente, a un empeoramiento de
su calidad. Resulta, por lo tanto, esencial procesar el alimento de manera que las
propiedades sensoriales y nutricionales se mantengan en los niveles más altos posibles
garantizando siempre la seguridad.
A este respecto, en las últimas décadas se han propuesto varias aproximaciones
para el diseño y/o optimización de este tipo de procesos (ver revisión por Banga
et al., 2003). Las metodologı́as más satisfactorias son aquellas que utilizan técnicas
computacionales basadas en modelos, especialmente aquellas que consisten en el
uso de modelos de principios fundamentales (Banga et al., 1991; Durance, 1997;
Ramaswamy et al., 1997; Balsa-Canto et al, 2002; Garcia et al. 2006).
Estos modelos matemáticos combinan las leyes de conservación con las propie-
dades fı́sicas de los productos, modelos de degradación de la calidad y cinéticas de
los microorganismos, permitiendo la predicción de las propiedades organolépticas y
la seguridad del producto final, sujetas a diferentes condiciones de procesamiento.
La degradación de la calidad se describe generalmente utilizando cinéticas de
cero, primer y segundo orden, y normalmente se asume que la dependencia de la
velocidad de degradación con la temperatura sigue la ecuación de Arrhenius (Sa-
95
96 Capı́tulo 10. Procesamiento térmico de alimentos
dC
− = f (T ) (10.2)
dt
donde C representa el factor de calidad considerado.
Con objeto de resolver las ecuaciones (10.1-10.2), deben imponerse condiciones
iniciales y frontera adecuadas en el dominio.
T (r, z, 0) = T0 (10.8)
Para este caso particular el estado medido corresponde a la retención media de
nutrientes que se calcula asumiendo una cinética de primer orden (10.2):
Z VT µ Z tf µ ¶ ¶
1 − ln 10 T (r, z, t) − TN,ref
y= exp exp ln 10 dt dV (10.9)
VT 0 DN,ref 0 ZN,ref
Nótese que este valor se mide a tiempo final por lo que sólo se dispone de una
medida por cada experimento.
Para transformar la ecuación diferencial parcial original (10.3) en un conjunto
de ecuaciones diferenciales ordinarias, se empleó el método numérico de las lı́neas
(NMOL) (Schiesser, 1991). El sistema de ODEs resultante, combinado con las ecua-
ciones cinéticas (una ecuación por posición espacial) se resolvió posteriormente con
ODESSA (Leis y Kramer, 1988) para permitir el cálculo de las sensibilidades pa-
ramétricas.
2
d
msqr
d
mabs
dmean
1.5
dmax
d
min
0.5
−0.5
p1 p2
Parámetros
Se consideró una sola medida por experimento a tiempo final y que éstas están
sujetas a un ruido gaussiano del 3 %.
3
Frecuencia
0
0.8 1 1.2 1.4 1.6 1.8 2 2.2
5
Funcion Objetivo x 10
DE
1
−10 SRES
SSm
2
−10
Función objetivo
3
−10
4
−10
5
−10
6
−10
0 5000 10000 15000
Tiempo CPU (s)
La Tabla 10.2 presenta los valores óptimos obtenidos por SSm mediante la ma-
ximización del criterio D considerando cinco, seis y ocho experimentos ası́ como los
valores del criterio E modificado para estos diseños.
5
4 x 10 4
3.5
3
3
E modificado
Criterio D
2 2.5
2
1
1.5
0 1
1 2 3 4 5 6
nexp
140 140
Experimentos: 2, 4, 5
Experimentos: 1, 3, 6
120 120
Temperatura (ºC)
Temperatura (ºC)
100 100
80 80
60 60
40 40
20 20
0 2000 4000 6000 8000 0 2000 4000 6000 8000
Tiempo (s) Tiempo (s)
Figura 10.5: Perfiles óptimos para los Figura 10.6: Perfiles óptimos para los
experimentos 1, 3 y 6 experimentos 2, 4 y 5
104 Capı́tulo 10. Procesamiento térmico de alimentos
1 140
1 140
Experimentos: 1, 3, 6
120 Experimentos: 2, 4, 5
120
0.8 0.8
100 100
T0 (º C)
T0
T0 (º C)
ret N
0.6 80
ret N
0.6 T0 80
retN 60 60
0.4 0.4
40 ret N 40
0.2 20 0.2 20
0 2000 4000 6000 8000 0 2000 4000 6000 8000
Tiempo (s) Tiempo (s)
0.8
0.6
Z N,ref
0.4
0.2
0.2
0.4
DN,ref
0.6
0.8
DN,ref Z N,ref
∗ ∗
DN,ref /DN,ref ZN,ref /ZN,ref
nexp Valor óptimo Int conf (95 %) Valor óptimo Int conf (95 %)
5 exp. 1.0 1.95e-2 1.0 2.02e-2
6 exp. 1.0 1.58e-2 1.0 1.58e-2
8 exp. 1.0 1.42e-2 1.0 1.44e-2
1.020
1.015
1.010
1.005
Zref/Z*ref
1.000
0.995
0.990
0.985
0.980
0.980 0.985 0.990 0.995 1.000 1.005 1.010 1.015 1.020
*
Dref /Dref
6
5
4
log(Jmc)
3
2
1
0
1.5 1.5
1 1
0.5 0.5 *
Dref /Dref
Zref /Z*
ref
10.8. Conclusiones
En este capı́tulo se consideró el diseño óptimo de experimentos para la estimación
de los parámetros cinéticos relativos al procesamiento térmico de bioproductos. El
empleo de un método local en modo multi-start detectó la presencia de óptimos loca-
les por lo que varios métodos estocásticos de optimización global fueron empleados
para la resolución del problema, siendo el método SSm el más rápido en converger.
Además, se ilustró el efecto del número de experimentos en el valor de distintos
criterios escalares de la FIM. De este modo, se ve como el valor del criterio D aumen-
ta con la cantidad de información a la vez que disminuyen los intervalos de confianza
para los parámetros. El empleo del criterio D dio lugar a diseños que consisten en
la repetición de dos tipos de experimentos: uno caracterizado por un procesamien-
to de corta duración a alta temperatura y otro por un procesamiento largo a baja
temperatura. Este resultado confirma la idea, ya apuntada por otros autores, de
que el criterio D tiende a repetir experimentos. Fisher (1935) ya señala las ventajas
de la repetición de experimentos argumentando que, si todos los experimentos son
diferentes, un error en uno de ellos disminuirá significativamente el rendimiento glo-
bal mientras que, si el experimento se repite varias veces, podrá probarse su validez
mediante la predominancia de las realizaciones con éxito.
Los resultados obtenidos demuestran que mediante un esquema experimental
óptimo no sólo se reducen los problemas de identificabilidad sino que además se
reducen los intervalos de confianza de los parámetros a la vez que disminuye el
esfuerzo experimental requerido, hasta un 50 % con respecto a las aproximaciones
tradicionales para el caso considerado.
Capı́tulo 11
11.1. Introducción
Este problema consiste en estimar cinco constantes de reacción (p1 , ..., p5 ) de
un sistema de reacción complejo estudiado originalmente por Box et al. (1973), que
forma parte de COPS (Collection of large-scale Constrained Optimization ProblemS)
(Dolan et al., 2004). La Figura 11.1 representa el esquema de reacción propuesto
para esta reacción quı́mica homogénea que describe la isomerización térmica del α-
pineno (y1 ) a dipenteno (y2 ) y allo-ocimeno (y3 ) que a su vez se convierte en α- y
β-pironeno (y4 ) y en un dı́mero (y5 ).
Este proceso fue estudiado por Fuguitt y Hawkins (1947), que proporcionaron
las concentraciones del reactante y de los cuatro productos en ocho intervalos de
tiempo (z̃ji ). Si los órdenes de las reacciones quı́micas son conocidos, se pueden
derivar modelos matemáticos que den las concentraciones de las distintas especies
109
110 Capı́tulo 11. Isomerización del α-pineno
dy1
= −(p1 + p2 )y1 (11.1)
dt
dy2
= p1 y1 (11.2)
dt
dy3
= p2 y1 − (p3 + p4 )y3 + p5 y5 (11.3)
dt
dy4
= p3 y3 (11.4)
dt
dy5
= −p4 y3 + p5 y5 (11.5)
dt
a01 = y1,0
a02 = y2,0
a03 = y3,0
a04 = y4,0
a05 = y5,0
a11 = −(p1 + p2 )y1
a12 = p1 y1
a13 = p2 y1 − (p3 + p4 )y3 + p5 y5
a14 = p3 y3
a15 = −p4 y3 + p5 y5
a11 ([p1 , p2 ]) = a11 ([pˆ1 , pˆ2 ]) =⇒ −(p1 + p2 )y1 = −(pˆ1 + pˆ2 )y1 (11.6)
a12 ([p1 ]) = a12 ([pˆ1 ]) =⇒ p1 y1 = pˆ1 y1 (11.7)
y esto es cierto sólo en el caso de que p4 = pˆ4 y p5 = pˆ5 siempre que la relación
entre y3 e y5 no sea constante, lo cual viene dado por el modelo de ODEs.
En este caso los coeficientes de Taylor de primer orden son suficientes para de-
mostrar que todos los parámetros del modelo son estructuralmente globalmente
identificables (s.g.i.).
2.5
dmsqr
2 d
mabs
dmean
1.5 dmax
dmin
1
Valor del criterio
0.5
−0.5
−1
−1.5
−2
−2.5
p1 p2 p4 p3 p5
Parámetros
5 X
X 8
J(p) = (zj (p, ti ) − z̃ji )2 (11.10)
j=1 i=1
Los métodos utilizados requieren un lı́mite superior e inferior ası́ como un punto
inicial para los parámetros. El lı́mite inferior para los cinco parámetros viene dado
por consideraciones fı́sicas, pi ≥ 0, y el lı́mite superior se considera pi ≤ 1, muy lejos
de la mejor solución conocida p1 = 5.93e-5, p2 = 2.96e-5, p3 = 2.05e-5, p4 = 27.5e-5,
p5 = 4.00e-5). Como punto inicial se elige pi = 0.5.
114 Capı́tulo 11. Isomerización del α-pineno
80
70
60
Frecuencia
50
40
30
20
10
0
3 3.2 3.4 3.6 3.8 4 4.2 4.4
Función objetivo 4
x 10
La Figura 11.4 muestra claramente que SSm convergió siempre a la solución global
en un tiempo de computación corto mientras que otros métodos de optimización
global fallaron o convergieron en un tiempo computacional mucho mayor. Con objeto
de favorecer la visualización, la curva correspondiente a SSm se representa en ejes
diferentes, ya que SRES y DE quedaron atrapados en soluciones locales cerca del
punto inicial mientras que SSm convergió al óptimo global lejos del primer valor.
Asimismo, la Figura 11.5 muestra una comparación entre los valores predichos
a partir del mejor vector de parámetros obtenido con SSm (linea continua) y los
datos experimentales para estas especies proporcionados por Fuguitt y Hawkins
(1947) (sı́mbolos), correspondientes a la concentración del reactante y de los cuatro
productos. Los parámetros estimados permiten reproducir los datos experimentales
y, como puede verse en la Figura 11.6, no existe correlación entre los residuos y el
tiempo lo que indica que el error de los datos experimentales es homocedástico.
11.6. Identificabilidad a posteriori 115
50000
SRES
Función objetivo
DE
40000
30000 0 1 2 3
10 10 10 10
Tiempo CPU (s)
5
10
SSm
Función objetivo
3
10
1
10
0 1 2 3
10 10 10 10
Tiempo CPU (s)
100 2
y1: alfa−pineno y1: alfa−pineno
90 y2: dipenteno y2: dipenteno
y3: allo−ocimeno 1.5 y3: allo−ocimeno
80 y4: pironeno y4: pironeno
Concentración (% peso)
50 0
40
−0.5
30
−1
20
−1.5
10
0 −2
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
Tiempo (min) x 10
4 Tiempo (min) x 10
4
p5 0.8
0.6
p4 0.4
0.2
p3 0
−0.2
p2 −0.4
−0.6
p1 −0.8
−1
p1 p2 p3 p4 p5
Las Figuras 11.8 y 11.9 muestran las lı́neas de contorno de la función objetivo
en el plano paramétrico para el par (p1 , p2 ) que presenta un valor del coeficiente de
correlación R1,2 = 0.13 y para el par más correlacionado (p4 , p5 ) con R4,5 = 0.82.
−5 −5
x 10 x 10
5.5
4
5
3.5
4.5
3
p2
p5
3.5
2.5
3
2
2.5
1.5 2
3 4 5 6 7 8 1.5 2 2.5 3 3.5 4
p −5
x 10 p x 10
−4
1 4
Las Figuras 11.10 y 11.11 muestran las regiones de confianza obtenidas por el
método de Monte Carlo. Como era de esperar, la correlación existente entre el par
(p4 , p5 ) dio lugar a una forma elı́ptica mientras que la correspondiente al par (p1 , p2 )
es prácticamente circular.
5
x 10
5 x 10
4.6
3.02
4.4
3
4.2
2.98
p2
p5
4
2.96
3.8
2.94
3.6
2.92
2.9 3.4
5.86 5.88 5.9 5.92 5.94 5.96 5.98 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9
5 p4
p1 x 10 x 10
4
11.8. Conclusiones
En este capı́tulo se consideró la calibración de un sistema de reacción para la
isomerización térmica del α-pineno. Mediante el método de series de Taylor se de-
mostró la identificabilidad estructural de los cinco parámetros del modelo. A pesar de
que el modelo es aparentemente sencillo, el problema inverso asociado es multimodal
debido, entre otros factores, a dependencias entre los datos. Esto se debe a que las
concentraciones proporcionadas por Fuguitt y Hawkins (1947) para algunos de los
productos no fueron obtenidas experimentalmente sino calculadas matemáticamente
a partir de las concentraciones de otros productos (Box et al., 1973).
Esta multimodalidad fue corroborada mediante un multi-start de un método
SQP motivando el uso de estrategias de optimización global para la estimación de
los parámetros del modelo. La consideración de un rango muy amplio para los lı́mites
de los parámetros y de un punto inicial muy alejado de la solución global hace que,
incluso algunos métodos de optimización global de probada eficacia como SRES y DE,
presenten problemas de convergencia. Sin embargo, la metaheurı́stica SSm probó ser
muy robusta y alcanzó la solución global muy rápidamente lo que la convierte en
una estrategia muy recomendable para la resolución de esta clase de problemas.
El análisis de identificabilidad a posteriori no revela grandes correlaciones en-
tre los parámetros que pudieron ser estimados con precisión como demuestran los
intervalos de confianza calculados.
Capı́tulo 12
12.1. Introducción
Este problema consiste en la estimación de un número de parámetros de un
modelo que describe el mecanismo de reacción para la inhibición irreversible de la
proteasa del HIV originalmente estudiado por Kuzmic (1996) (ver Figura 12.1).
119
120 Capı́tulo 12. Inhibición de la proteasa del HIV
d[M ]
= −2k11 [M ][M ] + 2k12 [E] (12.1)
dt
d[P ]
= k22 [ES] − k21 [P ][E] + k42 [EP ] (12.2)
dt
d[S]
= −k21 [S][E] + k22 [ES] (12.3)
dt
d[I]
= −k21 [I][E] + k52 [EI] (12.4)
dt
d[ES]
= k21 [S][E] − k22 [ES] − k3 [ES] (12.5)
dt
d[EP ]
= k41 [P ][E] − k42 [EP ] (12.6)
dt
d[E]
= k11 [M ][M ] + k12 [E] − k21 [S][E] + k22 [ES] + k3 [ES] (12.7)
dt
−k41 [P ][E] + k42 [EP ] − k51 [I][E] + k52 [EI]
d[EI]
= k51 [I][E] − k52 [EI] − k6 [EI] (12.8)
dt
d[EJ]
= k6 [EI] (12.9)
dt
12.3. Ranking de parámetros 121
500
dmsqr
d
mabs
400 dmean
dmax
dmin
300
Valor del criterio
200
100
−100
k42 k22 k3 k52 k6
Parámetros
Estos resultados reflejan grandes diferencias en los valores de δ msqr para los dis-
tintos parámetros lo que indica que la salida del modelo es muy sensible a unos y
poco sensible a otros. Las diferencias entre δ msqr y δ mabs , indican que existe cierta
variabilidad en las sensibilidades de los distintos estados con respecto a un mismo
parámetro (Sj ). Una comparación de δ max y δ min indica que solamente k6 presentan
sensibilidades tanto positivas como negativas mientras que para los demás paráme-
tros éstas son siempre positivas.
122 Capı́tulo 12. Inhibición de la proteasa del HIV
Para poder comparar los resultados se tomaron los lı́mites y valores iniciales para
los parámetros empleados por Mendes y Kell (1998) (ver Tabla 12.2). Una de las
dificultades añadidas de este problema son los amplios lı́mites considerados para las
constantes cinéticas.
Mediante la minimización de la suma de los cuadrados de los residuos entre
los datos medidos y los simulados, la mejor solución conocida hasta este trabajo
fue obtenida por Mendes y Kell (1998) utilizando el método Simulated Annealing,
con un coste computacional de tres millones de simulaciones. La siguiente mejor
solución fue obtenida utilizando un método Levenberg-Marquardt con un esfuerzo
computacional considerablemente menor (4000 simulaciones) aunque la convergencia
al óptimo global con este método sólo está garantizada si se inicializa en su vecindad.
En este trabajo se trató de resolver el problema mediante un método local tipo
SQP en modo multi-start. El histograma de frecuencias (Figura 12.3) muestra que
este método se quedó atrapado en soluciones locales la mayorı́a de las veces lo que
demuestra que este problema es multimodal por lo que se necesitarán métodos de
optimización global para poder asegurar la convergencia al óptimo global.
35
30
25
Frecuencia
20
15
10
0
0 10 20 30 40 50
Función objetivo
El método SSm convergió a mejores soluciones que las encontradas por Mendes y
Kell (1998) en menos de 1500 simulaciones lo que confirma el buen comportamien-
to de este método incluso en problemas complejos de estimación de parámetros.
Además, cuando se compara con otros métodos estocásticos de demostrada eficacia
124 Capı́tulo 12. Inhibición de la proteasa del HIV
como SRES o DE, SSm alcanzó mejores soluciones con una aceleración en el tiempo
de cálculo de casi tres órdenes de magnitud (ver Figura 12.4).
2
10
SRES
DE
SSm
1
10
Función objetivo
0
10
−1
10
−2
10 0 2 4 6
10 10 10 10
Tiempo CPU (s)
A pesar de que SSm convergió en todas las optimizaciones a valores muy buenos
de la función objetivo (siempre inferiores al mejor valor publicado hasta el momen-
to), los valores de los parámetros no siempre fueron los mismos lo que indica una
función objetivo muy plana en la región del espacio de parámetros cerca del óptimo.
En la Tabla 12.3 se ilustra este hecho mostrando el valor de los parámetros corres-
pondientes a dos soluciones cuyo valor de la función objetivo es muy próximo siendo
los valores de los parámetros significativamente diferentes.
La Figura 12.5 muestra el buen ajuste de los datos experimentales a las predic-
ciones del modelo obtenidas con el mejor vector encontrado por SSm (Solución I). La
Figura 12.6 representa los residuos y en ella se puede apreciar la falta de correlación
de los mismos con respecto al tiempo.
Tabla 12.3: Valor de los parámetros para dos resultados obtenidos con SSm
0.7
Experimento 1 0.03
Experimento 1
Experimento 2
0.6 Experimento 2
Experimento 3
0.02 Experimento 3
Experimento 4
Experimento 4
0.5 Experimento 5
Experimento 5
0.4 0.01
Residuos
Señal
0.3 0
0.2
−0.01
0.1
−0.02
0
−0.1 −0.03
0 500 1000 1500 2000 2500 3000 3500 4000 0 500 1000 1500 2000 2500 3000 3500 4000
Tiempo (s) Tiempo (s)
Figura 12.5: Datos experimentales versus Figura 12.6: Valores de los residuos en
valores predichos por el modelo función del tiempo
126 Capı́tulo 12. Inhibición de la proteasa del HIV
offset(5)
1
E0(5)
S0(5)
0.8
offset(4)
E0(4)
0.6
S0(4)
offset(3)
0.4
E0(3)
S0(3)
0.2
offset(2)
E0(2)
0
S0(2)
offset(1)
−0.2
E0(1)
S0(1)
−0.4
k6
k52
−0.6
k22
k42
−0.8
k3
−1
12.7. Conclusiones
En este capı́tulo se consideró el problema de estimación de parámetros y concen-
traciones iniciales para un modelo de la inhibición irreversible de la proteasa del HIV.
El método SSm demostró ser una estrategia muy eficaz para su resolución alcanzando
valores muy buenos de la función objetivo en un tiempo de cálculo muy razonable y
superando en más de dos órdenes de magnitud los requerimientos computacionales
de otros método de optimización global como SRES, DE o Simulated Annealing.
12.7. Conclusiones 127
13.1. Introducción
La apoptosis consiste en una cascada de reacciones enzimáticas que conduce a la
muerte celular programada o suicidio celular, un proceso que juega un importante
papel desde el desarrollo temprano hasta el envejecimiento. Las caspasas son unas
proteı́nas pertenecientes al grupo de las cisteı́n-proteasas, caracterizadas por presen-
tar un residuo de cisteı́na que media en la ruptura de otras proteı́nas. En el caso de
las caspasas el corte se produce al nivel de un residuo de aspartato de donde deriva
su nombre (cisteinil-aspartato proteasas). Estas enzimas son mediadoras esenciales
de los procesos de muerte celular programada y, una vez activadas, desmantelan la
célula mediante la ruptura selectiva de proteı́nas clave después de residuos de aspar-
tato. Los eventos que culminan en la activación de las caspasas están sujetos a un
intenso estudio debido a su papel en el cáncer y en enfermedades neurodegenerativas
y autoinmunes.
Fussenegger et al. (2000) presentaron un modelo matemático mecanı́stico, descri-
biendo los elementos clave de la activación de las caspasas por medio de receptores
y mecanismos inducidos por estrés (ver Figura 13.1). Este grupo utilizó principios
de conservación de masa junto con leyes sobre las velocidades cinéticas para for-
mular un sistema de ecuaciones diferenciales que describe la evolución temporal de
la activación de las caspasas. El modelo consiste en 19 estados (concentraciones de
proteı́nas) y 11 velocidades de reacción. Se simularon varias estrategias cualitativas
para la prevención de la activación de las caspasas mostrando concordancia con la
información disponible.
Gadkar et al. (2005) consideraron la identificación de este modelo y generaron
129
130 Capı́tulo 13. Función de las caspasas en la apoptosis
FAS/FASL
Activación inducida por el receptor
FADD FADD
Procaspasa-8 Procaspasa-8
FLIP
División proteolítica
Activación inducida por estrés
Caspasa-8 Caspasa-8
Citocroma - c
p53
Citocroma - c
Procaspasa-9 Procaspasa-9
IAPs
kl (x1 − x2 )L
r1 = (13.15)
KS−1 + L
· ¸
x3 x2 x4
r2 = ka − (13.16)
(1 + KA x3 + KA KB x23 ) KA KB x3
" #
x5 x6 x7
r3 = kh x19 − (13.17)
1 + KH x5 + KI 1+KJ x17
KH
k8za1 x28 x4
r4 = (13.18)
KC−1 KD
−1 −1
+ KD x8 + x28 + KF KC−1 KD −1 −1
x15 + KG KD x8 x15
132 Capı́tulo 13. Función de las caspasas en la apoptosis
k9za1 x29 x7
r5 = −1 −1 (13.19)
KK KL + KL−1 x9 + x29 + KN KK KL x16 + KO KL−1 x9 x16
−1 −1
donde
1.5
dmsqr
d
mabs
1 dmean
dmax
dmin
0.5
Valor del criterio
−0.5
−1
−1.5
10 3 16 4 5 19 21 23 25 8 26 27 17 14 13 11 9 7 2 1 15 12 18 20 22 24 6
Parámetros
Estos resultados reflejan grandes diferencias en los valores de δ msqr para los
distintos parámetros lo que indica que la salida del modelo es muy sensible a unos
134 Capı́tulo 13. Función de las caspasas en la apoptosis
y poco sensible a otros. Entre δ msqr y δ mabs no hay grandes diferencias lo que indica
que no existe mucha variabilidad en las sensibilidades de los distintos estados con
respecto a un mismo parámetro (Sj ). Una comparación de δ max y δ min indica que
todos los parámetros presentan sensibilidades tanto positivas como negativas.
18
16
14
12
Frecuencia
10
0
0 5 10 15 20 25
Función objetivo
Por este motivo, el problema de estimación fue también resuelto con los métodos
globales SRES, DE y SSm. El método DE no alcanzó la solución global mientras que
SRES y SSm convergieron a valores muy buenos de la función objetivo aunque el
tiempo de cálculo requerido por SSm fue de menos de 10 segundos, más de un orden
de magnitud inferior al requerido por SRES.
A pesar de que para la estimación solamente se consideraron los datos pseudo-
experimentales correspondientes a la concentración de siete proteı́nas, el ajuste de
todos los estados es bueno como muestra la Figura 13.6. El valor de las velocidades
de reacción predicho por el modelo también se ajusta bien a los valores teóricos a
pesar de que éstos no hayan sido empleados para la estimación (ver Figura 13.5).
136 Capı́tulo 13. Función de las caspasas en la apoptosis
2
10
DE
SRES
SSm
Función objectivo
1
10
0
10
0 1 2 3
10 10 10 10
Tiempo CPU (s)
0.2
0.15
0.01 0.02
0.15
r1
r3
r4
0.1
r
0.1
0.005 0.01
0.05
0.05
0 0 0 0
0 50 100 0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min) Tiempo (min)
−5 −5
0.015 x 10 x 10 0.1
1.4 3.5
3 0.08
1.2
0.01
2.5 0.06
r5
r8
6
r7
1
r
2 0.04
0.005
0.8 0.02
1.5
0 0.6 1 0
0 50 100 0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min) Tiempo (min)
0.04
0.15
0.02
0.03
10
r11
r9
0.1
r
0.02
0.01
0.05
0.01
0 0 0
0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min)
2.5 2.5 1 1
2 0.8
0.8
2
1.5 0.6
x1
x3
x4
0.6
x
1 0.4
1.5
0.4
0.5 0.2
1 0 0.2 0
0 50 100 0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min) Tiempo (min)
0.2
1.5 1.2
2
0.15
x5
x7
x8
1 1
x
0.1
1.5
0.5 0.8
0.05
0 1 0
0 50 100 0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min) Tiempo (min)
1.8 0.4
1 1
1.6 0.3
10
x11
x12
x9
1.4 0.2
0.5 0.5
1.2 0.1
1 0 0 0
0 50 100 0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min) Tiempo (min)
0.9 2
1 2
0.8 1.8
x13
14
x15
x16
x
0.7 1.6
0.5 1.5
0.6 1.4
0 0.5 1
0 50 100 0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min) Tiempo (min)
1.4
0.6 2
1.2
x17
18
x19
0.4 1.5
x
1
0.2 1
0.8
0 0.5
0 50 100 0 50 100 0 50 100
Tiempo (min) Tiempo (min) Tiempo (min)
27
26 0.8
25
21 0.6
19
0.4
16
15
0.2
13
12 0
11
10 −0.2
9
8 −0.4
5
4 −0.6
3
−0.8
2
1
1 2 3 4 5 8 9 10 11 12 13 15 16 19 21 25 26 27
13.7. Conclusiones
En este capı́tulo se consideró la estimación de parámetros en un modelo que
describe los elementos clave de la activación de las caspasas. El empleo de métodos
globales permitió realizar un buen ajuste de los datos experimentales. El método
SSm resultó ser mucho más rápido que otras técnicas globales de probada eficacia.
Sin embargo, el análisis de la identificabilidad a posteriori, demostró la existen-
cia de graves problemas de identificabilidad por lo que no se puede asegurar que los
parámetros estimados puedan reproducir el comportamiento del sistema en condi-
ciones experimentales diferentes.
Capı́tulo 14
14.1. Introducción
141
142 Capı́tulo 14. Ruta bioquı́mica en tres pasos
G1 G2 G3
E1 E2 E3
S M1 M2 P
dG1 V1
= ³ ´ni1 ¡ Ka ¢na1 − k1 G1 (14.1)
dt P
1+ Ki1
+ S
1
dG2 V2
= ³ ´ni2 ³ ´na2 − k2 G2 (14.2)
dt P Ka2
1 + Ki 2
+ M1
dG3 V3
= ³ ´ni3 ³ ´na3 − k3 G3 (14.3)
dt P Ka3
1 + Ki3 + M2
dE1 V4 G1
= − k4 E 1 (14.4)
dt K4 + G1
dE2 V5 G2
= − k5 E 2 (14.5)
dt K5 + G2
dE3 V6 G3
= − k6 E 3 (14.6)
dt K6 + G3
³ ´ ³ ´
1 1
dM1 kcat1 E1 Km1
(S − M1 ) kcat2 E2 Km3
(M1 − M2 )
= S M1
− M1 M2
(14.7)
dt 1+ Km1
+ Km2
1+ Km3
+ Km4
³ ´ ³ ´
1 1
dM2 kcat2 E2 Km3
(M1 − M2 ) kcat3 E3 Km5
(M2 − P )
= M1 M2
− M2 P
(14.8)
dt 1+ Km3
+ Km4
1+ Km5
+ Km6
14.3. Ranking de parámetros 143
Especie Concentración
G1 6.6667e-1
G2 5.7254e-1
G3 4.1758e-1
E1 4.0000e-1
E2 3.6409e-1
E3 2.9457e-1
M1 1.4190e+0
M2 9.3464e-1
1.5
dmsqr
dmabs
dmean
1 dmax
dmin
0.5
Valor del criterio
−0.5
−1
−1.5
5 11 17 4 10 21 16 6 1 12 7 27 24 18 13 19 22 25 28 34 20 26 23 29 35 8 14 2 31 30 32 3 33 9 36 15
Parámetros
Concentración de S Concentración de P
Experimento 1 0.1 0.05
Experimento 2 0.1 0.13572
Experimento 3 0.1 0.36840
Experimento 4 0.1 1.0
Experimento 5 0.46416 0.05
Experimento 6 0.46416 0.13572
Experimento 7 0.46416 0.36840
Experimento 8 0.46416 1.0
Experimento 9 2.1544 0.05
Experimento 10 2.1544 0.13572
Experimento 11 2.1544 0.36840
Experimento 12 2.1544 1.0
Experimento 13 10 0.05
Experimento 14 10 0.13572
Experimento 15 10 0.3684
Experimento 16 10 1.0
dentro de los lı́mites de los parámetros y realizar la búsqueda local desde cada uno
de ellos) dio lugar a un gran número de soluciones locales. En la Figura 14.3, se
puede observar el histograma de frecuencias para el método n2fb en modo multi-
start donde la mayor parte de las soluciones están lejos del óptimo global que no se
alcanzó en ninguna ocasión.
30
25
20
Frecuencia
15
10
0
0 200 400 600 800 1000 1200
Función Objetivo
Esto lleva a confirmar la idea de que sólo los métodos de optimización global son
adecuados para resolver esta clase de problemas.
mayor de la función objetivo) que el método SRES. A pesar de esta ventaja inicial,
puede verse como la fase local del hı́brido proporciona una convergencia mucho más
rápida a una solución mejor, dando lugar a una aceleración total de un orden de
magnitud.
Basándose en la observación de la Figura 14.8, serı́a natural argumentar que un
cambio más temprano de la búsqueda estocástica a la fase local del hı́brido podrı́a
dar lugar a una aceleración todavı́a mayor. Sin embargo, como ya se discutió en
la sección 7.3, si el cambio se realiza demasiado pronto, éste podrı́a dar lugar a la
convergencia a una solución local, es decir, un punto de cambio anterior tiene mayor
probabilidad de estar fuera de la zona de atracción de la solución global. Este efecto
se ilustra en la Figura 14.4 donde se representan tres búsquedas locales realizadas a
partir de diferentes puntos de la misma búsqueda global. Las flechas indican el punto
de cambio a lo largo de la curva de convergencia de SRES (lı́nea contı́nua), mientras
que la convergencia de n2fb para cada punto se representa por lı́neas discontinuas.
Las dos primeras, señaladas con (a), convergieron a soluciones locales, mientras que
la última, señalada con (b), alcanza la solución óptima global.
4
10
2
10
(a)
0
10
Función objetivo
−2
10
−4
10 (b)
−6
10
−8
10
0 2000 4000 6000 8000 10000 12000
Tiempo CPU (s)
que el hı́brido ajustado puede ser aplicado a otros conjuntos de datos, o incluso a
modelos ligeramente diferentes, sin ajuste adicional.
Para el conjunto de datos I, el mejor resultado (J=1.54e-7) se obtuvo después de
un tiempo de cálculo de 3.1 horas. Para el conjunto de datos pseudo-experimentales
II (3 % de error) el mejor resultado (J = 1.25) se obtuvo después de un tiempo de
cálculo de 3.4 horas y para el conjunto de datos pseudo-experimentales III (5 % de
error) se alcanzó un valor de la función objetivo J = 3.27 en un tiempo de cálculo
total de 3.5 horas. En la Tabla 14.5 se muestran los valores del tiempo computacional
y de la función objetivo en cada una de las etapas del método hı́brido (punto inicial,
punto de cambio y resultado final) para los tres conjuntos de datos.
25 25
20 20
15 15
10 10
error relativo (%)
0 0
−5 −5
−10 −10
−15 −15
−20 −20
−25 −25
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Parámetro Parámetro
Figura 14.5: Error relativo considerando Figura 14.6: Error relativo considerando
el conjunto de datos II (3 % de error) el conjunto de datos III (5 % de error)
La Figura 14.4 muestra los valores M1, M2, E1, E2, E3, G1, G2 y G3 teóricos
(lı́nea continua) frente a sus experimentales (marcador) considerando los mejores
parámetros estimados por el método hı́brido secuencial para el conjunto de datos
III en cada uno de los 16 experimentos. Nótese que existe una muy buena correlación
entre los datos experimentales y los predichos incluso para el conjunto de datos con
más ruido. El comportamiento para los otros dos conjuntos de datos es similar por
lo que se omitió su representación.
6
3
5
2.5
4 2
M1
M2
1.5
2 1
1 0.5
0 0
0 20 40 60 80 100 120 0 20 40 60 80 100 120
Tiempo (min) Tiempo (min)
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
E1
E2
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 20 40 60 80 100 120 0 20 40 60 80 100 120
Tiempo (min) Tiempo (min)
0.45 1.4
0.4 1.2
0.35
1
0.3
0.8
0.25
G1
E3
0.2 0.6
0.15
0.4
0.1
0.2
0.05
0 0
0 20 40 60 80 100 120 0 20 40 60 80 100 120
Tiempo (min) Tiempo (min)
1.4 0.8
1.2 0.7
0.6
1
0.5
0.8
G2
G3
0.4
0.6
0.3
0.4
0.2
0.2 0.1
0 0
0 20 40 60 80 100 120 0 20 40 60 80 100 120
Tiempo (min) Tiempo (min)
4
10
SRES
SRES
2 n2fb
10 SSm
Función objetivo
0
10
−2
10
−4
10
Método híbrido
−6
10
2 3 4 5
10 10 10 10
Tiempo CPU (s)
1
35
0.8
30
0.6
25 0.4
0.2
20
0
15 −0.2
−0.4
10
−0.6
5 −0.8
−1
5 10 15 20 25 30 35
1.5 1.5
1.4 1.4
1.3 1.3
1.2 1.2
1.1 1.1
p4
1 1
6
p
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.5 1 1.5 0.5 1 1.5
p p
1 1
puede observar, el valor de los intervalos de confianza es pequeño para todos los
parámetros y los tres conjuntos de datos lo que indica una buena estimación de los
mismos. Como era de esperar, estos intervalos son mayores a medida que aumenta
el error de los datos experimentales.
Aunque en todos los casos los valores de confianza son aceptables, éstos son
mayores para los parámetros que presentan altas correlaciones. Véase, por ejemplo,
que los intervalos para el parámetro p1 que presenta altas correlaciones con otros
parámetros, éstos son un orden de magnitud mayores que para el parámetro p2 que
está poco correlacionado. Esta situación puede ser mejorada mediante un diseño
experimental adecuado.
lı́mites para S y P
también podrı́a plantearse una formulación más (o menos) general del problema.
Por ejemplo, el número de nuevos experimentos N 2exp podrı́a ser considerado como
una variable de decisión, dando lugar a un problema de optimización no lineal en-
tero mixto (Mixed-Integer Non-Linear Programming, MINLP) con un problema de
valor inicial interno. El horizonte de tiempo y los tiempos de muestreo para cada
experimento también podrı́an ser considerados como variables de decisión. Obvia-
mente, aumentar la generalidad de la formulación implica resolver un problema de
optimización más complejo.
Para los objetivos de este trabajo, se consideró que el valor de N 2exp es fijo y
además:
los valores de P y S son invariantes con el tiempo para cada uno de los expe-
rimentos
0.018
fmincon
0.016 NOMADm
DE
0.01
0.008
0.006
0.004
0.002
0 0 1 2 3 4 5
10 10 10 10 10 10
Tiempo CPU (s)
14.8. Conclusiones
En este capı́tulo se consideró la estimación de parámetros y el diseño expe-
rimental óptimo para un modelo de una ruta bioquı́mica. Empleando el hı́brido
SRES+n2fb para el problema de estimación de parámetros se obtuvieron soluciones
mejores que con el método SRES a solas en un tiempo computacional mucho más
reducido. Además, esta metodologı́a demostró ser robusta cuando se maneja ruido
en las medidas. Sin embargo, empleando el método SSm el tiempo computacional
se reduce todavı́a más pasando de un rango de 35-40 horas con SRES a un par de
minutos, conservando la robustez.
El análisis de identificabilidad reveló correlaciones elevadas entre ciertos pares
de parámetros que estaban ocasionando un mal condicionamiento del problema.
Mediante el uso de nuevos diseños experimentales se demostró que esta situación
puede ser mejorada. Los parámetros estimados utilizando datos experimentales ob-
tenidos a partir del diseño óptimo presentaron intervalos de confianza menores que
los estimados con el diseño original.
Los resultados del diseño óptimo de experimentos confirman la utilidad de esta
metodologı́a mejorando la precisión de los parámetros estimados a la vez que se
reduce el esfuerzo experimental.
Capı́tulo 15
15.1. Introducción
La diabetes mellitus se define como un grupo de enfermedades metabólicas que
se caracterizan por altos niveles de glucosa en sangre (hiperglucemia) (Expert Com-
mittee on the Diagnosis and Classification of Diabetes Mellitus, 2003). Esta hiper-
glucemia es el resultado de defectos en la secreción de insulina, en la acción de la
insulina, o en ambas. En la diabetes de tipo 1 hay una deficiencia absoluta de secre-
ción de insulina debido a la destrucción de las células β. Las personas con diabetes
de tipo 1 son propensas a la cetoacidosis y son totalmente dependiente de insulina
exógena. Se estima que en el año 2000 habı́a 17.1 millones de personas con diabetes
de tipo 1 en todo el mundo (Wild et al., 2004; Eiselein et al., 2004).
En diabetes, la hiperglucemia crónica se asocia con complicaciones a largo plazo
debido a daños, disfunciones e insuficiencias en varios órganos, especialmente ojos,
riñones, nervios, corazón y vasos sanguı́neos. Las mayores complicaciones son en-
fermedades cardı́acas, ataques de apoplejı́a, retinopatı́as, nefropatı́as y neuropatı́as.
Éstas pueden eventualmente producir fallos renales, ceguera, amputación y otros
tipos de morbilidad. Los sujetos con diabetes tienen un alto riesgo de sufrir en-
fermedades cardiovasculares y se enfrentan a una mayor morbilidad y mortalidad
cuando son enfermos crı́ticos.
La eficacia de un tratamiento intensivo para la prevención de complicaciones
diabéticas ha sido demostrada por el Ensayo sobre el Control y las Complicacio-
nes de la Diabetes (Diabetes Control and Complications Trial, 1993) y el Estudio
Prospectivo de Diabetes del Reino Unido (United Kingdom Prospective Diabetes
Study, 1998). En ambos ensayos los regı́menes de tratamiento que redujeron la me-
161
162 Capı́tulo 15. Cinética de la glucosa en pacientes diabéticos
dia de la hemoglobina glicosilada A1C (medida clı́nica del control glucémico que
refleja los niveles medios de glucosa en sangre durante los 2-3 meses precedentes) a
un apróximadamente 7 % (el rango normal es de 4-6 %) fueron asociados con pocas
complicaciones microvasculares a largo plazo. A pesar de ello, evidencias recientes
sugieren que estos niveles objetivo no son suficientemente bajos (Khaw et al., 2001;
Muntner et al., 2005).
Los tratamientos intensivos requieren múltiples inyecciones diarias de insulina
(tres o más), o un tratamiento con una bomba de infusión de insulina (ver Figura
15.1). En cualquier caso, este control estricto (lo más cercano posible a la normali-
dad) debe mantenerse de por vida para aprovechar todos los beneficios que confiere.
Hay muchos factores que influyen en la dosis de insulina requerida a través del tiem-
po, incluyendo el peso, la condición fı́sica y los niveles de estrés. Debido a esto, se
requiere una monitorización frecuente de la glucosa en sangre. Basándose en estas
medidas, se puede modificar la dosificación de la insulina, implementar cambios en
la dieta (como alteraciones en los horarios, frecuencia y contenido de las comidas)
y variar las pautas de actividad y ejercicio.
muestran que el modelo es capaz de describir las dinámicas observadas para sujetos
de tipo 1 y por lo tanto puede ser empleado para simular el comportamiento del
paciente en estas condiciones.
dQ1 (t) c
= −F01 − x1 (t)Q1 (t) + k12 Q2 (t) − FR (15.1)
dt
+UG (t) + EGP0 [1 − x3 (t)]
dQ2 (t)
= x1 (t)Q1 (t) − [k12 + x2 (t)] Q2 (t) (15.2)
dt
dx1 (t)
= −ka1 x1 (t) + SIT ka1 I(t) (15.3)
dt
dx2 (t)
= −ka2 x2 (t) + SID ka2 I(t) (15.4)
dt
dx3 (t)
= −ka3 x3 (t) + kb3 I(t) (15.5)
dt
G(t) = Q2 (t)/VG (15.6)
164 Capı́tulo 15. Cinética de la glucosa en pacientes diabéticos
absorción G=Q1/VG
intestinal
UG
EGP0 Q1 k 12 Q2
Fc 01Q1/(GVG)-FR
k a1 k b1
absorción x1
insulina
UI /VI
k a2 kb2
I x2
ke
ka3 k b3
x3
donde
½
c F01 si G ≥ 4.5 mmol/L
F01 = (15.7)
F01 G/4.5 en otro caso
½
0.003(G − 9)VG si G ≥ 9 mmol/L
FR = (15.8)
0 en otro caso
DG AG te−t/tmax,G
UG = (15.9)
t2max,G
canal lento
ku
ka1
Q1a
1 Q2
LDa
ka1
canal rápido
ke
Figura 15.3: Estructura del modelo de infusión de insulina (Wilinska et al., 2005)
dQI1a (t)
= ku(t) − kia1 QI1a (t) − LDa (15.10)
dt
dQI1b (t)
= (1 − k)u(t) − kia2 QI1b (t) − LDb (15.11)
dt
dQI2 (t)
= kia1 QI1a (t) − kia1 QI2 (15.12)
dt
dI(t) 1
= (kia1 QI2 + kia2 QI1b (t)) − ke I(t) (15.13)
dt VI
donde
Vmax,ld QI1a
LDa = km,ld +QI1a
(15.14)
Vmax,ld QI1b
LDb = km,ld +QI1b
(15.15)
1
d
msqr
0.8 dmabs
d
0.6 mean
dmax
0.4 dmin
Valor del criterio 0.2
−0.2
−0.4
−0.6
−0.8
−1
SIT SID F01 tmaxG
Parámetros
80
70
60
50
Frecuencia
40
30
20
10
0
0 0.5 1 1.5 2 2.5 3
Función objetivo x 10
5
Por este motivo, se emplearon los métodos globales SRES, DE y SSm. Los tres
métodos convergieron a la solución óptima global en un tiempo reducido siendo SSm
el más rápido de los tres. En la Figura 15.6 se muestran las curvas de convergen-
168 Capı́tulo 15. Cinética de la glucosa en pacientes diabéticos
4
x 10
6
SRES
DE
5 SSm
4
Función objetivo
0
2 4 6 8 10 20 40 60
Tiempo CPU (s)
En la Tabla 15.3 se muestran los valores de los parámetros obtenidos para cada
uno de los cinco pacientes.
La Figura 15.7 representa el error entre los datos experimentales y los datos
predichos por el modelo para la concentración de glucosa como función del tiem-
po, mostrando buena concordancia para los cinco sujetos. Los valores medios de los
errores absolutos para la concentración de glucosa se dan en la Tabla 15.4. Se ob-
serva que los parámetros estimados son apropiados y que la dinámica del sistema es
capturada. Los errores medios son inferiores al 5 % para todos los sujetos a estudio.
15.4. Estimación de parámetros 169
Paciente 1
10 Paciente 2
Paciente 3
Paciente 4
Paciente 5
Error de Predicción (%)
5
−5
−10
Tabla 15.4: Valor medio de los errores de predicción para los niveles de glucosa
120 140
120
100
100
80
Glucosa (mg/dl)
Glucosa (mg/dl)
80
60
60
40
40
20
datos experimentales 20 datos experimentales
data filtrados datos filtrados
predicción del modelo predicción del modelo
0 0
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350 400
Tiempo (min) Tiempo (min)
Figura 15.8: Ajuste del modelo con los Figura 15.9: Validación del ajuste con
datos del experimento 1 los datos del experimento 2
0.8
tmaxG
0.6
0.4
F01
0.2
−0.2
SID
−0.4
−0.6
SIT
−0.8
−1
SIT SID F01 tmaxG
Cabe destacar que, a pesar de que los parámetros del sistema glucorregulatorio
varı́an entre sujetos y por eso su estimación se realiza por separado, la desviación
estándar de algunos de ellos es relativamente pequeña.
15.7. Conclusiones
En este capı́tulo se consideró la calibración de un modelo para la cinética de
la glucosa en pacientes con diabetes de tipo I. El método SSm demostró converger
en un tiempo computacional reducido al óptimo global proporcionando un buen
ajuste de los datos experimentales. La validación del modelo con otros conjuntos de
datos experimentales fue satisfactoria demostrando que el modelo evaluado es capaz
de describir las dinámicas observadas bajo el protocolo experimental del ensayo
hiperinulinı́mico-euglucémico incluyendo una comida.
Esto sugiere que el modelo considerado puede servir como punto de partida para
la incorporación de otros efectos que ningún otro modelo describe actualmente. Estas
otras dinámicas están relacionadas con la variación circadiana en la sensibilidad
a la insulina, cambios en el ritmo de flujo (debidos y no debidos a la insulina)
dependiendo de los niveles de actividad fı́sica, y respuestas contra-regulatorias a
hipoglucemia, estrés y otros.
Parte IV
Conclusiones
Conclusiones
175
176 Conclusiones
En una segunda parte se realizó un estudio sobre de los diferentes métodos para el
cálculo de sensibilidades y el análisis de identificabilidad obteniéndose las siguientes
conclusiones:
Los resultados obtenidos con los distintos criterios sugieren el uso de otras for-
mulaciones alternativas a ser estudiadas en trabajos futuros como, por ejemplo,
formulaciones multi-objetivo para considerar simultáneamente varios criterios
basados en la matriz de información de Fisher.
Inhibición de la proteasa del HIV: el método SSm resultó ser una estrategia
muy eficaz para la resolución del problema de estimación de parámetros aso-
ciado a este modelo. Sin embargo, el análisis de identificabilidad a posteriori
y los intervalos de confianza calculados demostraron que el buen ajuste de
los datos experimentales no es siempre garantı́a de una buena calibración del
modelo. De hecho, cuando los modelos presentan problemas de identificabi-
lidad, los parámetros estimados para un conjunto de datos no serán capaces
de reproducir el comportamiento del modelo en condiciones experimentales
diferentes.
Apéndices
Apéndice A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% USER MUST INTRODUCE HERE PROBLEM RELATED DATA: %
% --> SIMULATION DATA %
% --> OPTIMIZATION RELATED DATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User may change the name of the function according to his/her necessities
% Please, remember to save the file as name_function.m
function[problem_input,ivp_solver,opt_solver,results]=mendes_input_data;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ANALYSIS TO BE PERFORMED %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ’PE’ Parameter estimation
% ’OED’ Optimal experimental design
% ’prior_analysis’ Performs: Ranking of parameters and local a priori identifiability analysis
% ’post_analysis’ Performs practical identifiability analysis
% ’all’ Performs a priori analysis, OED and a posteriori analysis
problem_input.task=’PE’;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODEL RELATED DATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% results.n_states: number of state variables.
problem_input.n_states=8;
183
184 Apéndice A. Ejemplo de fichero de entrada para el entorno GOSBio
problem_input.ydot=char(...
’ydot(1)=par(28)*y(3)*(1.0d0/par(29))*(u(1)-y(1))/(1.0d0+(u(1)/par(29))+(y(1)/par(30))),...
-par(31)*y(4)*(1.0d0/par(32))*(y(1)-y(2))/(1.0d0+(y(1)/par(32))+(y(2)/par(33)))’,...
’ydot(2)=par(31)*y(4)*(1.0d0/par(32))*(y(1)-y(2))/(1.0d0+(y(1)/par(32))+(y(2)/par(33))),...
-par(34)*y(5)*(1.0d0/par(35))*(y(2)-u(2))/(1.0d0+(y(2)/par(35))+(u(2)/par(36)))’,...
’ydot(3)=par(19)*y(6)/(par(20)+y(6))-par(21)*y(3)’,...
’ydot(4)=par(22)*y(7)/(par(23)+y(7))-par(24)*y(4)’,...
’ydot(5)=par(25)*y(8)/(par(26)+y(8))-par(27)*y(5)’,...
’ydot(6)=par(1)/(1.0d0+(u(2)/par(2))**par(3)+(par(4)/u(1))**par(5))-par(6)*y(6)’,...
’ydot(7)=par(7)/(1.0d0+(u(2)/par(8))**par(9)+(par(10)/y(1))**par(11))-par(12)*y(7)’,...
’ydot(8)=par(13)/(1.0d0+(u(2)/par(14))**par(15)+(par(16)/y(2))**par(17))-par(18)*y(8)’);
% Initial conditions for the parametric sensitivities. This matrix will be zero unless the
% initial condition for a particular state depends on any parameter
problem_input.sens0=zeros(problem_input.n_states,problem_input.n_par);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT DATA TO CALCULATE THE OBJECTIVE FUNCTION. EXPERIMENTAL DATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% problem_input.n_exp: number of experiments
problem_input.n_exp=16;
for iexp=1:problem_input.n_exp
Apéndice A. Ejemplo de fichero de entrada para el entorno GOSBio 185
problem_input.t_in{iexp}=0;
problem_input.t_f{iexp}=120;
% Equidistant measurements
problem_input.t_m{iexp}=[0:6:120];
% measurement_type:
% ’sim’: the code will generate a pseudo_data, using par0, rel_error, and tm
% ’real’: the code will use real data, rel_error must be fixed to 1 in this case
%
problem_input.measurement_type=’sim’;
% For the case experimental data is provided this value should be fixed to 1.0
problem_input.rel_error{iexp}=0.05;
% problem_input.fobj_norm:
% 1: to normalize objective function with max experimental data.
% 0: no normalization is applied problem_input.
problem_input.fobj_norm=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OBJECTIVE FUNCTION FOR OPTIMAL EXPERIMENTAL DESIGN %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% problem_input.fobj_type: Function of the FIM to be minimized
% A_optimality = trace(inv(FIM))
% A_modified = -trace(FIM) (maximize trace(FIM))
% D_optimality = -det(FIM) (maximize det(FIM))
% E_optimality = -min(abs(eig(FIM))) (maximize min(eig(FIM)))
% E_modified = max(abs(eig(FIM)))/min(abs(eig(FIM)))
problem_input.fobj_type=’D_optimality’;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CONTROL RELATED DATA FOR EACH EXPERIMENT {i} %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% problem_input.n_u: number of control variables >=1 for experiment {i}:1-n_exp
problem_input.n_u=2;
problem_input.u{iexp}(1,:)=S_list(iexp);
problem_input.u{iexp}(2,:)=P_list(iexp);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ODE SOLVER RELATED DATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select the solver to be used from the following alternatives
% radau5: BDF based method suitable also for DAEs
% rkf45 : Runge-Kutta-Fehlberg ODE Solver
% ode15s: MATLAB code to be used when fcn.m is provided
% odessa: BDF based method suitable for parametric sensitivity analysis
ivp_solver.name= ’radau5’;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NLP SOLVER RELATED DATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Apéndice A. Ejemplo de fichero de entrada para el entorno GOSBio 187
opt_solver.name=’SSm’;
opt_solver.n_starts=100; % To be used in multistart
opt_solver.local_solver=’n2fb’; % To be used in multistart or only_local
opt_solver.par_min=[1.e-12 1.e-12 0.1 1.e-12 0.1 1.e-12 1.e-12 1.e-12 0.1 1.e-12 0.1 1.e-12 1.e-12,...
1.e-12 0.1 1.e-12 0.1 1.e-12 1.e-12 1.e-12 1.e-12 1.e-12 1.e-12 1.e-12 1.e-12 ,...
1.e-12 1.e-12 1.e-12 1.e-12 1.e-12 1e-12 1.e-12 1.e-12 1.e-12 1.e-12 1.e-12];
opt_solver.par_max=[1.e3 1.e3 10. 1.e3 10. 1.e3 1.e3 1.e3 10. 1.e3 10. 1.e3 1.e3 1.e3,...
10. 1.e3 10. 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3,...
1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3 1.e3];
% opt_solver.y0_guess/y0_in/y0_max:
% initial guess/lower/upper bounds for the initial conditions to be estimated
opt_solver.y0_guess{iexp}=[];
opt_solver.y0_min{iexp}=[];
opt_solver.y0_max{iexp}=[] ;
% opt_solver.t_f{iexp}/tf_min/tf_max
opt_solver.tf_guess{iexp}=[120];
opt_solver.tf_min{iexp}=[120];
opt_solver.tf_max{iexp}=[120];
end %iexp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT RELATED DATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Names for the figures. The name of opt and ivp solvers will be added to the names provided here
results.fit_plot=’fit_plot’;
results.convergence_curve=’conv_curve’;
results.histogram=’histogram’; results.residuals_plot=’residuals’;
results.correlation=’correlation_matrix’;
return
Parte VI
Bibliografı́a
Bibliografı́a
[9] M. Balaban. Effect of volume change in foods on the temperature and mois-
ture content predictions of simultaneous heat and moisture transfer models.
Journal of Food Process Engineering, 12:67–88, 1989.
191
192 Bibliografı́a
[22] J.R. Banga, C.G. Moles y A.A. Alonso. Frontiers In Global Optimiza-
tion., capı́tulo: Global optimization of bioprocesses using stochastic and hybrid
methods., páginas 45–70. Nonconvex Optimization and Its Applications, vol.
74. Kluwer Academic Publishers, 2003.
[39] M.J. Box. Some experiences with a non-linear experimental design criterion.
Technometrics, 12:569–589, 1970.
[48] A.K. Datta. Computer-aided engineering in food process and product design.
Food Technology, 52:44, 1998.
[51] J.E. Dennis, D.M. Gay y R.E. Welsch. Algorithm 573, NL2SOL - an
adaptive nonlinear least-squares algorithm. ACM Trans. Math. Software,
7:369–383, 1981.
[60] A.M. Dunker. The decoupled direct method for calculating sensitivity coef-
ficients in chemical kinetics. Journal of Chemical Physics, 81:2385–2393, Sep-
tember 1984.
[61] T. D. Durance. Improving canned food quality with variable retort tempe-
rature processes. Trends in Food Sci. & Technol., 8(4):113–118, 1997.
[63] B Efron. The jackknife, the bootstrap, and other resampling plans. Series
in Applied Mathematics, SIAM, 38:91, 1982.
[67] R. Fisher. Design of experiments. Oliver and Boyd Ldt., Edinburgh, 1935.
[69] C.A. Floudas. Deterministic Global Optimization: Theory, Methods and Ap-
plications. Kluwer Academics, The Netherlands, 2000.
[70] L.J. Fogel, A.J. Owens y M.J. Walsh. Artificial intelligence through
simulated evolution. Wiley, New York, 1966.
[73] K.G. Gadkar, R. Gunawan y F.J. Doyle III. Iterative approach to model
identification of biological networks. BMC Bioinformatics, 6:155, 2005.
[103] P. Kuzmic. Program dynafit for the analysis of enzyme kinetic data: appli-
cation to hiv proteinase. Analytical Biochemistry, 237:260–273, 1996.
[104] P.J.M. Laarhoven y E.H.L. Aarts. Simulated annealing: theory and ap-
plications. Reidel Publishing Company, Dordrecht, 1987.
[109] L. Ljung. System Identification: Theory for the User. Prentice Hall, Engle-
wood Cliffs, New Jersey, 1999.
[118] Z. Michalewicz y D.B. Fogel. How to Solve it: Modern Heuristics. Springer-
Verlag, Berlin, 2000.
Bibliografı́a 201
[125] H.B. Nahor, N. Scheerlinck, J.F. Van Impe y B.M. Nicolaı̈. Opti-
mization of the temperature sensor position in a hot wire probe set up for
estimation of the thermal properties of foods using optimal experimental de-
sign. J. Food Eng., 57(1):103–110, 2003.
[133] K. J. Park. Diffusional model with and without shrinkage during salted fish
muscle drying. Drying Technology, 16(3-5):889–905, 1998.
[134] R.S. Parker, F.J. Doyle III y N.A. Peppas. The intravenous route to
blood glucose control. IEEE Engineering in Medicine and Biology Magazine,
20(1):65–73, 2001.
[135] R.K. Pearson. Outliers in process modeling and identification. IEEE Tran-
sactions on Control Systems Technology, 10(1):55–63, 2002.
[144] Ph. Preux y E.-G. Talbi. Towards hybrid evolutionary algorithms. Intl.
Trans. in Op. Res., 6:557–570, 1999.
[146] H.S. Ramaswamy, G.B. Awuah y B.K. Simpson. Heat transfer and
lethality considerations in aseptic processing of liquid/particle mixtures: A
review. Crit. Rev. in Food Sci. Nutr., 37(3):253–286, 1997.
[157] T.P. Runarsson y Xin Yao. Stochastic ranking for constrained evolutionary
optimization. IEEE Transactions on Evolutionary Computation, 4:284–294,
2000.
[159] A. Saltelli, K. Chan y E.M. Scott. Sensitivity Analysis. John Wiley &
Sons, Inc., New York, 2000.
[160] W.E. Schiesser. The Numerical Method of Lines. Academic Press, New
York, 1991.
[162] H. P. Schwefel. Evolution and optimum seeking. Wiley, New York, 1995.
[164] G.A.F. Seber y C.J. Wild. Non Linear Regression. Wiley, New York, 1989.
[173] J. Timmer. Modeling noisy time series: physiological tremor. Int. J. Bifur-
cation Chaos, 8(7):1505–1516, 1998.
[185] K. Versyck. Dynamic Input Design for Optimal Estimation of Kinetic Para-
meters in Bioprocess Models. Tesis doctoral, Katholieke Universiteit Leuven,
Mayo 2000.
[190] F.S. Wang y J.P. Chiou. Optimal control and optimal time location pro-
blems of differential-algebraic systems by differential evolution. Ind. Eng.
Chem. Res., 36(12):5348–5357, 1997.
[192] M.E. Wilinska, L.J. Chassin, H.C. Schaller, L. Schaupp, T.R. Pieber
y R. Hovorka. Insulin kinetics in type–1 diabetes: continuous and bolus de-
livery of rapid acting insulin. IEEE Transactions on Biomedical Engineering,
52(1):3–12, 2005.
[193] D.H. Wolpert y W.G. Macready. No free lunch theorems for optimization.
IEEE Transactions on Evolutionary Computation, 1:67–82, 1997.
[194] Y. Ye. Interior algorithms for linear, quadratic and linearly constrained non-
linear programming. Tesis doctoral, Department of ESS, Stanford University,
1987.
Publicaciones
Publicaciones 211
Contribuciones a congresos
Balsa-Canto, E., Rodriguez-Fernandez, M., Alonso, A. A. y Banga,
J. R. Optimal Identification in Systems Biology. Applications in Cell Signa-
ling. AIChE Annual Meeting, 12-17 Noviembre, 2006, San Francisco, Estados
Unidos.
Banga, J.R., Balsa-Canto, E., Moles, C. G., Garcia, S., Sendin, O. H.,
Rodriguez-Fernandez, M. y Alonso, A. A. Advances in the Optimization
Publicaciones 213