Tesis Doctorado Quintero Final
Tesis Doctorado Quintero Final
Tesis Doctorado Quintero Final
DOCTOR EN INGENIERÍA
FILTROS DE PARTÍCULAS EN SISTEMAS
DINÁMICOS NO LINEALES NO
GAUSSIANOS:
IDENTIFICACIÓN, ESTIMACIÓN DE
ESTADOS Y APLICACIONES
Aprobada por la
Facultad de Ingeniería de la Universidad Nacional de San Juan
Para el otorgamiento del grado académico del grado académico de
DOCTOR EN INGENIERÍA
ISBN: 978-987-05-8642-5
In this thesis, the approach to be followed to solve the problem of the lack
information of some states in a dynamic system is the called “no lineal no
Gaussian filtering”. It means the well known Particle Filtering. This
technique is based on the Bayes Theory and the sequential methods of Monte
Carlo. The basis of this methodology is the approximation of the distributions
with randomly measurements built by particles, where the recursive
calculation of probability distribution functions (pdf) of the states, the
importance sampling and also the approximation of discrete randomly
measurements are involved. That is, the stochastic differential equations that
characterize the system in continuous time are discretized to obtain stochastic
equations in differences. At this point it is possible to propose a valid
probabilistic formulation using the Bayes Theory, where an optimal solution
is obtained from the probability density function of the states. But if the pdf
is not exactly known, must be used a recursive filter where the received data
must be processed sequentially. The continuous pdf of the states are
approached with discrete random measurements composed by weighted
particles and its weight. That is, the use of an independent number of
particles, sampled directly from the state spaces with the objective to
represent the posteriori probability and updating the posteriori of new
observations is posed. This methodology has the main advantage of not to
use linearization of the system to be approximated. Besides, has a big
potential to its parallel implementation over systems that require bigger
complexity and computational load. Algorithms incorporation into the
automatica main objective (the design of control systems that include the
algorithm in estimation schemes) depends on their success.
This work makes possible to apply this methodology to the automatic control
of processes and bio processes such as the continuous fermentation of the
bacteria Zymomonas mobilis (Z.m) for fuel production and the batch
fermentation of bacillus turingiensis (bt) to the sporulated cells production.
After the building of the state estimators for the two previously mentioned
bio process, where applied techniques of numerical methods based control
developed previous by Scaglia (Scaglia, 2006) and applied in trajectory
tracking of mobile robotics (Scaglia et al, 2007, 2008), over the Z.m
continuous alcoholic fermentation process.
Bibliografía………………………………………………………………..249
Capítulo 3
Tabla 3.1: Tabla resumen de los métodos de Monte Carlo y los métodos de
Filtrado no Lineal 43
Capítulo 5
Tabla 5.1: Parámetros Comportamiento Oscilatorio sostenido 140
Tabla 5.2: Parámetros del Bioreactor Parámetros usados para simulación del
bioreactor 140
Tabla 5.3. Tiempo de cómputo en minutos para cada prueba del proceso de
estimación experimento 1. 145
Tabla 5.4: Errores obtenidos en diferentes pruebas de simulación, ilustración
de los peores casos 148
Tabla 5.5: Parámetros usados para estimación con datos reales de
Zymomonas mobilis 155
Tabla 5.6 Tiempo de estimación de Biomasa en simulación (minutos) para
una simulación completa de 150 horas 156
Tabla 5.7. Tiempo de cómputo en minutos para cada ejecución del proceso
de estimación experiencia 5 160
Tabla 5.8. Tiempo de cómputo en minutos para cada ejecución del proceso
de estimación prueba 7 161
Tabla 5.9. Tiempo de cómputo en minutos para cada ejecución del proceso
de estimación experiencia 9 164
Tabla 5.10. Tiempo de cómputo en minutos para cada ejecución del proceso
de estimación experiencia 11 167
Tabla 5.11. Dinámica del Oxígeno Disuelto 177
Tabla 5.12: Resumen del error máximo de estimación de estados en las
pruebas iniciales para Bacillus Turingiensis 182
Capítulo 6
Tabla 6.1 Resumen de datos de las simulaciones para prueba de
controladores 217
Lista de Figuras
Capítulo 1
Figura 1.1 Esquema de metodología de filtrado 3
Capítulo 2
Figura 2.1 Ejemplo gráfico comparativo de esquemas de muestreo 41
Capítulo 3
Figura 3.1 Comparación del muestreo por Importancia y muestreo por
rechazo 49
Figura 3.2 Descripción esquemática del paso de re muestreo 60
Figura 3.3 Descripción gráfica del filtro de Partículas con re muestreo 61
Figura 3.4 Diagrama de Bloques Algoritmo de Filtro de Partículas con re
muestreo 62
Figura 3.5 Ilustración del paso de re muestreo en un filtro de partículas 67
Capítulo 4
Figura 4.1 Modelo de Incertidumbre de transformación lineal fraccional 118
Capítulo 5
Figura 5.1 Diagrama de proceso de Fermentación alcohólica en continuo 135
Figura 5.2 Esquema de Estimación en simulación para el proceso de
fermentación de Z.m 141
Figura 5.3 a. Estimación de Biomasa con 50 partículas y esquema de re
muestreo Multinomial 144
b. Errores en la estimación de los estados del bioreactor, durante una prueba
en simulación de 150 horas con un tiempo de muestreo de 0.1 horas 145
Figura 5.4 a. Concentración de biomasa estimada por el filtro SIR. La línea
punteada describe la estimación mientras que la línea continua describe el
valor del modelo considerado como real 147
b. Concentración de biomasa estimada por el filtro SIR. La línea punteada
describe la estimación mientras que la línea continua describe el valor del
modelo considerado como real 147
Figura 5.5 a. Dinámica de inhibición Z estimada por el filtro SIR. Línea
punteada representa las estimaciones y las líneas continuas son consideradas
las dinámicas reales del sistema Dinámicas de Inhibición estimadas por el
filtro SIR 149
b. Dinámica de inhibición I estimada por el filtro SIR. Línea punteada
representa las estimaciones y las líneas continuas son consideradas las
dinámicas reales del sistema 149
Figura 5.6 Concentración de biomasa estimada 150
Figura 5.7 a. Valor de Errores de estimación 150
b. Concentración de biomasa estimada 151
c. Estimación de Estados de Inhibición Z 151
d. Estimación de Estados de Inhibición I 151
Figura 5.8 Esquema de Estimación usado para la prueba con datos reales 153
Figura 5.9 Interpolación de datos reales.
a. Concentración de Sustrato 154
b. Concentración de Producto 154
c. 155
Figura 5.10 Concentración de biomasa estimada
a. Estimación de Biomasa experimento 5 usando 10 partículas y re muestreo
residual 158
b. Estimación de Biomasa experimento 5 usando 10 partículas y re muestreo
multinomial 159
c. Estimación de Biomasa experimento 5 usando 10 partículas y re muestreo
determinístico 159
d. Estimación de Biomasa experimento 5 usando 100 partículas y re muestreo
multinomial 159
Figura 5.11 a. Estimación de Biomasa experimento 7 usando 10 partículas y
emuestreo residual 162
b. Estimación de Biomasa experimento 7 usando 10 partículas y re muestreo
multinomial 162
c. Estimación de Biomasa experimento 7 usando 10 partículas y re muestreo
determinístico 162
d. Estimación de Biomasa experimento 7 usando 100 partículas y re muestreo
determinístico 163
Figura 5.12 a. Estimación de Biomasa experimento 9 usando 10 partículas y
Re muestreo multinomial 165
b. Estimación de Biomasa experimento 9 usando 50 partículas y re muestreo
determinístico 165
c. Estimación de Biomasa experimento 9 usando 100 partículas y re muestreo
multinomial 165
Figura 5.13 a. Estimación de Biomasa experimento 11 usando 10 partículas y
Re muestreo residual 167
b. Estimación de Biomasa experimento 11 usando 10 partículas y re muestreo
determinístico 168
c. Estimación de Biomasa experimento 11 usando 10 partículas y re muestreo
multinomial 168
d. Estimación de Biomasa experimento 11 usando 50 partículas y re muestreo
determinístico 168
Figura 5.14 Estimación de Biomasa, las líneas continuas son las de la
interpolación de los datos reales, las líneas punteadas corresponden a la
estimación del filtro 168
Figura 5.15 Resultados de la prueba 4 en la cual se usó un ruido de baja
potencia como perturbación de las medidas 170
a. Concentración de Biomasa estimada en la prueba 4 con baja potencia en la
componente de ruido 171
b. Estimación de la variable de Inhibición Z 171
c. Estimación de la variable de Inhibición I 172
Figura 5.16 a. Estimación de producto, en línea continua los datos de la
interpolación de los datos reales, en punteado las estimaciones del filtro 173
b. Estimación de producto, en línea continua los datos de la interpolación de
los datos reales, en punteado las estimaciones del filtro 173
Figura 5.17 Diagrama de Bloques del Estimador de Biomasa 181
Figura 5.18 Estimación de Biomasa con modelos en simulación y
componentes de difusión adicionadas. Tiempo de muestreo (t0) de 1 hora, re
muestreo residual. 183
Figura 5.19 Estimación de Biomasa con modelos en simulación y
componentes de difusión adicionadas. Tiempo de muestreo (t0) de 0.01 horas
re muestreo determinístico 183
Figura 5.20 Resultados normalizados de estimación de células vegetativas
para la prueba 7 del filtro SIR con esquema de re muestreo residual. Línea
continua el modelo considerado como real, puntos el resultado del filtro 184
Figura 5.21 Resultados normalizados de estimación de células esporuladas
para la prueba 7 del filtro SIR con esquema de re muestreo residual 185
Figura 5.22 Resultados normalizados de sustrato primario para la prueba 7
del filtro SIR con esquema de re muestreo residual 185
Figura 5.23 Resultados normalizados de sustrato secundario para la prueba 7
del filtro SIR con esquema de re muestreo residual 186
Figura 5.24 Resultados normalizados de oxígeno disuelto para la prueba 7 del
filtro SIR con esquema de re muestreo residual 186
Figura 5.25 Superficie generada por el máximo de los errores de estimación
de células vegetativas, usando los 3 tipos de re muestreo a 0.5 horas de
muestreo 187
Figura 5.26 Superficie generada por el máximo de los errores de estimación
de células vegetativas 187
Figura 5.27 Superficie generada por el máximo de los errores de estimación
de células vegetativas 188
Figura 5.28 Superficie generada por el máximo de los errores oxigeno
disuelto 188
Figura 5.29 Superficie generada por el máximo de los errores oxigeno
disuelto usando los 3 esquemas de re muestreo a 0.5 horas de muestreo 189
Figura 5.30 Resultados normalizados de Células Vegetativas para la prueba
del filtro SIR usando muestreo residual 189
Figura 5.31 Resultados normalizados de Células Esporuladas para la prueba
del filtro SIR usando muestreo residual 190
Figura 5.32 Sustrato primario para la prueba del filtro SIR usando muestreo
residual, a 0.2 horas de tiempo de muestreo 190
Figura 5.33 Resultados normalizados de Sustrato secundario para la prueba
del filtro SIR usando muestreo residual 191
Figura 5.34 Resultados normalizados de Oxígeno disuelto para la prueba del
filtro SIR usando muestreo residual 191
Figura 5.35 Células vegetativas para la prueba de el filtro SIR usando
muestreo residual 192
Figura 5.36 Sustrato primario para la prueba de el filtro SIR usando muestreo
residual 193
Figura 5.37 Oxígeno disuelto para la prueba de el filtro SIR usando muestreo
residual, a 0.01 horas de tiempo de muestreo 193
Figura 5.38 Diagrama de Bloques del Estimador de Biomasa 194
Figura 5.39 a. Modelado de Incertidumbre de la concentración de Biomasa
para la fermentación A 196
b. Modelado de Incertidumbre para la concentración de sustrato para la
fermentación A 197
c. Modelado de Incertidumbre para la concentración de oxígeno disuelto de la
fermentación A 197
Figura 5.40 a. Modelado de Incertidumbre de la concentración de Biomasa
para la fermentación B 198
b. Modelado de Incertidumbre para la concentración de sustrato para la
fermentación B 198
c. Modelado de Incertidumbre para la concentración de oxígeno disuelto de la
fermentación B 199
Figura 5.41 a. Modelado de Incertidumbre de la concentración de Biomasa
para la fermentación C 200
b. Modelado de Incertidumbre para la concentración de sustrato para la
fermentación C 200
c. Modelado de Incertidumbre para la concentración de oxígeno disuelto de
la fermentación C 201
Figura 5.42 Entradas para el estimador de Biomasa propuesto
a. Concentración de Sustrato primario (Sp) 202
b. Concentración de Oxígeno Disuelto (OD) 202
Figura 5.43 Estima de Biomasa para la fermentación B usando como modelo
de incertidumbre el construido a partir de la fermentación B 202
Figura 5.44 Estima de Biomasa para la fermentación B usando como modelo
de incertidumbre el construido a partir de la fermentación C 203
Figura 5.45 Estima de Biomasa para la fermentación C usando como modelo
de incertidumbre el construido a partir de la fermentación B 204
Figura 5.46 Estima de Biomasa para la fermentación C usando como modelo
de incertidumbre el construido a partir de la fermentación C 204
Capítulo 6
Figura 6.1 Diagrama del proceso continuo de fermentación 205
Figura 6.2 Desempeño del controlador de Sustrato y producto, para
fermentación 1 218
Figura 6.3 Acciones de control Ds y Sin para fermentación 1 219
Figura 6.4 Desempeño del controlador de Sustrato y producto, para
fermentación 2 220
Figura 6.5 a. Acción de control de sustrato de alimento para la fermentación
2 221
b. Acción de control de tasa de dilución para la fermentación 2 222
Figura 6.6 Desempeño del controlador de Sustrato y producto, para
fermentación 3 222
Figura 6.7 Acciones de control Ds y Sin para fermentación 3 223
Figura 6.8 Desempeño del controlador 2 con recirculación fija del 10% ante
perturbaciones Fermentación 1 226
a. Concentración de Producto en la fermentación controlada 226
b. Concentración de sustrato en la fermentación controlada 226
c. Acción de control concentración de sustrato de alimento 227
d. Acción de control de tasa de dilución 227
Figura 6.9 Desempeño del controlador 2 siguiendo trayectorias obtenidas de
datos experimentales Fermentación 2
a. Concentración de sustrato en la fermentación controlada 228
b. Concentración de Producto en la fermentación controlada 228
c. Acción de control de tasa de dilución 229
d. Acción de control concentración de sustrato de alimento 229
Figura 6.10 Desempeño del controlador 3 siguiendo trayectorias obtenidas de
datos experimentales 235
Figura 6.11 Acciones de Control para el seguimiento de trayectoria 236
Figura 6.12 Desempeño del controlador 3 ante perturbaciones siguiendo
trayectorias obtenidas de datos experimentales
a. Variables de control Biomasa, Susutrato y Producto 237
b. Variables manipuladas recirculación de microorganismos, sustrato de
entrada, tasa de dilución de sustrato y tasa de dilución total 237
Figura 6.13 Esquema general de estimación y control 239
Figura 6.14 Controlador de Sustrato y producto en lazo cerrado con
estimador de estados 240
Figura 6.15 Controlador de Sustrato y Producto en lazo cerrado con
estimador de Biomasa siguiendo trayectorias basadas en datos reales 241
Figura 6.16 Controlador mejorado y estimador de biomasa para control de la
población de Células de Z.m 242
Figura 6.17 Controlador mejorado y estimador de biomasa para control de
concentración de sustrato de alimento 243
Figura 6.18 Controlador mejorado y estimador de biomasa para control de
concentración de producto 243
1
Capítulo 1
INTRODUCCIÓN Y MOTIVACIÓN
Los procesos reales poseen, por naturaleza, dinámicas altamente no lineales debido a
las interacciones entre muchas de las variables internas del sistema. Estas dinámicas,
en algunos casos, pueden describirse mediante la representación del sistema en
ecuaciones diferenciales y, posteriormente, ser expresadas en términos de las variables
de estado, a través de las cuales, se tratan de establecer las relaciones entre las
entradas, los estados y las salidas del sistema. La forma en que la información se
propaga a través de la evolución del sistema determina los valores de los estados y su
repercusión en la salida global del sistema no lineal.
Los denominados Filtros de Partículas son una clase de filtros Bayesianos recursivos,
diseñados para la estimación de los estados de sistemas dinámicos no lineales con
perturbaciones no Gaussianas.
Para aplicar esta metodología, se plantea una formulación del problema en tiempo
discreto mediante ecuaciones en diferencias que describen los diferentes tipos de
interacciones. Éstas generalmente representan el modelo del sistema que describe la
evolución del estado en el tiempo y el modelo de incertidumbre, el cual relaciona las
medidas inciertas o ruidosas con el estado.
Ahora bien, es posible plantear una formulación válida en forma probabilística bajo la
aproximación Bayesiana, donde una solución óptima se obtiene a partir de la función
de densidad de probabilidad del estado (pdf). Pero si la pdf no se conoce exactamente,
como se verá más adelante, se debe usar un filtro recursivo donde los datos recibidos
pueden ser procesados secuencialmente. En otros términos, se plantea el uso de un
número independiente de partículas, muestreadas directamente del espacio de estado
con el fin de representar la probabilidad a posteriori, y actualizando a cada instante la
distribución con las nuevas observaciones. A partir del conocimiento de la
probabilidad a posteriori, es posible obtener información de la mejor estimación del
estado desconocido.
1
Variables aleatorias que son muestras desconocidas del espacio de estado
2
Masas de probabilidad calculadas con la teoría de Bayes
3
De esta forma, la estimación de los estados puede verse como un problema de filtrado
óptimo. Bajo la perspectiva Bayesiana, si las ecuaciones de estado son lineales y la
función de densidad de probabilidad de los estados a posteriori es gaussiana, el filtro
de Kalman proporciona esta solución óptima. Sin embargo, cuando no se verifican las
hipótesis acerca de las características probabilísticas de los estados, no existe una
solución analítica y, por consiguiente, deben realizarse aproximaciones subóptimas.
Un ejemplo de estas aproximaciones es el filtro de Kalman extendido, que asume que
la densidad a posteriori es gaussiana, y adopta una expansión en series de Taylor para
hacer una aproximación local del vector de estados actual. En la práctica, cuando las
ecuaciones de estado son altamente no lineales y la densidad a posteriori no es
gaussiana, el filtro de Kalman extendido puede presentar grandes errores de
estimación (Quintero et al, 2005). Alternativamente, uno de los algoritmos propuestos
en la literatura es el filtro aproximado basado en grillas, introducido en la ingeniería
de control aplicada a procesos por Terwiesch (1995). El costo computacional de este
filtro basado en grillas incrementa exponencialmente con la dimensión del espacio de
estados limitando la amplitud de su aplicación (Chen et al., 2003).
dxt
f xt Ruido1 (t ) (1.1)
dt
dyt
g xt Ruido2 (t ) (1.2)
dt
El estado xt y la observación yt están correlacionadas. Es posible que con las
componentes de ruido agregadas existan problemas en el planteamiento del problema
en tiempo discreto, considerando iteraciones ruidosas de algunas funciones f y g; sin
embargo con la buena elección de las características aleatorias el problema queda bien
caracterizado.
dxt f xt , t dt 1 xt , t dwt
(1.3)
dyt g xt , t dt 2 xt , t dvt
3
También se los denomina en la bibliografía Movimiento Browniano o procesos de Wiener Levy
5
t
t yt E g s Fsy ds (1.4)
0
t
E g ( xs ) 2 ds (1.5)
0
Las filtraciones producidas por la innovación y las observaciones son las mismas y
contienen la misma información. Sin embargo, (Fujisaki et al, 1972) demostró que
independientemente de los resultados que muestran la equivalencia de las filtraciones
de la innovación y las observaciones, la solución de filtrado puede expresarse en
términos de una integral estocástica respecto al proceso de innovación. Esta
aproximación se basa en el teorema de Girsanov (Girsanov, 1960) para transformar la
medida de probabilidad y así obtener una martingala respecto a Fsv . Puede usarse el
hecho que cualquier martingala continua de cuadrado integrable (respecto a Fsv ) y
puede representarse como una integral estocástica respecto al movimiento Browniano.
t t
t ( ) 0 ( ) t ( L )ds s (g ) s ( h) s ( ) dys s ( g) (1.6)
0 0
Donde la ecuación anterior permite una evaluación recursiva, pero tiene la dificultad
de que la media condicional depende de todos los momentos condicionales de orden
superior, lo que lleva a un sistema de dimensión infinita. Dado que esas ecuaciones no
poseen una solución cerrada, muchas investigaciones han tratado de determinar las
condiciones para la existencia de filtros no lineales, usando la teoría de Lie.
Afortunadamente, la ecuación puede ser mapeada a una ecuación diferencial parcial
estocástica (stochastic partial differential equation, SPDE), manteniendo la densidad
condicional; esta es la llamada ecuación que Zakai propuso en 1969 (Zakai, 1969).
1.1 Objetivos
Se investigó el estado del arte del filtrado no lineal y estudiaron las técnicas
más importantes desde la perspectiva Bayesiana que los fundamenta,
consiguiendo de este modo analizar las medidas comparativas para el
mejoramiento de los algoritmos y su implementación práctica.
sobre varios tipos de sistemas dinámicos, es posible plantear una metodología para
definir dichas funciones de distribución que, como se ha mencionado anteriormente,
representan la base de los algoritmos secuenciales de filtrado.
Capítulo 2
PROBLEMA DE FILTRADO EN GENERAL
Aunque en dicho algoritmo se asume que las distribuciones de los términos del error
en el modelo de espacio de estado son normales y/o las ecuaciones de las medidas y
las transiciones son lineales, no es posible obtener una ecuación explícita para el
algoritmo de filtrado. Dado que estas condiciones no siempre se cumplen, las
estimaciones de filtrado obtenidas resultan no ser reales.
Como posible solución se planteó una mejora (filtro de Kalman extendido) a través de
una aproximación de las ecuaciones no lineales de las medidas y de las transiciones,
usando las expansiones por series de Taylor y se aplicaron las ecuaciones linealizadas
directamente al algoritmo lineal recursivo del filtro de Kalman. Lo anterior, pudo
hacerse porque el hecho de aplicar las funciones linealizadas al filtro de Kalman
convencional implica que los términos no gaussianos del error con media no cero, se
aproximen como los términos del error con media cero.
Algunos estudios relacionados con el tema son los de Kitagawa e Higuchi (Kitagawa e
Higuchi, 1987) quien propuso una aproximación, representando las funciones de
densidad numéricamente por una función lineal por tramos, donde cada función de
densidad está representada por un número de segmentos, ubicación de nodos, el valor
de cada nodo, y su evolución a través de integración numérica. De acuerdo a esta
aproximación proporcionada por la integración numérica, el costo computacional
incrementa más que proporcionalmente al aumento de la dimensión del vector de
estados. Su programación es compleja si la dimensión es grande. Este tipo de filtrado
numérico de sistemas no lineales no normales se llama Filtro de integración numérica
(en inglés Numerical Integration Filter NIF).
Uno de los problemas en el filtro de muestreo por importancia, consiste en que, para la
aproximación de las funciones de densidades de filtrado (funciones de densidad
12
Como los métodos recursivos son el principal contenido de esta tesis, se tratarán a
profundidad posterior al estudio de las teorías básicas, en el capítulo 3. Ahora bien, en
este punto es necesario definir claramente las operaciones que se consideran dentro de
la teoría de filtrado en general y que definen y dan forma a los problemas que se van a
plantear más adelante.
13
Como manera de describir las cualidades de una o más variables, se les llaman
“estadísticas bayesianas suficientes” a las expresiones del tipo p(x/y), donde y son las
medidas u observaciones del sistema, y la variable x es la información que no se
dispone explícitamente. Dicha estadística cumple con la condición de que las
probabilidades condicionales sean de la forma p(x/y)= p(x/y’), para todo par de
observaciones y e y’. Para su uso, debe comenzarse con la selección de un modelo
donde se incluyan los datos tomados y los cálculos de sus probabilidades previas, y
posteriormente estimar los parámetros para incorporar los datos en el modelo a
conformar.
La diferencia principal entre ambos puntos de vista para la estimación de los estados
en sistemas dinámicos, radica en la aproximación Bayesiana al problema conocido
como método del máximo a posteriori (MAP), los estados y las medidas de las salidas
del sistema son considerados como variables aleatorias. Mientras que la teoría
desarrollada por Fisher llamada también como el método de máxima verosimilitud (en
14
inglés maximum likelihood ML), el estado está planteado como fijo pero desconocido
y no es una variable aleatoria.
Desde el punto de vista bayesiano que nos ocupa, el objetivo planteado es usar la
información obtenida a partir de las observaciones de las variables aleatorias que
representan las medidas de las salidas y(t) y las entradas u(t) con el fin de inferir
información de las otras variables aleatorias representadas por los estados del sistema
x(t).
4
Ver Apéndice A Estadísticas suficientes desde el punto de vista Fisheriano.
5
Ver Apéndice A definición de realización en un proceso estocástico.
15
Se asume que la evolución de cada uno de los estados del sistema obedece a la
siguiente ecuación diferencial estocástica:
dxt f xt , t dt 1 xt , t dwt
(2.1)
Teorema de Bayes
p( y ) p( y / x) p( x)dx (2.4)
n
Ecuación de Kushner
“El saber si es más útil conocer el modo condicional que la media condicional para
algún problema de filtrado, es un punto de controversia, dependiendo del problema
específico”.
xt f x t (2.6)
t g xt t (2.7)
da PPaa1 g a' 1
dy gdt Paa1 L P a
dt
1 2 1 (2.8)
P 2 Paa1 dy ' 1
g Paa1 g a' 1
dy P Paa u
aa 2
1
ui dy ' g a Paa1 Pai aa Paa1 g a' 1
dy (2.9)
1 (2.10)
L P fi P v P
i ai 2 i, j ij a a
i j
20
Q
Note que L Q es la ecuación de Kolmogorov6 para la función de densidad de
t
un proceso no condicionado.
Ecuación de Zakai
6
Ecuación de Kolmogorov: Según el teorema 4.3.1 en (Kushner, 1967) siendo aij(.) y fi(.) acotados y que
satisfagan la condición uniforme de Holder, y sea a(.) uniforme definida positva. Entonces existe un proceso
de difusión homogéneo con un generador diferencial L. Su densidad de transición p(x,t,y) (la densidad de la
función de transición de Markov P(x,t,y)) es la solución fundamenal de la ecuación diferencial (en variables
x,t donde x es igual al estado inicial) u / t x, t Lu x, t con condiciones
iniciales u x, t x y ,t 0 . El proceso es fuertemente Feller y Markoviano, y es una solución única
t t
de Y (t ) x f Y ( s), s ds Y (s), s dw(s)
0 0
21
Considere el problema de estimar una difusión (ver apéndice A) In valuada x(r), que
evoluciona en tiempo continuo de acuerdo a la siguiente ecuación diferencial
estocástica de Itô. (ver apéndice A).
Se asume que existen las matrices de funciones MxM valuadas A A(t ), A R MxM , y
el conjunto de las funciones del tiempo Bj=Bj(t) pertenecientes a MxM, para
j=1,2,…,M de modo que las ecuaciones diferenciales en derivadas parciales (2.13) y
(2.14), tengan una solución Rn valuada θ = θ(i,t) para todo i y t de tal forma que se
cumple que:
22
1
Qr T f A (2.13)
t x 2
T M
1
Q j Bj (2.14)
2 x x j 1
x, z, t CT z, t x, t (2.15)
En la cual
T
r log p x, t , 1 x, t ,..., M x, t ,Q GG T (2.16)
x
Y
2
T j
1 x, t ,..., M x, t , j Tr Q (2.17)
x2
Y definiendo
1 T 1
h x, t R t h x, t hT x , t R 1
t z (2.18)
2
d t
AT t t t C z, t (2.19)
dt
T
En el que 1 , 2 ,..., M con j
T
B j , para j=1,2,…,M y la condición
inicial de la ecuación (2.19) es ψ(t0)=0.
Teorema:
T
p x, t Z t p x, t exp x, t (t ) (2.20)
Se asume que la evolución de los estados del sistema en tiempo discreto obedece a la
siguiente ecuación estocástica en diferencias:
xk f k ( xk 1 ) 1k xk 1 , wk 1
(2.21)
yk hk ( xk ) 2k xk 1 , vk 1
(2.22)
Donde los estados xk son un proceso de Markov de orden 1, por eso cumplen
24
p( xk / x0:k 1 ) p( xk / xk 1 )
(2.24)
p( xk / y1:k 1 ) p( xk / yk 1 )
(2.25)
p( yk / xk ) p( xk / y1:k 1 )
p( xk / y1:k 1 )
p( yk / y1:k 1 )
(2.26)
2.2.1 Introducción
n
Este filtro, trata de estimar el estado x de un proceso controlado de tiempo
discreto, gobernado por la ecuación estocástica en diferencias:
xk Axk Buk wk
1 (2.28)
yk Hxk vk
(2.29)
w (t ) Q S
E w T ( ) vT ( ) (t ) (2.30)
v(t ) ST R
Q 0, R 0
o lo que es equivalente
E [ w (t ) w T ( )] Q (t ) , Q 0 (2.31)
E [v(t ) v T ( )] R (t ) , R 0 (2.32)
E [ w (t ) v T ( )] S (t ) (2.33)
1 para t
Donde (t ) es la delta de Kronecker: (t )
0 para t
Se define xˆk n
(súper menos) la estimación del estado a priori en cada tiempo k,
dado que se conoce el valor del estado previo en el tiempo k. Sea xˆk n
la
estimación del estado a posteriori en el tiempo k dadas las medidas yk . Ahora bien,
los errores de estimación a priori y a posteriori son:
ek xk xˆk ek xk xˆk
(2.35)
Pk E ek ek T
(2.36)
Pk E ek ekT
(2.37)
Debe encontrarse una ecuación que calcule el valor de la estimación del estado a
posteriori xˆk n
como combinación lineal de la estimación del estado a priori
xˆk n
y una diferencia ponderada entre la medida actual zk y la predicción de la
medida Hxˆk n
de la forma:
predicha Hxˆk n
y la medida actual o valor verdadero. La matriz Kk de dimensión n
x m en la ecuación (2.38) es la ganancia o factor de combinación que minimiza la
ecuación de la covarianza del error a posteriori. Dicha ganancia se calcula de la
siguiente forma:
1
Kk Pk H T HPk H T Rk
(2.39)
Con
lim K k H 1 , lim K k 0
Rk 0 Pk 0
(2.40)
E xk xˆk
T
E xk xˆk xk xˆk Pk
(2.41)
T
p xk zk N E xk , E xk xˆk xk xˆk
p xk zk N xˆk , Pk
(2.42)
El filtro de Kalman estima estados dentro de un proceso usando una forma de control
por realimentación: el filtro estima el estado del proceso y, al mismo tiempo, obtiene
la realimentación en forma de medidas con ruido. Entonces, las ecuaciones para el
filtro de Kalman se dividen en dos grupos: ecuaciones de actualización y ecuaciones
de medida. Las primeras son las que proyectan o adelantan en el tiempo el estado
actual y las estimaciones del error de covarianza para obtener las estimaciones a priori
para el nuevo tiempo de paso. Las ecuaciones de actualización de las medidas son las
30
dx(t ) f ( x (t ), t )dt dz
t (2.43)
yt h xs (t ) ds wt
0
Luego, nótese que las funciones de la forma P / Paa , Paa , P , f , V se evalúan sobre
at , t en el tiempo t . También para las funciones del tiempo de la forma
L at , t F t dF es un diferencial estocástico total que cumple
dL dF L at dt , t dt L at , t .
2
dqt d P / Paa dP / Paa PdPaa / Paa2 P dPaa / Paa3 dP dPaa / Paa2
(2.45)
Donde, en virtud del lema de Itô, los productos de los diferenciales se reemplazan por
los valores esperados de los productos.
El Filtro de Kalman Extendido (en inglés Extended Kalman Filter, EKF) se aplica en
sistemas dinámicos no lineales, y consiste básicamente en linealizar las ecuaciones del
sistema utilizando una metodología similar al desarrollo de primer orden de las series
de Taylor. La estimación de los estados puede obtenerse mediante las ecuaciones del
filtro de Kalman aplicadas al sistema lineal resultante de linealizar alrededor del
estado actual, usando las derivadas parciales de las funciones del proceso y de las
n
salidas medidas xk .
n
Se asume que los estados del sistema xk son gobernados por la ecuación
estocástica en diferencias
xk f xk 1 , uk , wk
(2.48)
33
m
con una salida medida yk obtenida a partir de
yk h xk , vk
(2.49)
xk f xˆk 1 , uk , 0
yk h xk , 0
(2.50)
Es importante decir que el defecto fundamental del Filtro de Kalman Extendido radica
en que las distribuciones (o las funciones de densidades en el caso continuo) de
algunas variables aleatorias, no son Gaussianas bajo sus transformaciones no lineales
respectivas. Es posible afirmar que el Filtro de Kalman Extendido es simplemente un
estimador que aproxima la optimalidad de la regla de Bayes por medio de la
linealización. Ahora bien, para estimar un proceso que se rige por medio de relaciones
no lineales, deben escribirse las ecuaciones que linealizan dichas expresiones:
xk xk A xk 1 xk 1 Wwk 1
yk y k H xk xk Vvk
(2.51)
fi
Ai , j xk 1 , uk , 0
xj
(2.52)
fi
Wi , j xk 1 , uk , 0
wj
(2.53)
34
hi
Hi, j xk , 0
xj
(2.54)
hi
Vi , j xk , 0
vj
(2.55)
Debe recordarse que en la práctica el valor real del vector de estados no está
disponible, ya que es la cantidad que se trata de estimar. Pero, por otro lado, se tiene
acceso a la medida real, la cual se usa para obtener el vector de estados estimado.
Se considera entonces un proceso del error con las siguientes ecuaciones de estado:
exk A xk 1 xk 1 k 1
ez k Hexk k
(2.56)
Note que las ecuaciones del proceso del error son lineales y se aproximan en su forma
a las ecuaciones usadas en el filtro de Kalman discreto, lo que motiva a usar este
modelo para estimar el error de predicción de los estados exk , usando el residuo de las
medidas ezk , y esta estimación denotada por eˆk , puede usarse para obtener la
estimación del estado a posteriori del proceso no lineal original
xˆk xk eˆk
(2.57)
p( k ) N 0,WQkW T (2.58)
p( k ) N 0,VRkV T
eˆk Kk eyk
(2.59)
La ecuación anterior puede usarse para actualizar las medidas en el filtro de Kalman
extendido a partir de la aproximación sin ruido e yk y la expresión de la ganancia de
Kalman del caso lineal. Resumiendo las ecuaciones de actualización en el tiempo
toman la siguiente forma:
xˆk f xˆk 1 , uk , 0
Pk Ak Pk 1 AkT Wk Qk 1WkT (2.61)
1
Kk Pk H kT H k Pk H kT Vk RkVkT
xki : i 1,..., N s
(2.63)
Donde x representa el vector de estados. Generalmente, el método de integración por
grillas, plantea que la aproximación se hace de la forma:
N
g ( xt )dxt g ( xk(i ) ) n
n i 1
(2.64)
Donde g ( xt ) es una función no lineal de los estados. Usando una grilla regular, con
el volumen ∆n, dentro del conjunto de muestras. El error de aproximación depende del
tamaño de la grilla, o la dimensión del espacio de estados. La grilla debe ser lo
suficientemente densa para dar una buena aproximación al espacio de estado continuo,
el cual debe estar predefinido.
p( xk(i ) yk ) k
1
pek ( yk h( xk( j ) )) p( xk(i ) yk 1 )
N
p( xk(i )1 yk ) pwk ( xk(i )1 f ( xk( j )1 )) p( xk(i ) yk ) n
j 1
(2.65)
Donde la notación para los estados x, las salidas y y las funciones no lineales de los
estados f y de las salidas h se mantiene. Y con el factor de normalización
i 1 (2.66)
N
xk k E ( xk yk ) xk(i ) p( xk(i ) yk ) n
i 1 (2.67)
37
Julier, Uhlmann, and Durrant-Whyte (Julier et al., 1995) argumentaron que con un
número fijo de parámetros debe ser más fácil aproximar una distribución Gaussiana
que una función no lineal arbitraria, e introdujeron una nueva generalización del filtro
de Kalman, la cual llamaron filtro de Kalman no perfumado, en inglés, Unscented
Kalman Filter (UKF) (se le llamará el filtro de Kalman Unscented, ya que la
traducción no es apropiada para su significado). Este filtro permite la propagación de
las medias y las covarianzas a través de las ecuaciones del sistema no lineal.
xk 1k 1
a
x k 1k 1
vk 1 (2.68)
k 1
xk 1k 1
a qx1
E x k 1k 1
0 (2.69)
0rx1
Donde 0qx1 y 0rx1 son vectores de ceros y q y r son las dimensiones de las matrices de
covarianza Q y R de los estados y la salida ruidosa del sistema, respectivamente. El
superíndice a significa vector ampliado.
Pk 1k 1 0nxq 0nxr
Pak 1k 1
0qxn Qk 1 Pk 1
vw
(2.70)
rxn wv
0 Pk 1 Rk 1
vw wv
Donde n es la dimensión del vector de estados original y Pk 1 y Pk 1 son las
correlaciones entre los ruidos del sistema y la salida medida.
*
En el paso siguiente, un conjunto de 2(n+q+r)+1 puntos sigmas simétricos se
calculan de la siguiente forma:
*
k 1k 1
0, n q r K Pka 1 k 1 , n q r K Pka 1 k 1
(2.71)
Para cualquier distribución previa simétrica con k, elegido tal que 0 < n+q+r+K ≤ k,
las predicciones de las medias y covarianzas son más confiables que las del Filtro de
Kalman extendido y, particularmente para distribuciones de tipo Gaussianas, debe
usarse n + q + r + K = 3. Se sabe que esta matriz de raíz cuadrada no tiene solución
única, sin embargo, las propiedades del algoritmo se mantienen para cualquier
elección de esa forma. Julier y sus colaboradores, recomiendan la descomposición de
39
a a
* *
k 1k 1 k 1 k 1,1
xk 1k 1 ,..., k 1 k 1,2 n a 1
xk 1k 1
(2.72)
kk 1
f k 1k 1
, uk 1
(2.73)
2 na 1
xk k 1 Wi k 1k 1
i 1 (2.74)
k
,i 1
na k
Wi (2.75)
1
, de otra forma
2 na k
2 na 1 T
Pk k 1
Wi k 1k 1
xk 1k 1 k 1k 1
xk 1k 1
i 1 (2.76)
yk h kk 1
, uk
a
2n 1
yk Wi yi , k (2.77)
i 1
2 na 1 T
Py Wi yi , k yk yi , k yk
i 1
2 na 1 T
Pxy Wi i,k k 1
xk k 1 yi , k yk (2.78)
i 1
xˆk k xˆk k 1
K k yk h xˆk , 0 (2.80)
y la matriz de covarianza
Pk k Pk k 1
K k Py K T k
(2.81)
Como parte del desarrollo de esta tesis, se desarrolló una aplicación de algunas de las
técnicas de estimación presentadas anteriormente, como lo son el filtro de Kalman y el
filtro de Kalman extendido, a un proceso de fermentación alcohólica en continuo
(Echeverry et al, 2003) (a definir específicamente en el Capítulo 6), superando los
resultados del trabajo de (Quintero et al, 2004). Los resultados obtenidos de la
aplicación del filtro de Kalman fueron reportados en (Quintero y di Sciascio, 2005)
del cual puede afirmarse: en el EFK la distribución de los estados es aproximada por
variables aleatorias de tipo Gaussianas, las cuales se propagan a posteriori a través de
la linealización de primer orden del sistema no lineal, este procedimiento puede
introducir errores grandes en las media y varianza a posteriori y verdaderas de las
variables aleatorias, lo que generalmente lleva a un desempeño subóptimo y algunas
veces a la divergencia o inestabilidad del filtro (Anderson and Moore, 1979; Maybeck,
1982; Sorenson, 1985).
42
Capítulo 3
MÉTODOS SECUENCIALES DE MONTE
CARLO Y FILTROS DE PARTÍCULAS
Estado del Arte
En las últimas décadas las técnicas de Monte Carlo han atraído muchísima atención de
algunos investigadores y han sido aplicadas en diferentes áreas, en las cuales se usan
el muestreo estadístico y técnicas de estimación para evaluar soluciones a problemas
matemáticos intratables analíticamente. A continuación se presentan las tres categorías
en las cuales estas técnicas se encuentran clasificadas:
E p ( xt / yt ) g ( xt ) p ( xt / yt )dX t
(3.1)
Donde g ( xt ) una función de los estados, que componen el vector X t . Esta función
que define el valor esperado, contiene la información de la mejor estimación de los
43
f ( x ) dP ( x ) (3.2)
X
Para estudiar a fondo la teoría de filtrado no lineal, se citará cada uno de los métodos
de muestreo y su respectiva inclusión en los algoritmos de filtrado, así mismo, se
analizarán las ventajas y desventajas que posee cada uno de ellos, según la literatura
analizada. Luego de cada etapa de estudio se realizará la escritura del código de cada
algoritmo planteado.
La idea del muestreo por importancia es elegir una función de distribución propuesta
q(x) (en la región de interés) denominada función de importancia en lugar de la
función de densidad de probabilidad (de masa) real p(x), la cual es difícil de
muestrear. Se asume que el soporte de la función q(x) debe cubrir a la función p(x)
(ver apéndice A).
Np
1
I W ( xi ) f ( xi )
Np i 1
(3.4)
Donde los pesos que ponderan dicha suma están determinados por:
p( xi )
W ( xi )
q( xi )
(3.5)
Np
p( xi )
(3.6)
i 1 q( xi )
Np
W ( xi ) 1
i 1 (3.7)
Np
1
W ( xi ) f ( xi ) Np
Np i 1
I Np
W ( xi ) f ( xi ) (3.8)
1 j i 1
W (x )
Np j 1
W ( xi )
W ( xi ) Np
(3.9)
j
W (x )
j 1
1
Var I Varq f ( x )W ( x )
Np
(3.10)
1 p ( x)
Var I Varq f ( x)
Np q( x)
(3.11)
2 2
1 f ( x) p ( x) E p f ( x)
Var I dx
Np q( x) Np
(3.12)
El tipo de muestreo de importancia señalado como (3.8) incluye una desviación, pero
a su vez es consistente7, ya que la desviación desaparece rápidamente a una velocidad
proporcional al número de muestras o de partículas Np.
Eq W ( x ) f ( x )
I
Eq W ( x )
(3.13)
7
Un estimador es consistente si el estimador converge al valor real casi seguro cuando el número de
observaciones se aproxima a infinito.
47
N p f Ep f N 0, f
f Varq W ( x ) f ( x ) E p f ( x )
(3.14)
f
medida de efectividad
Varp f
(3.15)
Por otro lado, la función de importancia, debe tener un gran soporte y ser insensible a
los errores. Usualmente las distribuciones súper gaussianas tienen la función de
activación de la forma:
d log q ( x)
( x)
dx (3.16)
Si ésta es acotada, q tiene una cola larga, de otro modo no. El muestreo por
importancia puede mezclarse con métodos de grilla o con el algoritmo de Metrópolis
Hastings para mejorar las técnicas.
Las muestras tomadas del muestreo por rechazo son exactas y la probabilidad de
aceptación de una variable aleatoria es inversamente proporcional a la constante C.
Ahora bien, en la práctica, la elección de la constante C es muy crítica, dado que se
basa en el conocimiento de la función p(x): si C es muy pequeña, las muestras no son
confiables ya que hay una tasa de aceptación alta; si C es muy grande, el algoritmo se
vuelve ineficiente ya que la tasa de aceptación es baja.
p( y / x) p( x) Cq( x)
p( x / y ) C ' q( x)
p( y ) p( y )
(3.17)
Donde la tasa de aceptación para cada x que pertenece al vector de estados, se puede
calcular como
p( x / y ) p( y / x)
C ' q( x) C
(3.18)
Como se mencionara anteriormente, las muestras tomadas del muestreo por rechazo
son exactas, pero el prerrequisito de este tipo de muestreo es el conocimiento previo
de C, el cual en ciertas ocasiones no se conoce. El algoritmo toma algún tiempo
cuando la tasa p(x)/(Cq(x)) es cercana a cero. En la siguiente Figura (3.1) se observa
una comparación de las técnicas anteriormente descritas. Puede observarse que para el
muestreo por importancia (izquierda), la función q(x) no necesariamente cubre en su
totalidad a la función p(x) lo que puede representar para algunas de las partículas que
se encuentren en la zona de diferencia sea asignado un peso no apropiado; en cambio,
para el muestreo por rechazo, la constante C asegura el cubrimiento de la función p(x)
49
logrando de este modo que las partículas se ubiquen lo mas cerca posible de la
distribución objetivo:
Figura 3.1 Comparación del muestreo por Importancia y muestreo por rechazo.
Ventajas y desventajas del muestreo por rechazo: La ventaja de este método es que la
partícula resultante es una muestra verdadera e independiente, obtenida de la función
de densidad a aproximar. El método compensa naturalmente los errores que se
pudieran generar cuando la densidad propuesta no sea similar a la densidad objetivo.
En este caso, el promedio de la probabilidad de aceptación será pequeño, y el
algoritmo generará más puntos a partir de la densidad propuesta. Una desventaja de
los métodos de aceptación y rechazo es que para problemas de tiempo real, es mas
ventajoso tener un método que no posea tiempos fluctuantes para completar un paso
completo. Usando estos métodos, pueden obtenerse resultados con pocos números de
partículas aceptadas en el tiempo que el filtro deba generar las muestras del próximo
tiempo de muestreo.
Para lograr la eficiencia del muestreo de importancia es necesario tener una buena
distribución propuesta, de ahí que la elección correcta de la función de distribución de
importancia q(x), es la clave para aplicar exitosamente el muestreo por importancia.
Sin embargo, es muy difícil encontrar una función de importancia adecuada para
espacios de grandes dimensiones. Una forma natural de atenuar este problema es
construir una distribución propuesta secuencialmente, los cual es la idea de este tipo
de muestreo por importancia secuencial, en inglés Sequential Importance Sampling
(SIS) que se expondrá a continuación
k
q( x0:k / y0:k ) q( x0 ) q( xi / x0:i 1 , y0:i ) (3.19)
i 1
p ( x0:k ) p ( x0 ) p ( x1 / x0 )... p ( xk / x0 , xk 1 )
(3.20)
q ( x0:k ) q0 ( x0 )q1 ( x1 / x0 )...qk ( xk / x0 , xk 1 )
Para pasos de muestreo 0:k. Cuyos pesos pueden ser escritos de la forma:
51
p( x0 ) p( x1 / x0 )... p( xk / x0 , xk 1 )
W ( x0:k ) (3.21)
q0 ( x0 )q1 ( x1 / x0 )...qn ( xk / x0 , xk 1 )
p( xk / x0 , xk 1 )
Wk ( x0:k ) Wk 1 ( x0:k ) (3.22)
qk ( xk / x0 , xk 1 )
Dado que la varianza no condicional de los pesos de importancia crece con el tiempo,
el algoritmo presenta el llamado problema de degeneración, el cual hace que en pocas
iteraciones del algoritmo, algunos de los pesos W(x(i)) sean cero (se busca que sean
todos diferentes de cero). Esto representa una desventaja si se gastan grandes
esfuerzos computacionales para actualizar pesos triviales. Para evitar estas situaciones
se usa el paso de re muestreo antes de la normalización de los pesos.
8 p( xt )
La ecuación de Focker Plank Kolmogorov (FPK) está dada por Lt p( xt ) la cual puede ser
t
interpretada de la siguiente forma:
52
Siempre que la pdf pueda ser aproximada por un histograma de puntos de masa, por
medio de un muestreo aleatorio del espacio de estados, se obtendrá un número de
partículas que representan la pdf que está evolucionando. Además, en caso que el
modelo de la densidad de probabilidad a posteriori sea desconocido o difícil de
muestrear, puede elegirse otra distribución para lograr un muestreo eficiente.
Np
p( x0:k / y1:k 1 ) Wki ( x0:k x0:i k )
i=1 (3.23)
Donde los pesos Wki se eligen por muestreo por importancia. Las suposiciones acerca
de las funciones de densidad de probabilidad involucradas son las siguientes: la
densidad de probabilidad de evaluación de las muestras es proporcional a la función
de importancia elegida p( x) q( x) y el vector x i de las muestras serán generadas a
partir de la función de densidad de importancia x i q( x) .
El primer término es la ecuación de movimiento para una nube de partículas con distribución p( xt ) , cada
dx
punto de la cual obedece la ecuación de f ( xt , t ) .
dt
El segundo término describe la perturbación provocada por un movimiento Browniano. La solución de la
ecuación de FPK puede ser encontrada usando la transformada de Fourier, ya que invirtiéndola, puede
obtenerse una distribución gaussiana de una trayectoria determinística.
53
Np
p( xk / y1:k ) Wki ( xk xki )
i=1 (3.25)
Np Np
xki , Wki SIS xki 1 ,Wki 1 , yk
i 1 i 1
For i 1: Np
xki q xk / xki 1 , yk
p ( yk / xki ) p( xki / xki 1 )
Wki Wki 1
q ( xki / xki 1 , yk )
endfor
Degeneración: este fenómeno ocurre cuando una partícula, luego de unas cuantas
iteraciones no tiene un peso considerable, lo que implica que la varianza de los pesos
pueda crecer. La degeneración de los pesos proporciona aumento en el costo
computacional, ya que se siguen considerando partículas (de poco peso) que no
aportan a la estimación de la pdf a posteriori. Para evitarla, existe una medida factible
de la degeneración de los pesos, en función de la varianza de los pesos verdaderos, la
cual puede ser modificada por una estimación cercana del tipo (Arumpalam et al.,
2002):
ˆ 1
Neff Np
2
Wki
i 1 (3.26)
54
9
Como se verá en este capítulo, los Jump Markov Linear Systems, son usados en seguimientos de objetivos
maniobrables (Doucet, Gordon, and Krishnamurthy, 2001) donde en el estado en modo discreto (que define
el índice de maniobrabilidad) se sigue usando filtros de partículas y el estado continuo (condicionado al
índice de maniobrabilidad) se sigue usando el filtro de Kalman.
55
q xk xki 1 , yk p ( xk xki 1 )
(3.30)
Y esta, es una de las formas más sencillas e intuitivas de elegir esta función tan
importante en el diseño del filtro de partículas.
3.3 Re muestreo
El paso de re muestreo esta orientado a eliminar las muestras con pesos de menor
importancia y duplicar las muestras con mayores pesos. Una breve descripción es la
siguiente:
Np
i
Se toman N p muestras aleatorias x a partir de la distribución propuesta
i 1
q(x).
i i
Se calculan los pesos de importancia W p( x) / q( x) de cada muestra x .
i
Se normalizan los pesos de importancia W
56
Np
i
Se remuestrea N veces, a partir de el conjunto discreto muestreado x ,
i 1
El paso de re muestreo puede ser aplicado en cada paso del algoritmo de muestreo de
importancia, si se considera necesario; en el cual las partículas y los pesos de
importancia asociados se reemplazan por muestras de igual importancia. Pero la
inclusión en el algoritmo de muestreo por importancia fue justificada por Liu. El paso
de re muestreo juega un papel importante en el muestreo por importancia, ya que si los
pesos de importancia están irregularmente distribuidos al propagar los pesos triviales
gasta muchísimo la potencia computacional; y por otro lado cuando los pesos de
importancia están sesgados, el re muestreo proporciona una oportunidad de
seleccionar las muestras importantes y de algún modo rejuvenecer el conjunto de
muestras para el uso siguiente. Es probable que el re muestreo no mejore
necesariamente la estimación del estado actual, ya que introduce una variación de
Monte Carlo extra (Liu, et al 2001).
En este tipo de filtro, que se abrevia con sus siglas en inglés SIR, se aplica el paso del
re muestreo para reducir los efectos de la degeneración. Este procedimiento debe
hacerse observando dicha función de degeneración N̂eff y comparándola con un
valor base.
Aplicando esta técnica, debe tenerse en cuenta que el hecho de hacer un re muestreo,
Np
implica generar un conjunto nuevo de muestras xni* por re muestreo o
i 1
Np
p( xn / y1:n ) Wni ( xn xni )
i=1 (3.32)
1
Wni
Np
(3.33)
Np Np
xkj * ,Wk j , i j RESAMPLE xki ,Wki
j 1 i 1
Inicializa el contador
c1 0
For i 2 : Np
Construye la CDF
ci ci 1 Wki
endfor
Comenzar al inicio de la CDF
i 1
Generar el punto de inicio
1
u1 U 0,
Np
For j 1: Np
Moverse a través de la CDF
( j 1)
u j ui
Np
While u j ci
i i 1
end While
Asigna la muestra
xki* xki
Asigna peso
1
Wk j
Np
Asigna parentezco
ij i
endfor
partículas con grandes pesos, se replican y las de pesos pequeños se eliminan. Por
ejemplo, la partícula azul con el mayor peso, se replica 3 veces, la partícula amarilla
se replica 2 veces, mientras que las partículas verdes que tiene pesos pequeños se
eliminan en el paso de re muestreo. Además, después del re muestreo los círculos que
representan las partículas tienen el mismo tamaño, lo que significa, que todos los
pesos se hacen iguales a 1/Np =1/10.
Continuando con la definición del filtro con re muestreo SIR, puede afirmarse que es
un filtro de tipo SIS con la característica que la densidad de importancia, es
idénticamente igual a la densidad de probabilidad a priori del vector de estado.
Np
Ni N p , Ni 0 i 1,..., N p 1
i 1 (3.35)
E Ni N pW (i )
(3.36)
Var N i N pW ( i ) 1 W ( i )
(3.37)
3.3.3.2 Re muestreo Residual: este método fue sugerido por Liu y Chen (Liu y Chen,
1998) y tiene como característica que es menos costoso computacionalmente que el
SIR convencional, garantiza una menor varianza de muestreo y no introduce una
desviación adicional; además de ello, cada una de las partículas del re muestreo
residual se replica en el algoritmo.
64
Var N i E Ni E Ni
(3.38)
Una descripción en pseudocódigo de este algoritmo es la siguiente:
Sin pérdida de generalidad, se supone que ( x)q x, x ' ( x ')q x ', x , (donde la
función π indica una función de distribución no normalizada). La ecuaciçon anterior
significa que la probabilidad de moverse de x a x ' es mayor (es decir, es mas
frecuente este movimiento) que la probabilidad de moverse de x ' a x .
( x ')q x, x '
min ,1 , si ( x ')q x, x ' 0
x, x ' ( x ) q x, x ' (3.40)
1 , de otro modo
10
Ver apéndice A
66
Resumiendo el algoritmo:
Por otro lado, las partículas muestreadas son recordadas como las muestras tomadas a
partir de la densidad objetivo y solo después de que la cadena ha pasado la fase
transitoria, la convergencia de la distribución invariante ocurre bajo malas condiciones
de regularidad como la irreductibilidad y la no-periodicidad. Una de las características
de este algoritmo es que su eficiencia está determinada por la tasa de muestras
aceptadas respecto al número total de muestras. Una varianza muy grande o muy
pequeña, puede generar un muestreo ineficiente.
Np Np
p( x) ci pi ( x), ci 1
i 1 i 1 (3.43)
Np
Mi N p , (N p Mi )
i 1 (3.44)
p( xk / xk(i )1 ) p( yk / xk )dxk
ci Np p ( xk / xk( i )1 ) p ( yk / xk )
p( xk / x (i )
) p( yk / xk )dxk pi ( xk )
k 1 p ( xk / xk( i )1 ) p ( yk / xk )dxk
i 1 (3.45)
69
Para el i-ésimo estrato, los pesos de importancia asociados con las partículas se
actualizan recursivamente por la siguiente expresión:
Intuitivamente usaron una medida ponderada antes del re muestreo, en vez que usar
una medida no ponderada; debido a que se espera que las muestras ponderadas
contengan más información que el mismo numero de muestras sin ponderar. En el
instante de re muestreo se selecciona una muestra de tamaño Np a partir de los 10 x Np
valores predichos con el fin de hacer que el tamaño del conjunto de partículas no
cambie. El numero 10, fue sugerido por Rubin (Rubin, 1987) donde Np>>10. Se
asume que el número de conjunto de partículas no cambia.
Como referencia, puede tenerse en cuenta que (Carpenter et al, 1999) desarrollaron un
algoritmo de complejidad de orden O(Np) vía muestreo estratificado, aprovechando las
ventajas del método de simulación de las estadísticas de distintos órdenes. Muchos
filtros de partículas mejorados, se deben a la implementación de un paso de re
muestreo. A continuación un tipo de algoritmo sugerido:
Np
Dado un conjunto de Np partículas xk( i ) , Wki , se sugiere crear que nuevo
i 1
Np
conjunto de partículas independiente del anterior xk( j ) ,Wk j , en el instante
j 1
de re muestreo
Para j=1,…,Np, los valores de xk( j ) reeemplacen a xk(i ) con probabilidad
proporcional a a(i ) .
Los nuevos pesos asociados Wk j se actualizan de la siguiente forma :
Wk j Wki / a(i ) .
Nx
f ( x)
(3.47)
' 1
f (fa f b) (3.48)
2
Donde
2Vara f 2Varb f
Vara f Varb f (3.50)
Np Np
N pVar f Var f
2
Vara f Varb f Ea f Eb f
N pVar f
2 4
2
' Ea f Eb f
N pVar f N pVar f
4
'
N pVar f N pVar f
(3.51)
'
Donde la varianza estratificada Var f nunca es mayor que la varianza convencional
Este tipo de filtro se plantea como una solución al problema de pérdida de diversidad
en las partículas debido al re muestreo (donde las muestras se obtienen de una
distribución discreta). Este re muestreo puede llevar a un colapso en las partículas, o a
un posible caso severo de Sample Impoverishment, donde las Np partículas ocupan un
punto en el espacio. Este tipo de filtro es idéntico al SIR, excepto por la fase de re
muestreo donde se re muestrea a partir de una densidad posterior continua aproximada
del tipo:
Np
p( xk / y1:k ) Wki K h xk xki
i 1 (3.52)
Donde los Wki son los pesos normalizados y la función:
1 x
Kh X K
h Nx h (3.53)
2
xK x dx 0, x K x dx
(3.54)
d ( K h v)
( x) K h ( x u )v(du )
dx (3.55)
y con (*) el operador convolución. Si la densidad de probabilidad v puede expresarse
en términos de las deltas de Dirac, quiere decir que es una densidad de probabilidad
Nx
discreta en y los valores de los estados son muestras tomadas a partir de ella.
Luego,
Np Np
d ( K h v) 1 1
( x) w(i ) K ( ( x x (i ) )) w( i ) K h ( x x ( i ) )
dx h nx i 1 h i 1 (3.56)
K se elige para minimizar el error integral medio cuadrático (en inglés mean
integrated square error MISE) entre la densidad posterior verdadera y la
representación empírica regularizada. En caso que las muestras tienen el mismo peso,
73
existe un kernel óptimo y un ancho de banda óptimo, de la misma forma que una
densidad Gaussiana. De esta forma pueden definirse:
nx 2 2
(1 x ), si x 1
K opt 2cnx
0, de otro modo
(3.57)
Nx
Donde cnx es el volumen de la hiper esfera unitaria en . Cuando la densidad
adyacente es gaussiana, con matriz de covarianza unitaria, el ancho de banda óptimo
esta dado por
hopt AN p1/( nx 4)
(3.58)
1/( nx 4)
1 nx
A 8cnx (nx 4)(2 )
(3.59)
Luego, se aplica una transformación lineal para lograra una covarianza unitaria (este
proceso se llama blanqueamiento, en inglés whitening). Es decir, la muestra x(i ) , se
cambia a A 1 x(i ) , donde AAT S . La matriz de covarianza del nuevo sistema de
partículas es la matriz unitaria, y el ancho de banda se puede usar directamente!. Con
esto, el problema se reduce a usar el siguiente Kernel de regularización:
(det A) 1 1
K ( A 1 x)
hnx h (3.60)
Para poder manejar densidades multi modales, se elige h=hopt/2. Existen dos tipos de
filtros de partículas regularizados, dependiendo en que momento del algoritmo de
filtros de partículas se ubique el paso de regularización: si antes o después del paso de
corrección.
Fue propuesto por Musso y Oudjane en (Musso et al., 2001). Una simple
aproximación a este tipo de filtro es la que se presenta a continuación en tres pasos
que resumen el algoritmo:
74
Básicamente cumple con la base de los filtros de regularización, solo que se realiza
este paso en un momento diferente del algoritmo. Un análisis teórico de este tipo de
filtro puede encontrarse en la publicación de Legland y colaboradores (Legland et al,
1998), los pasos de este algoritmo son los siguientes:
Ambos filtros convergen al filtro óptimo en sentido débil, con una tasa de h2 1/ N .
El término h2 corresponde al error debido a la regularización. Cuando h=0, se alcanza
la tasa de convergencia del filtro clásico de regularización. En el sentido fuerte del
espacio L1, el error estimado es proporcional a h2 1/ Nhnx para los filtros post y
pre regularizados. De igual forma que para el filtro clásico, estos estimativos del error,
generalmente no están acotados en el tiempo. A continuación, se presenta el modelo
de algoritmo propuesto para el filtro de partículas regularizado (Arumpalam et al.,
2001).
76
Np Np
xni , Wni RPF xni 1 ,Wni 1 , yn
i 1 i 1
For i 1: Np
xni q xn / xni 1 , yn
p ( yn / xni ) p( xni / xni 1 )
Wni Wni 1
q ( xni / xni 1 , yn )
end
Np
PesoTotal Wni
i 1
For i 1 : Np
Wni
Wni
PesoTotal
end
^ 1
Neff Np
2
Wni
i 1
^
if Neff NT
Sk Cov xni ,Wni
se calcula Dk tal que S k Dk DkT
REMUESTREAR : usando el método más apropiado
For i 1: Np
calcular el kernel K
xni* xni hopt Dk K
end
endif
procedimientos de muestreo por importancia. La cualidad de esta técnica que mas nos
interesa es su propiedad intrínseca de reducción de la varianza, la cual se usa en
filtrado de partículas para mejorar su desempeño.
fˆ ( y) E p ( x ) fˆ Y (Y ) ( y)
(3.61)
Varp ( x ) fˆ ( y) Varp ( x ) fˆ Y
(3.62)
P fˆ Y fˆ ( y) 1
(3.63)
Definiendo x0:j k x0j ,..., xkj , puede lograrse una expresión para el posteriori de la
siguiente forma:
78
f k x10:k , x0:2 k p y0:k / x0:1 k , x0:2 k p x0:1 k /, x0:2 k d x10:k p x10:k d x10:k
E( fk )
p y0:k / x0:1 k , x0:2 k p x0:2 k / x0:1 k d x0:2 k p x0:1 k d x0:1 k
Debe asumirse que la probabilidad condicional dada una realización de x10:k g x0:1 k y
1
p y0:k / x0: k puedan ser evaluados analíticamente, y de este modo pueden hacerse las
dos estimaciones de la densidad de probabilidad a posteriori basadas en el muestreo
de importancia.
Np
g x1,( i)
w* x1,( i)
Eˆ ( f k ) i 1 0:k 0:k
Np
w* x1,(
0:k
i)
i 1 (3.67)
Donde los pesos son
w* x1,(
0:k
i)
p x10:k / y0:k / q x10:k / y0:k
(3.68)
La siguiente proposición muestra que “si se puede integrar analíticamente una de las
componentes de la varianza del estado obtenido, es más débil que la de la estimación
simple”.
varp x1 / y0:k
w* , x10:k varp x1 , x0:2 k / y0:k
w* , x0:1 k , x0:2 k
0:k 0:k
(3.69)
puede ser usada para estimar la distribución marginal p x10:k / y0:k , pero
Np
f k x1,( i)
w x1,( i)
Eˆ ( f k ) i 1 0:k 0:k
Np
w* x1,(
0:k
i)
i 1 (3.70)
que
Np
i 1
Ep x0:2 k / y0:k . x10:k
f k x0:2 k w x0:1,(ki )
Eˆ ( f k ) Np
w x1,(
0:k
i)
i 1 (3.71)
En todos los casos, es posible usar los Métodos de Monte Carlo desarrollados
anteriormente (Doucet, 1998).
80
end
Paso de selección
Multiplicar o suprimir las muestras x0:i k con altos o bajos pesos de importancia wk( i )
respectivamente, para obtener N p muestras aleatorias x0:i k distribuidas aproximadamente
de acuerdo a p x0:i k / y1:k
Paso de Monte Carlo y Cadenas de Markov
Aplicando un kernel de transición con una distribución
dad por p x0:i k / y1:k para obtener x0:i k
k
q x10:k / y1:k q x0 q x1m / y1:m , x1:m 1
m 1 (3.72)
p yk / y1:k 1 , x0:k p xk / xk 1
wk
q xk / y1:k 1 , x1:k 1
(3.73)
Brevemente se dirá de esta técnica que se utiliza la técnica de Jittering para aliviar el
problema de empobrecimiento de las muestras es decir, en cada tiempo de muestreo se
agrega una cierta cantidad de ruido gaussiano a cada partícula remuestreada. Este
procedimiento es equivalente a usar un kernel gaussiano para suavizar la distribución
a posteriori, previene el problema de la divergencia. Más adelante se encontrarán mas
detalles de posible aplicación en los filtros de partículas.
2
p p f 1 p
f ptr tr Q 2
t x x 2 x
Q GGT , matriz de coeficientes de ruido (3.74)
Los filtros que resuelven las ecuaciones de densidad de probabilidad en tiempo real
son los llamados filtros exactos de dimensión finita, porque el vector de estados en si
mismo tiene una dimensión finita fijada previamente. Debe considerarse pues, que las
estadísticas suficientes propuestas por Fisher son las que son propagadas a través del
algoritmo y que la predicción de las se hace mediante la integración numérica de las
ecuaciones diferenciales ordinarias en la familia exponencial de las densidades de
probabilidad. Nuevamente, se aplica la regla de Bayes para la corrección de las
características estadísticas.
Los Filtros Recursivos Exactos No lineales pueden definirse brevemente así: Un filtro
es exacto si proporciona exactitud de la estimación óptima, con un algoritmo que solo
requiera la solución de una ecuación diferencial ordinaria (ODE) en tiempo real, pero
no exige la solución de las densidades de probabilidad de los estados en tiempo real.
Por ejemplo, el filtro de Kalman es exacto para sistemas lineales gaussianos.
Este tipo de filtro, esta basado en el uso de una familia exponencial de distribuciones
de probabilidad del tipo:
T
p xt / yt a x, t b y, t exp x, t y, t
(3.75)
Así pues, para solucionar este tipo de problemas, se asume que la función p no
desaparece en ningún punto, es doblemente continuamente diferenciable en x y
continuamente diferenciable en t; y además, la función p se aproxima a cero lo
suficientemente rápido como x y tiene la propiedad de satisfacer la ecuación
de Focker Planck. Ahora bien, para una condición inicial dada de la función p x, tk ,
la ecuación de Focker Planck tiene una única solución acotada para todo x, en el
intervalo tk t tk 1 .
T
p xt / yt p x, t exp x, t y, t (3.76)
Fuera de línea
2
p p f 1 p
f ptr tr Q 2
t x x 2 x
Q GGT , matriz de coeficientes de ruido (3.77)
T 1
Q f A
t x 2 (3.78)
En línea
d
AT (t ) (t ) (t )
dt
(tk ) (tk ) C ( yk , tk )
(3.79)
85
Este algoritmo fue introducido por Pitt y Shepard (Pitt y Shepard, 1997) en donde
plantean que el algoritmo típico de filtro de partículas no es completamente robusto.
Para respaldar esta afirmación inician su trabajo sustentando su falta de robustez
debido a dos razones principales: la primera es el diseño de los simuladores y la
segunda es el uso de los soportes discretos de las funciones de probabilidad, para
representar la distribución previa actualizada secuencialmente.
Esta propuesta de algoritmo, tiene como objetivo aumentar las partículas buenas
existentes dado que las probabilidades predictivas son grandes; y para ello usa la
distribución óptima en los casos que esta sea exacta, en caso contrario, usa la mejor
aproximación. Esta metodología, agrega un componente nuevo al algoritmo de
filtrado: se trata de un índice de componentes mixtas, llamado variable auxiliar, la
cual es el centro de su propuesta.
Este último paso implica que el paso de re muestreo no se hace, ya que se hace el tipo
de filtrado de un paso adelante basado en un punto de estimación. Cuando el ruido del
proceso es pequeño el desempeño es mejor que el SIR. Ahora bien, cuando el ruido del
proceso es grande, el punto estimado ki no da suficiente información de p ( xk / xki 1 ) .
Para este tipo de filtro se introduce una función de importancia q ( xk ,i , y1:k ) , que
Ns
muestrea el par xk , i j , donde ij es el índice de la partícula en el instante k-1.
j 1
p( xk ,i , yk ) p( yk / xk ) p( xk / xki 1 )Wki 1
(3.80)
p( xk / yk ) .
q ( xk ,i , y1:k ) p xk / xki 1
(3.81)
i
k E xk / xki 1
(3.83)
Np
p( xk / y0:k ) p( yk / xk ) p( xk / xk 1 ) p( xk 1 / y0:k 1 )dxk 1 Wk(i ) p ( yk / xk ) p ( xk / xk(i )1 )
i 1
(3.84)
Donde
Np
p( xk 1 / y0:k 1 ) Wk(i )1 ( xk / xk(i )1 )
i 1 (3.85)
Ahora bien, el producto Wk(i ) p ( yk / xk ) se trata como la probabilidad combinada,
atribuida a la densidad de filtrado. Al momento de introducir la variable auxiliar,
1,..., Np que juega el papel de índice del componente mixto, la densidad
conjunta aumentada p( xk , / y0:k ) se actualiza de la siguiente manera:
p( xk , i / y0:k ) p( yk / xk ) p( xk , i / y0:k 1 )
p( xk , i / y0:k ) p( yk / xk ) p( xk , i / y0:k 1 ) p(i / y0:k 1 )
(i )
p( xk , i / y0:k ) p( yk / xk ) p( xk / x )Wk(i )1
k 1 (3.86)
Np
p( xk / y0:k ) Wk(i )1 p ( yk / xk( i ) . i ) p( xk / xk( i )1 ) (3.87)
i 1
i
Donde denota el índice de la partícula xk(i ) en el instante de tiempo k-1,
Np
llámese i
i . La distribución propuesta usada para muestrear xk( i ) . i
se
î 1
Donde
i (i )
q / y0:k p yk / k Wk(i )1
q xk . i / y0:k p xk / xk(i )1
(3.89)
Np
p( xk / y0:k ) Wk(i )1 p( yk / (
k
i)
) p( xk / xk( 1
i)
)
i 1 (3.90)
De las ecuaciones anteriores, puede obtenerse una expresión para la obtención de los
pesos importantes, los cuales serán actualizados recursivamente.
q xk . y0:k
(i )
p ( yk / x )
Wk(i ) k
p( yk / k( i ) )
(3.91)
for i 1,..., Np
calcular
(i )
k E p xk xk( i )1
end
for i 1,..., Np
calcular los pesos del primer paso
Wk( i ) Wk( i )1 p yk (i )
k
end
Usar el procedimiento de remuestreo del SIR
i Np
para obtener las nuevas partículas xn( i ) ,
i 1
for i 1,..., Np
muestree
xn( i ) p xn( i ) xn(i )1 , i
Vale la pena hacer una comparación entre el filtro de partículas auxiliar y el filtro de
partículas de muestreo por importancia, sobre la eficiencia estadística en el contexto
de la medida aleatoria E W 2 xk . Los pioneros Pitt y Shephard mostraron que
cuando la probabilidad no varia sobre diferentes valores de , entonces la varianza
del ASIR es más pequeña que la del filtro tipo SIR. El ASIR puede entenderse como un
filtrado de un paso adelante en el cual la partícula xk( i )1 se propaga a k(i ) en el
próximo tiempo, para poder permitir que se haga el muestreo del posterior.
Por otro lado, el ASIR hace un re muestreo de p( xk 1 / y0:k ) en vez del re muestreo de
p( xk / y0:k ) que se hace normalmente en el SIR, lo que hace que usualmente se alcance
90
menores varianzas porque el estimado pasado es mas confiable. De este modo, el ASIR
toma ventaja de la anticipación de la información a partir del modelo de probabilidad,
para evitar el muestreo ineficiente porque las partículas con bajas probabilidades son
consideradas menos informativas. En otras palabras, las partículas a ser muestreadas
son empujadas intuitivamente a la región de más alta probabilidad.
xk x1k , xk2
(3.92)
Donde x1k y xk2 se usan para denotar las variables de estado lineales y no lineales,
respectivamente. La idea básica debajo del filtro de partículas marginalizado es partir
la probabilidad p x10:k , x0:2 k / yk de acuerdo a la siguiente expresión:
f k x10:k , x0:2 k p y0:k x0:1 k , x0:2 k p x0:1 k , x0:2 k d x0:1 k p x0:1 k d x0:1 k
E( fk )
p y0:k x0:1 k , x0:2 k p x0:2 k x0:1 k d x0:2 k p x0:1 k d x0:1 k
Si el filtro de partículas estándar se usa para todos los estados, la dimensión del
espacio en el cual las partículas residen será dimensión (n x t) =dimxt, mientras que si
se usa el filtro de partículas marginalizado, la dimensión correspondiente será
dimensión (n x nt) = dimxnt. Intuitivamente, como dimxnt < dimxt. Nótese que para
tener una mejor estimación, aplicando el filtro de partículas tradicional, deben usarse
más partículas, que las que se usarían con el filtro de partículas marginalizado. En el
93
Este tipo de filtro de partículas, el marginalizado ha sido usado con éxito en muchas
aplicaciones, por ejemplo, en seguimiento de objetivos con auto movimiento (Eidehall
et al., 2005), posicionamiento de agentes autónomos (Svenzén, 2002), aeronavegación
[Nordlund, 2002], navegación subacuática (Karlsson y Gustafsson, 2003),
comunicaciones (Chen et al., 2000), (Wang et al., 2002), identificación de sistemas no
lineales (Schön y Gustafsson, 2005), (Li et al., 2003), (Daly et al., 2005) y separación
de fuentes de audio (Andrieu y Godsill, 2000). Además, en el trabajo H en la tesis
doctoral de Schön, se describió el filtro de partículas desde el punto de vista de las
aplicaciones, utilizando muchas de ellas.
zn xn ,
(3.96)
Asumiendo que la variación de los parámetros se rige por una caminata aleatoria, el
problema de identificación de sistemas puede tratarse como un caso de estimación de
estados no lineales, lo que hace posible el uso de varios algoritmos en este tipo de
problema. El modelo dinámico resultante es:
xt 1 f1 xt , ut , t f 2 xt , ut , t wt
t 1 t wt
yt h1 xt , ut , t h2 xt , ut , t et
(3.97)
El cual es un caso especial del modelo 5 de Schön, lo que implica que el filtro de
partículas marginalizado puede usarse. Aquí, este algoritmo puede utilizarse para
obtener la solución al problema de identificación de parámetros. Los detalles de esta
aproximación se dan en el trabajo E de la tesis de Schön (Schön, 2005). También, una
aproximación similar fue desarrollada por (Li et al., 2003), (Andrieu y Godsill, 2000)
y también fue empleada por (Daly et al., 2005). Esta idea fue explorada anteriormente
por (Ljung, 1979), donde existe la salvedad que el problema de estimación resultante
fue manejado usando el filtro de Kalman. El trabajo de (Kitakawa, 1998) también es
94
Johansen basa su propuesta en los modelos de Markov ocultos (Hidden Markov Model
HMM), el estado y las medidas en un tiempo k son dos vectores, posiblemente con
dimensiones diferentes y que pertenecen a espacios E y F respectivamente. Esos
vectores evolucionan aleatoriamente a través del tiempo, pero sus dimensiones son
fijas. El propósito del filtrado es obtener un cálculo recursivo en el tiempo de la
distribución del estado oculto, dadas todas las observaciones que se hallan hecho. En
el filtrado multi objetivo, que ha sido introducido recientemente y estudiado por las
comunidades de fusión de datos y estudiosos en seguimiento (Goodman et al, 1997),
(Mahler, 2003); el objetivo es realizar un filtrado cuando las variables de estado y
observación son subconjuntos finitos de E y F. Conceptualmente, este problema puede
95
ser pensado como que se realiza un filtrado cuando los espacios de estado y
observación son uniones disjuntas, respectivamente. Es importante anotar que
desarrollar herramientas eficientes computacionalmente en este aspecto es
extremadamente difícil. (Gentil et al., 2005). Una alternativa que es fácil de aproximar
computacionalmente es el filtro de densidad de hipótesis de probabilidad (Probability
Hypothesis Density PHD), propuesto recientemente por (Mahler et al, 2003), el cual es
un algoritmo recursivo que propaga el primer momento, llamado también la intensidad
del posterior multiobjeto; el primer momento, es una medida definida apropiadamente
en el espacio E (puede usarse este término para definir también la derivada de Radon-
Nikodym de esta medida, respecto alguna medida dominante apropiadamente definida
en el mismo espacio). Mientras que el primer momento es una función en el espacio
E, la dimensión del espacio esta fija, la recursividad del filtro PHD aún envuelve el
uso de integrales múltiples que no tiene una expresión en forma cerrada en general.
Así pues, el filtro de densidad de hipótesis de probabilidad (PHD Filter), es un
método de actualización de una medida, dado un conjunto aleatorio de observaciones
que puede ser interpretada como la aproximación de primer momento de la ecuación
usual de filtrado Bayesiano. Dentro de este contexto, la cantidad de interés es la
medida de un punto de proceso. Aunque pueda ser descrita por una medida, no es una
medida de probabilidad en general, y es necesario mantener una estimación de la masa
total y la distribución de esa masa. A continuación se dará una explicación:
k ( A) : E N Ek ( A) y1 ,..., yk , k 0
(3.98)
ˆ k (dx) k k 1 dx k 1 Kk dx k dx
( x)
(dx) ˆk dx vk ( x)
ky
ˆ k (dx)
k k
y Yk kk ( y ) ˆk ( ky )
(3.99)
Kk x, dy ek ( x) f k x, dy bk x, dy
(3.100)
En el año de 2005, Klaas y sus coautores, generaron una propuesta que mejora el
trabajo realizado por Pitt y Shepard (Pitt y Shepard, 1997), en ese trabajo, ellos
presentaron un algoritmo secuencial de Monte Carlo llamado, el Filtro de Partículas
11
Uno de los cuales describe la ecuación de actualizacion del filtrado bayesiano.
97
Usando esta idea, llegaron a una versión mejorada del filtro de partículas auxiliar,
mencionado como Filtro de Partículas Marginal; mostraron resultados teóricos y
empíricos con lo cual demostraron la reducción en la varianza sobre el filtro de
partículas convencional, y presentaron técnicas para reducir el costo del filtro de
partículas marginal con N partículas del orden O(N2) al orden O(N logN).
p xk / y1:k p yk / xk p xk / y1:k 1
(3.102)
p xk / y1:k p yk / xk p xk / x1:k 1 p xk 1 / y1:k 1 dxk 1
(3.103)
98
El algoritmo que propone (Klaas et al., 2005) puede ser resumido en el siguiente
código:
end
for i 1,..., Np
Evaluar los pesos de importancia
N
p yk xk( i ) wk( j )1 p xk( i ) x1:( kj ) 1
(i) j 1
w k N
wk( j )1q xk yk , xk( j )1
j 1
end
Normalizar los pesos de importancia
wk( i )
wk( i )
wk( i )
j
Np
pˆ xk / y1:k p yk / xk wk( j )1 p xk / xk( j )1 (3.104)
j 1
Np
pˆ xk / y1:k wk( j )1 p yk / xk( j )1 p xk / xk( j )1 , yk (3.105)
j 1
99
Para llegar a la aproximación del filtro marginal, debemos retomar, del filtro de
partículas auxiliar donde se usa la siguiente distribución conjunta
(m)
wk( m1) p yk / (m)
k
k 1 Np
(3.109)
wk( j )1 p yk / ( j)
k
j 1
Donde
(m)
q m / y1:k k 1 , y, q xk / y1:k , m q xk / xk( m1) , yk
(3.111)
pˆ xk / y1:k , m
w m, xk (3.112)
q xk / y1:k , m
wk( m1) p yk / xk( m1) p xk / xk( m1) , yk
w m, xk (m)
(3.113)
k 1 q xk / xk( m1) , yk
100
N
wk( j )1 p yk / xk( j )1 p xk / xk( j )1 , yk
pˆ xk / y1:k j 1
w xk , w m, xk N
(3.114)
q xk / y1:k ( j)
k 1 q xk / xk( j )1 , yk
j 1
(i )
(i )
k 1 wk( i )1 p yk (i )
k , (i )
k 1 Np
k 1
( j)
k 1
j
end
for i 1,..., Np
Muestrear de la distribución propuesta
Np
xk( i ) ( j)
k 1 q xk yk , xk( j )1
j 1
end
for i 1,..., Np
Evaluar los pesos de importancia
Np
p yk xk( i ) wk( j )1 p xk( i ) xk( i )1
j 1
wk( i ) Np
( j)
k 1 q xk( i ) yk , xk( i )1
j 1
end
Normalizar los pesos de importancia
wk( i )
wk( i ) Np
wk( j )
j 1
Debido que
E var w m, xk xk 0
(3.117)
var w m, xk var w m, xk xk
(3.118)
El filtro de partículas marginal mejora los filtros SIR y ASIR, siempre y cuando una
partícula tenga una probabilidad marginal mayor, pero con poca probabilidad conjunta
a través de la trayectoria. Esto puede ocurrir debido a las distribuciones de
importancia con colas largas o modelos con previos de transición angostos o sin
especificaciones. Desde otro punto de vista, las mejoras del MPF sobre el SIR no
serán muy pronunciadas si el modelo de observación tiene comportamiento de pico en
ciertas regiones (por ejemplo, si la probabilidad esta bastante concentrada en algún
punto), como esto puede influenciar los pesos de importancia mucho mas que el efecto
de muestrear en el espacio conjunto. En esos casos, debería usarse el AMPF.
Finalmente, debe notarse que los métodos secuenciales de Monte Carlo aplican a
dominios fuera del filtrado Bayesiano, y, puede derivarse como consecuencia de ello,
en un algoritmo marginal de Monte Carlo secuencial análogo dentro del contexto de
los sequential monte carlo methods (SMC).
Una anotación importante es la de citar el caso en el cual los filtros SIR y MPF son
equivalentes y este es cuando el previo de la transición se usa como distribución
propuesta, entonces la ecuación de actualización de los pesos en el filtro marginal se
convierte en la siguiente:
Np
p yk xk(i ) wk( j )1 p xk( i ) xk( j )1
j 1
wk(i ) Np
(3.119)
( j)
wk 1 p xk(i ) xk( j )1
j 1
(i )
w k p yk xk( i ) (3.120)
Respecto a el costo computacional del algoritmo, fue mostrado por Klaas que puede
ser reducido de O Np 2 a ordenes de O Np log Np y en algunos casos del orden de
103
Nx
Para las siguientes definiciones debe considerarse un vector de estados x en un
espacio de probabilidad , F , P , y el kernel de transición en ese espacio, K(-,-) , el
cual representa la probabilidad que los estados se muevan de un punto x a un punto en
el conjunto S B (siendo B el álgebra de sigma álgebra de Borel de RNx). Una cadena
de Markov es una secuencia de variables aleatorias xk k 0 de modo que se
cumpla:
Pr xk B x0 ,..., xk 1 Pr xk B xk 1 (3.121)
K xk 1 , xk p xk xk 1 (3.122)
Markov es el teorema de ergodicidad, que establece bajo que condiciones, una cadena
de Markov puede ser analizada para determinar su comportamiento en estado estable.
Donde x' S Nx
y es la densidad que determina que la medida de Lebesgue de
Q sea tal que Q(dx’) = (x’)dx’. La n-ésima iteración está dada por
Kn 1
x, dx ' K ( x ' , S )
x (3.125)
(3.127)
Sin pérdida de generalidad se supone que ( x)q x, x ' ( x ')q x ', x , (donde la
función π indica una función de distribución no normalizada) lo que significa que la
probabilidad de moverse de x a x ' es mayor (es decir, es mas frecuente este
movimiento) que la probabilidad de moverse de x ' a x . Lo que se desea es cambiar
esta situación, para reducir el numero de movimientos de x a x ' . Haciéndolo, se
considera la probabilidad de movimiento 0 x, x ' 1 . Si este movimiento no se
realiza, el proceso entrega el valor de x proporcionado por la distribución objetivo.
Luego, la transición de x a x ' se convierte en:
( x ')q x, x '
min ,1
x, x ' ( x ) q x, x ' , si ( x ')q x, x ' 0 (3.129)
1, de otro modo
Resumiendo el algoritmo:
106
Una de las características de este algoritmo es que su eficiencia está determinada por
la tasa de muestras aceptadas respecto al número total de muestras. Una varianza muy
grande o muy pequeña, puede generar un muestreo ineficiente.
Es una herramienta poderosa para generar muestras aleatorias que pueden ser usadas
para calcular estimaciones estadísticas y probabilidades marginales y condicionales.
Los métodos de Monte Carlo y Cadenas de Markov (MCMC) se basan en secuencias
de Markov dependientes, que poseen una distribución limitada correspondiente a la
distribución de interés. El término “Cadena de Markov” en las técnicas de MCMC no
es completamente consistente con su uso estándar en los procesos estocásticos. Este
término se reserva generalmente para su uso con procesos que tienen salidas discretas.
E f x f x p x dx
(3.132)
n
1
f xk
n M k M 1 (3.133)
dx dy
xk 1 f xk , Ak 1 ,Wk , yk g xk , Ak ,Vk
(3.134)
control. En los modelos tipo Modelos de Markov Ocultos (Hidden Markov Models) es
primordial resolver el problema de estimar los estados ocultos, lo cual se hace
propagando la distribución posterior (o densidad de filtrado). Con una elección
correcta de la secuencia de control, la evolución del estado y de la observación del
proceso puede ser orientada con el fin de obtener densidades de filtrado que den
estimaciones mas certeras de los estados de proceso. Este problema se conoce como
planificación de sensores.
x0 p0 , xk 1 p . xk , Ak 1
(3.135)
yk q . xk , Ak
(3.136)
N
N n 2
min J A1:N E 0 , A1:N
Xn n ,
A1:N A
n 1
(3.137)
da N
Donde A es un conjunto abierto, y para cualquier 1 n N y la función
n n dy n
dx da
integrable h : .Y 0,1 es un factor de descuento y
su inclusión favorece mejor desempeño en las tareas de seguimiento en las últimas
épocas de simulación.
E 0 , A1:k
h x1:k , A1:k , y1:k :
k
(3.138)
h x1:k , A1:k , Y1:k i 1
q yi xi , Ai p xi xi 1 , Ai 0 x0 dx0:n dy0:n
da N
Con el fin de minimizar el criterio antes visto, cuando A es un conjunto
abierto y continuo sin posibilidades de discretización, la técnica de la simulación con
Aproximación Estocástica es un algoritmo de gradiente descendiente que requiere
únicamente estimaciones ruidosas del gradiente de la función de costo respecto a la
secuencia de acciones:
L
a. Inicialice con x x,i L 1
i 1
t x s x
q min i 1
,1
t x s x
p w q xw
p x, w min ,1
p x q wx
(3.140)
El método está basado en una concatenación de m algoritmos de M-H, uno para cada
variable en el vector de interés, de la siguiente forma:
px / y , z w y, z q x y , z
p x, w min ,1 min 1,1 1
px / y , z x y, z q w y , z
(3.141)
Capítulo 4
APROXIMACIÓN AL PROBLEMA DE
DISCRETIZACIÓN DE MODELOS NO
LINEALES Y MODELADO DE
INCERTIDUMBRES PARA FILTRADO NO
LINEAL
“Most control designs are based on the use of a design model. The relationship
between models and the reality they represent is subtle and complex. A
mathematical model provides a map from inputs to responses. The quality of a
model depends on how closely its responses match those of the true plant. Since no
single fixed model can respond exactly like the true plant, we need, at the very least,
a set of maps. However, the modeling problems is much deeper - the universe of
mathematical models from which a model set is chosen is distinct from the universe
of physical systems. Therefore, a model which includes the true physical plan can
never be constructed. It is necessary for the engineer to take a leap of faith
regarding the applicability of a particular design based on a mathematical model.
To be practical, a design technique must help make this lead small by accounting
for the inevitable inadecuacy of models. A good model should be simple enough to
facilitate design, yet complex enough to gives the engineer confidence that designs
based on the model will work on the true plant” Kemin Zhou with John C. Doyle and
Keith Glover. Robust and optimal control, Prentice Hall, NJ
Motivación:
En los capítulos anteriores se ha hecho una revisión de las técnicas de filtrado usadas
en estimación de estados e identificación de sistemas. Asimismo se ha proporcionado
una descripción detallada de las técnicas de filtrado no lineal y filtrado de partículas
como aplicación de los métodos secuenciales de Monte Carlo en una aproximación a
las posibles soluciones al problema de filtrado óptimo. No debe olvidarse que todas
estas técnicas que han sido aplicadas desde hace décadas y que serán utilizadas en esta
tesis para la obtención de información de sistemas dinámicos con características no
lineales no Gaussianas, tiene como objetivo la reconstrucción de la información
faltante en un sistema.
115
Debe tenerse en cuenta que aquella información recolectada, puede ser usada en la
mayoría de los casos para fines de monitoreo, optimización y/o finalmente dentro de
lazos de control para múltiples aplicaciones.
Recordado este objetivo, entonces puede hacerse hincapié en que para la ingeniería de
sistemas de control, es necesario recolectar el mayor conocimiento posible del sistema
que se va a tratar y para este propósito se ha dedicado mucho esfuerzo en la
construcción de modelos y una teoría alrededor de ello.
Guergachi y Patry (Guergachi y Patry, 2001) propusieron prestar una mayor atención
al modelado del sistema físico en sí mismo, y de esta forma modelar la incertidumbre
que rige su comportamiento.
dxt f xt , t dt 1 xt , t dwt
(4.1)
dyt g xt , t dt 2 xt , t dvt
Este modelo básico se describe simplemente en los casos que se tiene las funciones
1 xt , t y 2 xt , t iguales a uno, y el modelo de incertidumbre esta descrito
únicamente por las componentes dwt y dvt para los estados y las salidas,
respectivamente. Esto indicaría que se modela la diferencia entre la realidad observada
y el modelo como únicamente un proceso de Wiener caracterizados. Simulando que en
los estados no hay conocimiento preciso de la forma en la que los estados del modelo
difieren de la realidad y respecto a la salida o medidas, se considere imprecisión de
sensores o altos ruidos de señales aditivas.
∆
w z
u P0 y
z w P11 P12 w
P0 (4.2)
y u P21 P22 u
1
Fu P0 , : P22 P21 I P11 P12 (4.3)
De modo que
y Fu P0 , u (4.4)
diag 1 ,..., N
1
LTI i es nxn y i i (4.5)
i es LTI y causal
119
Para presentar este enfoque es necesario pensar en las incertidumbres en las variables
de estado, como series de tiempo. Éstos representan cómo se diferencia el modelo del
sistema (planta real) en un tiempo determinado, del cual se disponen datos para la
construcción del modelo de incertidumbre.
Yˆ t Y t 1 Y t 1 Y t 2 e t 1 (4.6)
Mínimos cuadrados lineales versus no lineales: Los modelos ARIMA que incluyen
solamente términos (auto regresivos, AR) son casos especiales de modelos de
regresión lineales, ya que pueden ser aproximados usando mínimos cuadrados
ordinarios. Las predicciones (AR) son una función lineal de los coeficientes así como
una función lineal de los datos pasados.
En principio, las estimaciones de los coeficientes auto regresivos (AR) por mínimos
cuadrados pueden ser calculados exactamente a partir de correlaciones en una sola
iteración. En la práctica, se puede aproximar un modelo (AR) en el procedimiento de
regresión múltiple (solo hacer regresión del diferencial en los retardos en sí mismo).
Pero deberían tenerse resultados levemente diferentes del procedimiento ARIMA.
Los modelos ARIMA que incluyen términos de media móvil (MA) son similares a los
modelos de regresión, pero no pueden ser aproximados de ninguna manera por
mínimos cuadrados ordinarios: las predicciones son una función lineal de los datos
pasados, pero ellos son funciones no lineales de los coeficientes, por ejemplo un
modelo ARIMA sin constante es un promedio móvil exponencialmente ponderado en
el cual las predicciones son funciones no lineales de los términos de media móvil MA
llamados θ.
Para identificar correctamente el modelo ARIMA para una serie de datos, en nuestro
caso, una componente del vector de estados, es necesario comenzar por identificar el
orden o los órdenes de diferenciación necesarios para estacionarizar las series y
remover las características de algún tipo de etapas, en conjunción con la
transformación de estabilización de la varianza como el llamado “logging” o el
“deflating”.
El uso de los modelos de datos muestreados radica en la relación que existe entre los
modelos de tiempo discreto que describen las muestras tomadas del sistema y los
modelos de tiempo continuo originales. Lo que se suele hacer es simplemente
muestrear rápidamente y después reemplazar las derivadas de los modelos de tiempo
122
Sin embargo, la situación es más compleja que en los sistemas lineales. De hecho,
cuando Yuz hizo su investigación de los modelos de datos muestreados lineales,
encontró que la caracterización explícita de las dinámicas de muestreo de ceros no
había sido resuelta. No obstante, en el trabajo de Mónaco y Normand-Cryrot (Monaco
y Normand-Cryrot, 1988) ya se había reportado una caracterización de el muestreo no
lineal de ceros. Por otro lado, en el trabajo de Kazantzis y Kravaris (Kazantzis y
Kravaris, 1997) se estudiaron las propiedades de los modelos de datos muestreados
para sistemas no lineales.
Así pues, se presentarán los modelos de datos muestreados para representar sistemas
no lineales estocásticos. Debe afirmarse que encontrar la solución exacta de una
ecuación diferencial (estocástica) es usualmente imposible. Como consecuencia,
serán obtenidos modelos de datos muestreados, sin embargo la confiabilidad del
modelo propuesto será caracterizada en términos de un orden de convergencia muy
específico.
El siguiente desarrollo fue propuesto por Yuz en su tesis doctoral, luego de una base
matemática y de procesos estocásticos que puede ser encontrada en el Apéndice A de
esta tesis.
dxt f xt , t dt 1 xt , t dwt
(4.7)
dyt g xt , t dt 2 xt , t dvt
dxt
f xt , t 1 xt , t wt
dt (4.8)
yt c xt , t
Se asume que existe una transformación Z = Φ(X), de modo que en las coordenadas
nuevas, el sistema pueda ser expresado como:
dZ t
AZ B * f Z,t 1 Z , t wt
dt (4.9)
yt z1
0 0
In 1
A B (4.10)
0 0
0 0 1
dx1 x2 dt
(4.11)
dxn 1 xn dt
dxn f ( x)dt ( x)dwt
Donde las salida es y = x1. Las funciones f(·) y σ(·) son analíticas. Nótese que si se
expande la ecuación de estados se puede usar la primera expansión de Itô-Taylor:
125
T
donde X 0 x1 (0), x2 (0),..., xn (0) y ( X 0 ) wn es el término que agrupa los demás de
la serie.
T
W w1 w2 wn
n l
wn dw (4.13)
0
n l !
l 1,..., n
Y se obtiene
2
xn 1 ( ) xn 1 (0) xn ( ) d xn 1 (0) xn (0) f (X0) ( X 0 ) wn 1 (4.14)
0
2
wl I (4.15)
1,0,...,0
n l ceros
n 1 n
x1 ( ) x1 (0) x2 (0) ... xn (0) f (X0) ( X 0 ) w1 (4.16)
n 1! n!
X X0
X AX B f (X0 ) ( X 0 )W (4.17)
Donde
126
n 1
n 1
x1 0 1 n!
n 1! n 2
x2 1
X A 0 B W W (4.18)
n 1!
1
xn
0 0 0
1
X k AX k B f( k ) (k )Wk (4.19)
1
Donde A y B están definidas anteriormente y Wk Wk corresponde a
T
W w1 w2 wn
n l (4.20)
wn dw ; l 1,..., n
0
n l !
2 n
E y0 y 0 C0 para alguna constante C0. (4.21)
Capítulo 5
IMPLEMENTACIÓN DE LA TÉCNICA DE
FILTRADO NO LINEAL EN SISTEMAS
DINÁMICOS NO LINEALES NO
GAUSSIANOS. ESTIMACIÓN DE ESTADOS,
IDENTIFICACIÓN DE SISTEMAS
“Bayesian filtering based on particle algorithms have recently emerged as an active
research field because, for many application areas it is becoming important to
consider non linearity and non Gaussianity. Bioprocesses are examples of such non
linear and non Gaussian systems. However, to my knowledge no use of particle
algorithms on biotechnological processes has been reported so far. Consequently,
the paper is novel in the sense that proposes to use for states estimation the above
mentioned algorithms”.
Motivación
12
Debe recordarse que la teoría bayesiana está formulada en un dominio discreto.
130
Otra aplicación sobre la cual se hace uso de las técnicas de filtrado no lineal, es un
proceso de fermentación por lotes, para la producción de δ-endotoxinas de Bacillus
131
thuringiensis (Bt). Esta es una bacteria formadora de esporas, la cual durante el proceso
de esporulación produce uno o más cuerpos cristalinos de naturaleza proteica
adyacentes a la espora conocidos como δ-endotoxinas. Posee dos etapas en su ciclo de
vida, una donde sólo presenta crecimiento y una segunda etapa de esporulación. Al
final de crecimiento vegetativo, el agotamiento del medio induce el inicio de la
esporulación. Normalmente la esporulación está acompañada por la síntesis de δ-
endotoxina. Después de la esporulación se completa el proceso con la ruptura de la
pared celular y la liberación de las esporas y los cristales al medio, etapa que se conoce
como lisis celular.
Existen muchos aspectos del problema de estimación en sistemas por lotes (batch),
que son relevantes a nuestra discusión. Por ejemplo, el hecho que los procesos por
lotes sean procesos de transición sin estados estacionarios, bien definidos significa que
es muy frecuente el uso de una técnica no lineal de estimación. Como resultado,
comúnmente se afirma que no es válido asumir que las distribuciones de las variables
son Gaussianas. Adicionalmente, la complejidad del modelo usado para el cálculo en
línea debe mantenerse en un nivel razonable, ya que debe existir una correspondencia
estructural suficiente entre el modelo y la planta.
Un punto muy importante a ser resaltado es que los procesos batch operan de forma
muy dependiente a sus condiciones iniciales y su desempeño cronológico esta
determinado en mayor parte a ellas. Como consecuencia, la característica que
diferencia la estimación de estados aplicada a los procesos por lotes de la estimación
de estados en los procesos continuos, como el presentado anteriormente, es el efecto
de las condiciones iniciales.
132
5.1.1 Justificación
Actualmente el principal proceso biológico que puede ser usado para producción de
líquidos que transporten energía, son las fermentaciones para la producción de etanol,
y la mezcla de acetona, butanol y etanol (ABE). Otro líquido de interés es el biodiesel.
El biodiesel, es un combustible renovable para motores diesel (diesel no es un
combustible exclusivo para motores, también sirve para turbinas, calentadores,
hervidores, etc), que se deriva de los aceites vegetales como la soja. Esos aceites se
convierten en éteres de metanol o de etanol. El biodiesel no es el tópico de esta parte
del trabajo dado que no se produce a partir de la fermentación (Rogers et al., 2007).
Zymomonas mobilis es una bacteria Gram negativa que atrae la atención de los
investigadores debido a que usa la trayectoria de Entner-Doudoroff para producir
energía a partir del catabolismo de glucosa. Esta trayectoria es ineficiente al producir
una molécula de ATP (trisfosfato de adesosina), por molécula de glucosa consumida.
Esta ineficiencia se compensa por la capacidad de la bacteria en metabolizar glucosa a
altas proporciones (Parker et al., 1997). El interés industrial en el uso de la
Zymomonas mobilis, es su capacidad de producción de etanol y sorbitol (Oliveira et
al., 2005). Para la industria de producción, esto es un requisito fundamental, en la
búsqueda de microorganismos etanologénicos más competitivos.
X, S, F2 X, S, P, F2 P
X, S, F6
X, S, F4,
R
Desde el punto de vista del control, se cuenta con un sistema altamente no lineal con 5
variables de estado, planteado como un sistema de 5 ecuaciones diferenciales de las
siguientes variables: Concentración de Biomasa (X), Sustrato (S), Producto (P),
tiempo de retardo del efecto de inhibición que es modelado como el efecto de la media
ponderada de la tasa de cambio de la concentración de etanol (Z) y la variable
intermedia auxiliar para la determinación del efecto de inhibición (I). Adicionalmente,
el modelo fenomenológico se completa con 3 ecuaciones algebraicas en términos de
136
dX F4 F6 (3.1)
R DS 4 1 X
dt V F2
dS 1 (3.2)
Qp X DSin DS
dt Yp / s
dP (3.3)
Qp X DP
dt
dZ (3.4)
I Z
dt
dI (3.5)
Qp X DP I
dt
Z Z
e e
f 1/ 2 1 Z Z
(3.6)
e e
a b
P P Pob
max S 1 1
Pma Pmb Pob (3.7)
e
S S Si
Ks S
K i Si
Donde Pma y Pob son factores de inhibición de etanol para la tasa específica de
crecimiento expresados en (g/l), Pmb el factor de máxima inhibición de etanol para el
crecimiento de las células en (g/l). Sin es la concentración de sustrato a la entrada,
Yp/s es el rendimiento producto sustrato, a y b son variables para el modelado de la
tasa de crecimiento estática e , Ks Ki y Si son constantes usadas para el mismo
propósito; max es la tasa de crecimiento máxima. Con las condiciones de:
137
P Pob
0, P Pob
Pmb Pob
(3.8)
S Si 0, S Si
P Pob
1, P Pmb
Pmb Pob
f * e
(3.9)
Qp Qpmax
S
1
P (3.10)
K mp S Pme
Donde:
En este tipo de proceso, es importante lograr mediante filtrado una buena estimación
de los estados no medibles del sistema como Z e I . También es importante lograr una
buena estimación de la concentración de microorganismos en el bioreactor, puesto que
su medición es muy costosa y es posible que no pueda ser realizada en línea. Los
experimentos de (Bravo et al, 2000) han mostrado también que la síntesis de etanol
está asociada con su tasa de crecimiento. Debido a esto, es necesario conocer los
estados del sistema para diseñar posteriormente una estrategia de alimentación
138
adecuada para mantener el cultivo a una alta tasa constante que pueda mejorar la
productividad de etanol.
Para poder evidenciar mejor las interacciones entre las variables definidas
previamente, se puede escribir el sistema anterior en términos de las variables de
estado. De las ecuaciones (5.1) a (5.6) se tiene que:
F4 F6
x1 R DS 4 1 x1
V F2
1
x2 Qp x1 D * Sin Dx2
Yp / s
x3 Qp x1 Dx3 (3.11)
x4 x5 x4
x5 Q p x1 Dx3 x5
a b
x3 x3 Pob
max x2 1 1
1 e x4
e x4 Pma Pmb Pob F4 F6
x1 1 x4 x4
R DS 4 1 x1
2 e e x2 x2 Si V F2
Ks x2
K i Si
1 x2 x3
x2 Q pmax 1 x1 D * Sin Dx2
Yp / s K mp x2 Pme
(3.12)
x2 x3
x3 Q pmax 1 x1 Dx3
K mp x2 Pme
x4 x5 x4
x2 x3
x5 Q pmax 1 x1 Dx3 x5
K mp x2 Pme
a b
x3 x3 Pob
max x2 1 1
1 e x4
e x4 Pma Pmb Pob
x1 1 x4 x4
x1
2 e e x2 x2 Si
Ks x2
K i Si
RD / 4 Ds R 1 x1
1 x2 x3
x2 Q pmax 1 x1 D * Sin Dx2 (3.13)
Yp / s K mp x2 Pme
x2 x3
x3 Q pmax 1 x1 Dx3
K mp x2 Pme
x4 x5 x4
x2 x3
x5 Q pmax 1 x1 Dx3 x5
K mp x2 Pme
140
Parámetros Oscilatorio
8.77
0.0366
0.824
21.05
a 0.3142
b 1.415
Qpmax 2.613
max
0.41
Yp / s 0.495
Tabla 5.1. Parámetros Comportamiento Oscilatorio sostenido
inclusión como parte de los algoritmos de estimación y control que serán descritos en
el Capítulo 6. Los tiempos de muestreo considerados fueron desde 0.05 a 0.1 hora es
decir, entre 3 y 6 minutos. Estos tiempos fueron elegidos ya que se consideran
apropiados para este tipo de sistema dinámico con base en un análisis de las dinámicas
involucradas. Debe anotarse también que normalmente las fermentaciones de tiempos
de duración largas como 150 horas, son muestreadas cada hora ya que los métodos de
medida son fuera de línea.
u t MODELO DE PROCESO
y t
dxt f xt , t dt 1 xt , t dwt
dyt g xt , t dt 2 xt , t dvt
xt 1 f xt , t dt 1 xt , t dw,t
yt g xt , t dt 2 xt , t dv ,t
xˆt
Figura 5.2: Esquema de Estimación en simulación para el proceso de fermentación
de Z.m
Las señales del vector de entrada u(t) de la Figura 5.2 corresponden a la concentración
de sustrato a la entrada S in , la tasa de dilución D y el término de recirculación de
microorganismos R ; y las señales del vector de salida y(t) hacen referencia a la
142
Entonces, es importante recalcar que las mediciones usadas para alimentar el filtro,
como información externa, son las tomadas del modelo que corresponde en este caso a
la “planta real” y están ya contaminadas con la incertidumbre asociada 1 xt , t dwt
(porque el sustrato y el producto son también estados del sistema). A estas mediciones
(sustrato y el producto) se les adiciona su término correspondiente de incertidumbre
de medida 2 xt , t 1 y dvt , mediante el ruido gaussiano y su término de difusión
para simular los efectos del ruido electrónico en sensores, que fue discutido
anteriormente.
Experimento 1:
La estructura del filtro tiene como planta real el modelo de (Daugulis et al, 1997)
validado en la tesis de Li, 1999, con los parámetros dinámicos encontrados y
reportados por Raposso et al., 2005. Presentados en las Tablas 5.1 y 5.2. Como
modelo del filtro, se uso el mismo modelo con ruidos aditivos descritos a
continuación.
El modelo utilizado para el filtrado de partículas es del tipo tipo 3 de Schon (Schon,
2005) definido como: Espacio de estados no linear discreto con ruido aditivo
(Nonlinear discrete state-space model with additive noise) de la forma
x t 1 f x t w t
y t h x t v t
(5.14)
2
estado1 0 0 0 0 0.0625 0 0 0 0
2
0 estado 2 0 0 0 0 0.04 0 0 0
2
R1 0 0 estado 3 0 0 0 0 .04 0 0
2
0 0 0 estado 4 0 0 0 0 .01 0
2
0 0 0 0 estado 5
0 0 0 0 .01
2
salida1 0 0.0625 0
R2 2
0 salida 2
0 0.0625
144
En las Figura 5.3 se encontrará que el filtro sigue bien el modelo, que estando en
concurrencia con los parámetros de la fermentación experimental de (Raposso et al.,
2005), se convierte en otro modelo diferente al de (Daugulis et al., 1997), (Li, 1999).
Estas fermentaciones en continuo tienen esas características cuando se les cambian las
condiciones iniciales, y se reflejan en estos parámetros de modelo. Los datos
marcados con asteriscos en la Figura 5.3.a son los datos reales de esas
experimentaciones. Con estas difusiones a medida que aumenta el número de
partículas, disminuye el error al comienzo de la estimación. Nótese de la Figura 5.3a.
que hay una diferencia entre el modelo del filtro (Filtro) y el usado para simular el
bioproceso (Modelo), debido a la diferencia de los valores de ruido aditivo del modelo
de Schon w t y v t para el filtro en contraste con los de 1 xt , t , dwt , 2 xt , t 1 y
dvt con los que se simuló la planta real.
2.5
1.5
1
Filtro
0.5
Modelo
Dato Real como referencia
0
0 10 20 30 40 50
Tiempo [Horas]
Figura. 5.3. a. Estimación de biomasa con 50 partículas y esquema de re muestreo
Multinomial
145
Figura. 5.3. b. Errores en la estimación de los estados del bioreactor, durante una
prueba en simulación de 150 horas con un tiempo de muestreo de 0.1 horas.
La importancia de las Figura 5.4 radica en el hecho que las variables biotecnológicas
no son siempre muy fáciles de muestrear y ésta aproximación entrega un valor de
estimación suficientemente aproximado, que puede llegar a reemplazar la medida real
en un intervalo de tiempo teniendo en cuenta el modelo de ruido y difusión de gran
magnitud agregado. De la misma manera, el desempeño del filtro nos permite asumir
las incertidumbres en las medidas y perturbaciones en los mismo métodos de
medición.
0.5
0.5
0
Filtro SIR
Simulado como Real
-0.5
-1
10 20 30 40 50 60 70 80 90
Tiempo Horas
Figura. 5.4b. Concentración de biomasa estimada por el filtro SIR. La línea punteada
describe la estimación mientras que la línea continua describe el valor del modelo
considerado como real.
Tabla 5.4. Errores obtenidos en diferentes pruebas, ilustración de los peores casos
En la Figura 5.4a 5.4b, puede observarse el resultado de una prueba realizada para 500
partículas, usando el esquema de re muestreo determinístico para el bioproceso
continuo. Se consideró como salida medible del sistema el Sustrato y el Producto (S y
P). Para el modelo de ruido aditivo a las mediciones, se utilizó ruido blanco normal,
con media y varianza de adecuada y considerable magnitud.
En las Figura 5.5 a y b pueden observarse las estimaciones de las variables Z y I que
corresponden a las variables de inhibición. Cabe aclarar que este par de variables
hacen parte del vector de estados, pero no son medibles físicamente. Lo que se
persigue es lograr una buena estimación de dichas variables a partir de los datos de
Sustrato y Producto, ya que en caso de un implementación en línea, la buena
estimación de estas variables hará mucho más viable su monitoreo con fines de
optimización e inclusive de producción.
1
Variable de Inhibición Z
Filtro SIR
0.5 Real
-0.5
0.6
0.4
Filtro SIR
0.2 Real
0
-0.2
-0.4
-0.6
-0.8
-1
70 80 90 100 110 120 130 140 150
Tiempo [Horas]
Figura. 5.5b. Dinámica de inhibición I estimada por el filtro SIR. Línea punteada
representa las estimaciones y las líneas continuas son consideradas las dinámicas
reales del sistema.
1.5
Error
1
Porcentaje %
0.5
-0.5
-1 Biomasa
Inhibición Z
-1.5 Inhibición W
-2
0 50 Tiempo [Horas] 100 150
Biomasa
0.5
-0.5
Filtro SIR
Real
-1
134 136 138 140 142 144 146 148 150
Tiempo [Horas]
Figura. 5.7b Estimación de Biomasa
0.3
0.25
0.2
0.15
53 53.2 53.4 53.6 53.8 54 54.2
Tiempo [Horas]
Por esta razón, solamente se lograrían aproximaciones. Debe recalcarse que en este
caso particular no estamos interesados en la solución exacta de la ecuación diferencial
estocástica, pero si estamos interesados en las muestras que puedan ser tomadas del
modelo para calcular la solución aproximada de la función de densidad de
probabilidad de los estados (pdf) a través de simulaciones recursivas de Monte Carlo,
para finalmente estimar los estados.
El modelo de incertidumbre para los estados, usado en el primer caso para aproximar
los datos reales consiste de un proceso estocástico definido por las siguientes
componentes:
153
z (t ) p ( z (t ) z (t 1)), z (t ) U a, b
1 xt , t p ( x(t ) z (t t 1), x(t 1)) (5.14)
,
1 xt , t p( x(t ) z (t t 1), x(t 1))dw t
ut yt
DATOS REALES
PROCESADOS
xt 1 f xt , t dt 1 xt , t dw,t
yt g xt , t dt xt , t dv ,t
2
xt
Figura. 5.8 Esquema de estimación usado para la prueba con datos reales
Para la configuración se usa el modelo dentro del filtro y se toman datos reales de
sustrato y producto, de los datos de las fermentaciones reportadas por (Raposso et al,
2005). Para el modelo interno del filtro se le agregaron modelos de perturbación y de
incertidumbre al modelo base bajo las consideraciones de incertidumbres de modelado
y medición de diversa naturaleza. El modelo se convierte de esta forma en un modelo
de ecuaciones diferenciales estocásticas.
La Figura 5.8 presenta el esquema de estimación para la primera prueba con datos
reales. Las señales de entrada u(t) (concentración de entrada de sustrato Sin , tasa de
dilución D y el término de recirculación de microorganismos R ), las señales de
salida y(t) (salida de sustrato S y de producto P ) alimentan el filtro que hace la
estimación de los estados. La simulación se hace en condiciones de lazo abierto
sostenido oscilatorio, con entradas constantes determinadas en la experimentación.
Los resultados se comparan con los datos reales. A partir del modelo patrón que
considera las incertidumbres, se observaron las dinámicas de inhibición, que hacen
que el proceso sea altamente no lineal.
154
Los datos reales reportados por (Raposso et al, 2005) se han muestreado a la velocidad
de 2-3 muestras por hora. Para propósitos de filtrado, los datos fueron preprocesados
y el vector de datos original fue completado usando una herramienta de interpolación
de Matlab®. Es importante resaltar que los datos de sustrato y producto son mucho
más factibles de medir, que los datos de concentración de biomasa que se obtienen
generalmente usando métodos fuera de línea. Los valores de sustrato y producto
usados como señal de entrada del filtro para los experimentos en estado oscilatorio se
pueden ver en la Figura 5.9
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90
a. Tiempo [Horas]
Datos Reales de Sustrato
Concentración de Glucosa [g/L]
120
100
Interpolación Datos Reales
80 Datos Reales
60
40
20
0 10 20 30 40 50 60 70 80 90
b. Tiempo [Horas]
155
110
100
70
60
50
40
0 10 20 30 40 50 60 70 80 90
c. Tiempo [Horas]
Figura 5.9. Interpolación de datos reales. a. Concentración de Biomasa real b.
Concentración de sustrato. c. Concentración de producto.
Parámetro Valor
α 0.1000*100
β 0.430457*0.1
δ 0.10000*0.01
λ 0.241458*10
a 0.1000000*10
b 1.72620*10
Qpmax 0.3828880*10
Yps 0.514353
µmax 0.41
Tabla 5.5. Parámetros usados para estimación con datos reales de Zymomonas
mobilis.
Para aproximar los datos reales mediante el filtrado de partículas se usaron un número
de partículas desde 10 hasta 1000, usando tres tipos de re muestreo diferentes:
residual, determinístico y multinomial, mencionados en el Capítulo 3. Se realizaron
pruebas haciendo no solo variaciones en el número de partículas sino también en las
componentes aleatorias de los modelos de incertidumbre.
Otro de los objetivos de hacer pruebas de varias mediciones era probar la efectividad
del filtro, respecto a su costo computacional, analizando la eficiencia de los esquemas
de re muestreo utilizados en esta configuración del filtro de partículas. Para ello, se
156
conFiguraron cinco números de partículas con las cuales hacer la prueba, en tres
esquemas de re muestreo mencionados anteriormente. Se fijaron las componentes de
difusión y no se hicieron las demás variaciones en las componentes aleatorias, ya que
el objetivo era medir el tiempo de cálculo. El procesador principal de prueba fue un
Intel Core duo 3.0 GHz 1GB RAM, cuyos resultados en tiempo de simulación se
reportan en la Tabla 5.6. También se probaron los mismos algoritmos en un
procesador AMD Sempron 2200 192 RAM; los resultados obtenidos fueron cercanos
a 2 o tres veces el tiempo de cálculo que presentó el primer procesador.
Posteriormente se realizaron test en un procesador Intel Core 2 duo T5600@ 1.83 Ghz
2GB RAM reduciendo el tiempo de cálculo del primero a aproximadamente la mitad.
posible que la estimación generada sea adecuada, y puede ser usada para propósitos
reales de control.
Resultados:
Resultados:
3
Estimación de células de Z.m
Concentración Biomasa [g/L]
Filtro
2.5 Interpolación de Datos Reales
Datos Reales
2
1.5
0.5
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
Figura. 5.10 a. Estimación de Biomasa experimento 5 usando 10 partículas y
re muestreo residual
159
3.5
Filtro
3 Interpolación de Datos Reales
Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90 100
Tiempo [Horas]
Figura. 5.10 b. Estimación de Biomasa experimento 5 usando 10 partículas y
re muestreo multinomial
3
Estimación de células de Z.m
Concentración Biomasa [g/L]
Filtro
2.5 Interpolación de Datos Reales
Datos Reales
2
1.5
0.5
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
Figura. 5.10 c. Estimación de Biomasa experimento 5 usando 10 partículas y
re muestreo determinístico
Estimación de células de Z.m
Concentración Biomasa [g/L]
3.5
Filtro
3 Interpolación de Datos Reales
Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
Figura. 5.10 d. Estimación de Biomasa experimento 5 usando 100 partículas
y re muestreo multinomial
Figura. 5.10. a. b. c. d. Experimento 5. La línea discontinua describe la concentración
de biomasa estimada por el filtro SIR y la línea punteada el valor de la interpolación
160
de los datos reales. Los datos reales se muestran como asteriscos. Esta prueba contiene
desviación de la media de los términos de difusión.
'
1 0, 0, 0, 0, 0
2
estado1 0 0 0 0 0.04 0 0 0 0
2
0 estado 2 0 0 0 0 0.04 0 0 0
2
R1 0 0 estado 3 0 0 0 0 0.04 0 0
2
0 0 0 estado 4 0 0 0 0 0.01 0
2
0 0 0 0 estado 5
0 0 0 0 0.01
'
2 0, 0
2
salida1 0 0.0625 0
R2 2
0 salida 2
0 0.0625
161
En la Tabla 5.8 se presentan los tiempos de cómputo del algoritmo (Intel Core 2 duo
T5600@ 1.83 Ghz 2GB RAM), para distintos números de partículas y esquemas de re
muestreo, ilustrando los casos de rendimiento del algoritmo. Igualmente se usaron dos
tipos de series de datos reales para las pruebas de estimación por filtro.
3
Filtro
Interpolación de Datos Reales
2.5 Datos Reales
1.5
10 20 30 40 50 60 70 80
Tiempo [Horas]
a. Estimación de Biomasa experimento 7 usando 10 partículas y re muestreo
residual
Estimación de células de Z.m
Concentración Biomasa [g/L]
3.5 Filtro
Interpolación de Datos Reales
3 Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90 100 110
Tiempo [Horas]
b. Estimación de Biomasa experimento 7 usando 10 partículas y re muestreo
multinomial
Estimación de células de Z.m
Concentración Biomasa [g/L]
2.5
1.5
Filtros
1 Interpolación de Datos Reales
Datos Reales
10 20 30 40 50 60 70 80
Tiempo [Horas]
c. Estimación de Biomasa experimento 7 usando 10 partículas y re muestreo
determinístico
163
3
Filtro
Interpolacion de Datos Reales
2.5 Datos Reales
1.5
10 20 30 40 50 60 70 80
Tiempo [Horas]
d. Estimación de Biomasa experimento 7 usando 100 partículas y re muestreo
determinístico
Resultados:
3.5 Filtro
Interpolación de Datos Reales
3 Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90 100
Tiempo [Horas]
a. Estimación de Biomasa experimento 9 usando 10 partículas y re muestreo
multinomial
Estimación de células de Z.m
Concentración de Biomasa [g/L]
3.5
Filtros
3 Interpolación de Datos Reales
Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
b. Estimación de Biomasa experimento 9 usando 50 partículas y re muestreo
determinístico
Estimación de células de Z.m
Concentración de Biomasa [g/L]
3
Filtro
Interpolación de Datos Reales
2.5 Datos Reales
1.5
10 20 30 40 50 60 70 80
Tiempo [Horas]
c.
Estimación de Biomasa experimento 9 usando 100 partículas y re muestreo
multinomial
Figura. 5.12. a. b. c. Experimento 9. La línea discontinua describe la concentración de
biomasa estimada por el filtro SIR y la línea punteada el valor de la interpolación de
los datos reales. Los datos reales se muestran como asteriscos.
166
'
1 0, 0, 0, 0, 0
2
estado1 0 0 0 0 1 0 0 0 0
2
0 estado 2 0 0 0 0 1 0 0 0
2
R1 0 0 estado 3 0 0 1.0e 4 0 0 1 0 0
2
0 0 0 estado 4 0 0 0 0 1 0
2
0 0 0 0 estado 5
0 0 0 0 1
'
2 0, 0
2
salida1 0 0.0625 0
R2 2
0 salida 2
0 0.0625
multivariable de media 2
2
y matriz de covarianza R2 2 2
”
La Tabla 5.10 resume los tiempos de cómputo del algoritmo, para diferente número de
partículas y esquemas de re muestreo. Según la tabla, puede verse que se demora,
aproximadamente entre segundos (10 partículas) a segundos (1000 partículas), por
cada tiempo de muestreo de 0.1 horas en simulación. Las pruebas fueron realizadas
con un procesador inferior a las de los experimentos 1-10, las características del
sistema son: 192 MB de RAM, AMD sempron 2200+. En las pruebas marcadas con
no, se presentaron interrupcoines en el algoritmo dado que no se terminó la prueba por
limitaciones de hardware. Consecuente con esto, los tiempos de cómputo del
algoritmo aumentaron:
3.5 Filtro
Interpolación de Datos Reales
3 Datos Reales
2.5
1.5
0 10 20 30 40 50 60 70 80
Tiempo [Horas]
a. Estimación de Biomasa experimento 11 usando 10 partículas y re muestreo
residual
168
3.5 Filtro
Interpolación de Datos Reales
3
Datos Reales
2.5
1.5
10 20 30 40 50 60 70 80
Tiempo [Horas]
b. Estimación de Biomasa experimento 11 usando 10 partículas y re muestreo
determinístico
Estimación de células de Z.m
Concentración de Biomasa [g/L]
3.5 Filtro
Interpolación de Datos Reales
3
Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90 100
Tiempo [Horas]
c. Estimación de Biomasa experimento 11 usando 10 partículas y re muestreo
multinomial.
Estimación de células de Z.m
Concentración de Biomasa [g/L]
3.5 Filtro
Interpolación de Datos Reales
3 Datos Reales
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
d. Estimación de Biomasa experimento 11 usando 50 partículas y re muestreo
determinístico
Figura. 5.13. a. b. c. d. Experimento 11. La línea discontinua describe la
concentración de biomasa estimada por el filtro SIR y la línea punteada el valor de la
interpolación de los datos reales. Los datos reales se muestran como asteriscos.
169
'
1 0, 0.5, 0.5, 0, 0
2
estado1 0 0 0 0 0.04 0 0 0 0
2
0 estado 2 0 0 0 0 0.25 0 0 0
2
R1 0 0 estado 3 0 0 0 0 0.25 0 0
2
0 0 0 estado 4 0 0 0 0 0.01 0
2
0 0 0 0 estado 5
0 0 0 0 0.01
'
2 0, 0
2
salida1 0 0.0625 0
R2 2
0 salida 2
0 0.0625
multivariable de media 2
2
y matriz de covarianza R2 2 2
”
Las condiciones sobre las cuales se diseño este experimento, permiten predecir que la
estimación no será satisfactoria. Este experimento se realizó con el fin analizar la
sensibilidad del filtro a las características de la componente estocástica. En la Figura
5.14 pueden verse los resultados.
170
Estimación de Z.m
Concentración de Biomasa [g/L]
3
Filtro
2.5 Interpolación Datos reales
Datos Reales
2
1.5
0.5
0
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
Figura. 5.14. Estimación de Biomasa, las líneas continuas son las de la interpolación
de los datos reales, las líneas punteadas corresponden a la estimación del filtro. Los
modelos de perturbación tienen media cero en el término de difusión y una potencia
alta de ruido correspondiente a la prueba 3.
Observaciones generales:
Las Figuras 5.15 muestran los resultados de la prueba 4 en la cual se usó un ruido de
baja potencia como perturbación de las medidas. Es importante recalcar que los filtros
de partículas pueden no ser aplicados a sistemas con perturbaciones de bajo ruido. En
comparación con los demás métodos de estimación reportados para este proceso
fermentativo, los filtros implementados satisficieron las expectativas y se encontró que
fueron superiores, inclusive en lazo abierto bajo las condiciones más altas de
oscilación.
Estimación de Z.m
Concentración de Biomasa [g/L]
3.5 Filtro
Dato Real
3
2.5
1.5
0.5
0 10 20 30 40 50 60 70 80 90 100 110
Tiempo [Horas]
a. Concentración de Biomasa estimada en la prueba 4 con baja potencia en la
componente de ruido.
Inhibición Z
0.5
Filtro SIR
0.4 Real
0.3
0.2
0.1
10 20 30 40 50 60 70 80 90 100 110
Tiempo [Horas]
b. Estimación de la variable de Inhibición Z
172
Inhibición I
1.5
Filtro SIR
1 Real
0.5
-0.5
-1
0 10 20 30 40 50 60 70 80 90 100 110
Tiempo [Horas]
c. Estimación de la variable de Inhibición I
Figura 5.15. Resultados de la prueba 4 en la cual se usó un ruido de baja potencia
como perturbación de las medidas
Estimación de Producto
Concentración de Etanol [g/L]
110
100 Filtro
Interpolación de Datos reales
90 Datos reales
80
70
60
50
40
0 10 20 30 40 50 60 70 80 90
Tiempo [Horas]
a. Estimación de producto, en línea continua los datos de la interpolación de los
datos reales, en punteado las estimaciones del filtro
Estimación de Sustrato
Concentración de Glucosa[g/L]
120
Filtro
100 Interpolación de Datos reales
Datos reales
80
60
40
20
0
0 10 20 30 40 50 60 70 80 90 100 110
Tiempo [Horas]
b. Estimación de sustrato, en línea continua los datos de la interpolación de los
datos reales, en punteado las estimaciones del filtro
Figura. 5.16. Estimación de entradas. a) Estimación de producto, en línea continua los
datos de la interpolación de los datos reales, en punteado las estimaciones del filtro b)
Estimación de sustrato, en línea continua los datos de la interpolación de los datos
reales, en punteado las estimaciones del filtro.
174
Por su parte la planta consiste en un reactor de tanque agitado con volumen nominal
de 20 litros. Las fermentaciones se desarrollaron en 11 litros de medio de cultivo
(volumen efectivo) y se inocularon al 10% (v/v) con el cultivo de Bt. Mediante lazos
de control tipo ON/OFF se mantiene la temperatura alrededor de 30ºC y el pH entre
6.5– 8.5. Se evita la formación de espuma con el agregado de antiespumante. El flujo
de aire se fija en 22 (litros.min-1). Los valores de porcentaje de OD, pH y temperatura
(T) se registran durante toda la fermentación. La fermentación se suspende cuando
presenta 90% o más de lisis celular.
dV
Fin F fresh.means FR Fout (3.14)
dt
dX v FR Fin F fresh.means FR
X v, R Xv ks ke X v (3.15)
dt V V
dX s FR Fin F fresh.means FR
X s,R Xs ks X v kl X s (3.16)
dt V V
dX FR Fin F fresh.means FR
XR X Xv ke X v kl X s (3.17)
dt V V
Sp
max
(3.18)
Ks Sp
En (Amicarelli et al, 2005) se hizo el modelado de las dinámicas del proceso desde el
punto de vista fenomenológico, en este trabajo se parte del balance general de OD en
el medio de cultivo. La ley de conservación puede aplicarse a la masa total del sistema
o a la de cualquier componente individual que pertenezca a éste.
d (V COD ) (3.20)
Fin CODin Fout COD V rg V rc
dt
dV (3.21)
Fin Fout
dt
b. Proceso por lotes Alimentado (fed - batch): el caudal de salida es nulo Fout=0
por lo que V aumentará en el tiempo en función del caudal de entrada
dV/dt=Fin. La ecuación (5.21) queda:
d (V C OD ) (3.22)
Fin C ODin V (rg rc )
dt
V Permanece dentro del operador diferencial porque varía con el tiempo. El proceso
por lotes alimentado, a diferencia del caso anterior, tiene duración limitada ya que el
volumen no puede incrementarse más allá del volumen útil que posee el bio reactor.
c. Por lotes (batch): ambos caudales son nulos Fin=Fout=0 por lo que los dos
primeros términos de la ecuación (5.21) desaparecen y dV/dt=0.
d (COD ) (3.23)
(rg rc )
dt
La duración del cultivo por lotes es también limitada y depende de las condiciones
iniciales del cultivo. Una vez inoculado el medio, la concentración de biomasa
aumenta a expensas de los nutrientes y cuando el sustrato que limita el crecimiento se
agota, finaliza el lote.
Sólo una parte del caudal de aire (N2 y O2 pero el de interés es el Oxígeno),
introducido en los reactores se disuelve en el medio, dependiendo de la eficiencia de
los elementos de aireación y de la proximidad a la concentración de saturación del O2
en el medio líquido (COD* ) . Esto se evidencia en el modelo de la forma:
177
dXv (3.25)
rc K1 K2 X v
dt
Continuo dCOD K1 dX v K2
V F (CODin COD ) V [ K3 QA (COD * COD )] [ ( ) Xv]
dt Yx / O2 dt Yx / O2
FedBatch d (V C OD ) K dX v K2
Fin C ODin V [ K 3 Q A (C OD * C OD )] [( 1 ) Xv]
dt Yx / O2 dt Y x / O2
Batch d (COD ) K dX v
[( K 3 QA (COD * COD ))] [ 1 K2 X v ]
dt Yx / O2 dt
Tabla 5.11. Dinámica del Oxígeno Disuelto
Consideraciones:
Desde el punto de vista estocástico (Rossi y Vila, 2005) propusieron una formulación
como un problema de filtrado. Los términos estocásticos se introducen en las
dinámicas, en las medidas y en las condiciones iniciales. El objeto de interés es la ley
de probabilidad condicional del estado completo dadas unas medidas ruidosas de
algunos componentes. En los trabajos de (Rossi et al, 2005), (Rossi et al 2006) se
presenta el modelamiento de incertidumbre de las dinámicas, mediante un sistema
diferencial estocástico consistente con la noción de invariante. Luego, diseñaron un
conjunto de observadores asintóticos, que usaron en conjunto con las componentes
observadas, para aproximar la ley de probabilidad de las no observadas. Los autores
mencionados, además establecieron la relevancia del modelado estocástico en la
biotecnología, lo que permite usar el “know how” existente en biología y en
optimización. Además de ello tuvieron éxito mostrando la factibilidad de usar esta
aproximación en simulaciones numéricas.
El uso de observadores de estado en procesos por lotes debe atender sobre todo a la
variabilidad e incertidumbre (estructural o paramétrica) de los mismos. Russell,
asegura que el problema de control de calidad en procesos por lotes y semi continuos
se puede reducir a la estimación de los estados del proceso a partir de las condiciones
iniciales y de las de proceso (Russell et al 2000). En el trabajo de (Dondo y Marqués,
2002) se enumeran algunas características típicas de estos procesos a tener en cuenta.
Se pueden mencionar trabajos en donde se realiza estimación de estados en
bioprocesos como (Quintero et al, 2004, 2007, 2008a), (Silveira y Molina, 2005), (Leal
et al, 1997) entre otros. Dochain, propone un observador de estados para un bioproceso
con cinética de estructura conocida y de parámetros cinéticos desconocidos (Dochain,
2002). A tal fin propone el uso de estos parámetros como parámetros extras de diseño
para garantizar un error de estado estacionario igual a cero. El autor divide a los
observadores en dos clases: los clásicos como el de Luenberger y el de Kalman y los
asintóticos. La primera clase depende completamente de la precisión del modelo, la
segunda supone que la incertidumbre del modelo radica en la cinética de los
bioprocesos. En los trabajos de (Leal et al, 1997), (Leal et al, 2001) se realiza la
estimación de biomasa en un bioproceso usando como observador diferentes redes
neuronales artificiales (RNA) y provocando luego la fusión sensorial. (Dorsey y Lee,
2003) abordaron la predicción de variables en procesos discontinuos empleando
técnicas de identificación por subespacios. De forma general, tanto en procesos
continuos como discontinuos, la incertidumbre en el modelo y los parámetros hace que
surjan los observadores asintóticos los cuales no requieren el conocimiento perfecto del
180
Resultados:
En esta sección, se presentan los resultados obtenidos y las gráficas muestran las
estimaciones de las variables de estado que resumen el rendimiento de los filtros. Se
realizan algunas observaciones acerca del número de partículas usadas, el esquema de
re muestreo, el costo computacional y de los criterios de calificación de la efectividad
de la solución propuesta. Finalmente, se sugieren mejoras posibles de realizar sobre
esta aproximación usando técnicas avanzadas de filtrado no lineal. Estos resultados
presentan una innovación al utilizar las técnicas de filtrado no lineal no gaussiano en
bioprocesos por lotes.
Primera Aproximación:
En la Figura. 5.17 se presenta el diagrama de bloques del estimador. Donde las señales
u(t) y y(t) correspondientes al bloque denominado “Modelo de Proceso”, alimentan el
“Modelo del Filtro” quien realiza la estimación de la biomasa (células esporuladas
mas vegetativas). Estas estimaciones son las que se comparan con las entregadas por el
modelo patrón (considerado el proceso real) (Quintero et al, 2007). Se realizaron
pruebas con diferentes tiempos de muestreo para la fermentación por lotes, de la cual
se consideran aproximadamente las 20 primeras horas.
Puede verse que al igual que en la sección 5.1.3, los términos de difusión tanto del
filtro como del modelo del proceso, fue la expresión de 1 xt , t como un ruido
gaussiano dependiente de una variable interna distribuida uniformemente, actualizada
cada instante de muestreo y partícula a partícula; en combinación de distintos ruidos
gaussianos en dwt y dw,t . Para el valor de la componente de las mediciones, se usó el
valor de 2 xt , t 1 igualmente combinado con muestras de distintos ruidos
gaussianos en las componentes de dvt y dv ,t .
xt 1 f xt , t dt xt , t dw,t
,
yt g xt , t dt dv t
xt
Estimación deBiomasa
10
9
Concentración de Células [g/l]
6
Datos de Modelo Sin incertidumbre
Filtro SIR Residual
5
3 4 5 6 7 8 9
Tiempo [Horas]
11
10
9
Concentración de Células [g/l]
1
2 4 6 8 10 12 14
Tiempo [Horas]
Figura. 5.19. Estimación de Biomasa con modelos en simulación y componentes de
difusión adicionadas. Tiempo de muestreo (t0) de 0.01 horas re muestreo
determinístico
184
De la Figura 5.19 puede observarse que los resultados son satisfactorios y que el
algoritmo converge para este tipo de estimación. Los esquemas de re muestreo tienen
una influencia en el tiempo de simulación, pero respecto a la calidad de los resultados,
sus variaciones no son significativas. Los resultados de las pruebas realizadas con
tiempos de muestreo más cortos (de 0.2 horas), para esquemas de re muestreo
determinístico, arrojaron el menor tiempo de ejecución del algoritmo de filtrado, para
una caracterización del ruido de medición de mayor magnitud. Las estimas logradas
son adecuadas y satisfactorias de acuerdo con la magnitud de los errores de estimación,
frente a los obtenidos con mayores tiempos de muestreo. Los siguientes resultados
mostrados en las Figuras 5.20-5.24. se obtuvieron a partir de una simulación con
tiempos de muestreo de 1 hora, usando 500 partículas
Celulas Vegetativas
0.5
Real
Filtro SIR
-0.5
-1
0 5 10 15 20 25 30
Tiempos de Muestreo
Figura. 5.20. Resultados normalizados de estimación de células vegetativas para la
prueba 7 del filtro SIR con esquema de re muestreo residual. Línea continua el modelo
considerado como real, puntos el resultado del filtro.
185
Celulas Esporuladas
Real
-0.4 Filtro SIR
-0.5
-0.6
-0.7
-0.8
-0.9
-1
5 10 15 20 25 30
Tiempos de Muestreo
Figura. 5.21. Resultados normalizados de estimación de células esporuladas para la
prueba 7 del filtro SIR con esquema de re muestreo residual. Línea continua el modelo
considerado como real, puntos el resultado del filtro.
Sustrato Primario
Real
0.8
Filtro SIR
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-1
5 10 15 20 25 30
Tiempos de Muestreo
Figura. 5.22. Resultados normalizados de sustrato primario para la prueba 7 del filtro
SIR con esquema de re muestreo residual. Línea continua el modelo considerado
como real, en puntos el resultado del filtro
186
Sustrato Secundario
-0.2 Real
Filtro SIR
-0.3
-0.4
-0.5
-0.6
-0.7
-0.8
-0.9
-1
5 10 15 20 25 30
Tiempos de Muestreo
Figura. 5.23. Resultados normalizados de sustrato secundario para la prueba 7 del
filtro SIR con esquema de re muestreo residual. Línea continua el modelo considerado
como real, en puntos el resultado del filtro.
Oxigeno Disuelto
0.6
0.4
0.2
-0.2 Real
Filtro SIR
-0.4
-0.6
-0.8
-1
0 5 10 15 20 25 30
Tiempos de Muestreo
Figura. 5.24. Resultados normalizados de oxígeno disuelto para la prueba 7 del filtro
SIR con esquema de re muestreo residual. Línea continua el modelo considerado
como real, en puntos el resultado del filtro
Comentarios:
5.28se muestran las graficas de las superficies generadas por los errores de las pruebas
con tiempos de muestreo de 0.5 horas y 1 hora, con características de las varianzas de
las mediciones correspondientes a 0.04 y 0.25 usando 100 y 500 partículas:
Celulas Vegetativas
-3
x 10
2.5
2
Error Máximo
1.5
0.5
500
400 0.25
0.2
300 0.15
200 0.1
0.05
Particulas 100 0
Varianzas de Medida
12
Error Máximo
10
0.4
0.3
0.2
0.1 200 150 100
350 300 250
500
0 450 400
Varianzas de Medida Particulas
-3
x 10
1.6
1.4
Error Máximo
1.2
0.8
0.6
500
400 0.25
0.2
300 0.15
200 0.1
0.05
Particulas 100 0
Varianzas de Medida
Figura. 5.27. Superficie generada por el máximo de los errores oxigeno disuelto
usando los 3 tipos de re muestreo a 1 hora de muestreo.
Oxigeno Disuelto
-4
x 10
14
12
Error Máximo
10
600
8
400
6
200 0.25
0.2
0.1 0.15
40 0.05
Varianzas de Medida
Figura. 5.28. Superficie generada por el máximo de los errores oxigeno disuelto
usando los 3 esquemas de re muestreo a 0.5 horas de muestreo
3.5
Celulas Vegetativas
Celulas Esporuladas
3 Sustrato Primario
Sustrato Secundario
Oxgeno Disuelto
2.5
1.5
0.5
0
0 50 100 150
Tiempos de Muestreo
Figura. 5.29. Error en la estimación de los estados para la prueba de el filtro SIR
usando muestreo residual, a 0.2 horas de tiempo de muestreo.
Celulas Vegetativas
0.6
0.4
0.2
-0.2
-0.4
Real
-0.6 Filtro SIR
-0.8
-1
10 20 30 40 50 60 70 80
Tiempos de Muestreo
Figura. 5.30. Resultados normalizados de Células Vegetativas para la prueba del filtro
SIR usando muestreo residual, a 0.2 horas de tiempo de muestreo. Línea continua el
resultado del filtro, en puntos el modelo considerado como real.
190
Celulas Esporuladas
Real
-0.8 Filtro SIR
-0.85
-0.9
-0.95
-1
-1.05
10 15 20 25 30 35
Tiempos de Muestreo
0.8 Real
Filtro SIR
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-1
-1.2
5 10 15 20 25 30 35 40 45 50 55
Tiempos de Muestreo
Figura. 5.32. Sustrato primario para la prueba del filtro SIR usando muestreo residual,
a 0.2 horas de tiempo de muestreo. Línea continua el resultado del filtro, en puntos el
modelo considerado como real.
191
Sustrato Secundario
-0.8 Real
Filtro SIR
-0.82
-0.84
-0.86
-0.88
-0.9
-0.92
-0.94
-0.96
-0.98
30 35 40 45 50 55 60
Tiempos de Muestreo
Figura. 5.33. Resultados normalizados de Sustrato secundario para la prueba del filtro
SIR usando muestreo residual, a 0.2 horas de tiempo de muestreo. Línea continua el
resultado del filtro, en puntos el modelo considerado como real.
Oxigeno Disuelto
0.4 Real
Filtro SIR
0.2
-0.2
-0.4
-0.6
-0.8
-1
10 20 30 40 50 60
Tiempos de Muestreo
Figura. 5.34. Resultados normalizados de Oxígeno disuelto para la prueba del filtro
SIR usando muestreo residual, a 0.2 horas de tiempo de muestreo. Línea continua el
modelo considerado como real, en puntos el resultado del filtro.
Resultados:
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-1
0 500 1000 1500 2000 2500 3000
Tiempos de Muestreo
Figura. 5.35. Celulas vegetativas para la prueba de el filtro SIR usando muestreo
residual, a 0.01 horas de tiempo de muestreo. Línea continua el resultado del filtro, en
puntos el modelo considerado como real.
193
Sustrato Primario
0.8
Real
0.6 Filtro SIR
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-1
-1.2
200 250 300 350 400 450 500 550 600
Tiempos de Muestreo
Figura. 5.36. Sustrato primario para la prueba de el filtro SIR usando muestreo
residual, a 0.01 horas de tiempo de muestreo. Línea continua el resultado del filtro, en
puntos el modelo considerado como real.
Oxigeno Disuelto
0.8
0.6
Real
0.4 Filtro SIR
0.2
-0.2
-0.4
-0.6
-0.8
-1
-1.2
0 500 1000 1500 2000 2500 3000
Tiempos de Muestreo
Figura. 5.37. Oxígeno disuelto para la prueba de el filtro SIR usando muestreo
residual, a 0.01 horas de tiempo de muestreo. Línea continua el resultado del filtro, en
puntos el modelo considerado como real.
Resultados:
Nótese que los resultados no son satisfactorios como se desearía y que es necesario
hacer ajustes para lograr unos resultados más confiables y de mayor precisión, es por
esto que se realizó más cuidadosamente una nueva aproximación, más compleja de
acuerdo a lo planteado en el capítulo 4.
194
Para el uso del filtro de partículas como estimador de Biomasa, fue necesario utilizar
la información muestreada de Sustrato Primario y Oxígeno Disuelto, obtenidas del
experimento de fermentación “A” de 18 horas de duración y del experimento “B” de 8
horas de duración. Los datos de Oxígeno disuelto fueron tomados con una placa de
adquisición de datos cada 0.01 horas. El periodo de muestreo seleccionado
previamente como el más representativo en este tipo de sistema dinámico, es 6
minutos (0.1 hora), ya que aunque las medidas de oxígeno, son posibles de obtener a
menores tiempos de muestreo, el comportamiento dinámico en ese intervalo, no tiene
cambios significativos.
Para el modelo del sistema dinámico dentro del filtro, se usó una aproximación
discreta al comportamiento dinámico del bio reactor. En la Figura 5.38 se muestra el
esquema de información utilizado para la aproximación comparativa con los datos
reales.
u t DATOS EXPERIMENTALES y t
FERMENTACIONES BATCH
A-B-C
x̂ t
Figura 5.38. Diagrama de Bloques del Estimador de Biomasa
El modelo usado para los términos de difusión del modelo del filtro, fue una
componente 1 xt , t INC xt , t
con la estructura de la ecuación 4.6 como un modelo
ARIMA en Matlab®, construidos a partir de los errores entre los datos reales y los
datos simulados por el modelo base previamente descrito. Usando como banco de
datos tres fermentaciones “A”, “B” y “C” de diferentes duraciones y con condiciones
iniciales alrededor del mismo punto de operación, con esto se lograría tener una mayor
cercanía entre las trayectorias del proceso fermentativo. En las Figura 5.39 a 5.41
puede observarse claramente los modelos de incertidumbre desarrollados para cada
una de las fermentaciones anteriormente mencionadas.
Los modelos ARIMA, que incluyen términos de promedio móvil, (en inglés Moving
Average, MA), son similares a los modelos de regresión, pero no pueden ser ajustados
por métodos de mínimos cuadrados ordinarios las predicciones son funciones lineales
de los datos pasados, pero son funciones no lineales de los coeficientes, por ejemplo,
un modelo sin términos constantes es un promedio móvil ponderado
exponencialmente, en el cual las predicciones son funciones no lineales de los
parámetros de la MA.
10
-2
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
a. Modelado de Incertidumbre de la concentración de Biomasa para la
fermentación A
13
Un spline es una curva definida en porciones mediante polinomios.
197
Modelo de Sustrato
40
Concentración de Sustrato [g/L]
35
30
25
Modelo sin incertidumbre
20 Dato Real
15
Con Incertidumbre Modelada
10
-5
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
7
Modelo sin Incertidumbre
Dato Real
Oxígeno [g/L]
0
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
Es de notar que los datos reales del oxígeno disuelto fueron obtenidos mediante una
tarjeta de adquisición de datos, su característica ruidosa se debe a la adquisición de la
señal (ver Amicarelli, 2009).
12
10
8
Modelo sin incertidumbre
6 Dato Real
Con Incertidumbre Modelada
0
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
a. Modelado de incertidumbre para la concentración de Biomasa de la
fermentación B.
Modelo de Sustrato
35
Concentración de Sustrato [g/L]
30
25
10
0
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
b. Modelado de incertidumbre de la concentración de sustrato de la fermentación
B
199
-3
x 10
Modelo Oxígeno
10
-2
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
Modelo Sustrato
7
Concentración de Sustrato [g/L]
0
0 1 2 3 4 5 6 7 8
Tiempo [Horas]
x 10
-3 Modelo Oxígeno
8
6
Oxígeno [g/L]
0
0 1 2 3 4 5 6 7 8
Tiempo [Horas]
Figura 5.41. Modelos de incertidumbre para las fermentaciones C
En función de los datos arrojados por las simulaciones de las tres diferentes
fermentaciones (“A”, “B” y “C”) del modelo fenomenológico base en Simulink y los
datos reales se crearon tres tipos de modelos de incertidumbre para cada una de las
variables principales. Con dicha información se pudieron construir los modelos
ARIMA y se combinaron para encontrar el mejor modelo conjunto que pudiera ser
combinado con el modelo fenomenológico establecido.
mismo valor. Tuvo que tenerse en cuenta que el modelo fenomenológico estaba siendo
simulado en tiempo continuo, y que el modelo de incertidumbre debería acoplarse a la
información que éste entregaba. Luego, se tuvo en cuenta este acoplamiento de
modelos con desempeño del filtro y de los datos “muestreados” reales de las
fermentaciones.
z (t ) p( z (t ) z (t 1)), z (t ) U a, b
2 xt , t p( x(t ) z (t t 1), x(t 1)) (3.27)
,
2 xt , t p( x(t ) z (t t 1), x(t 1))dw t
35
Concentración de Sustrato [g/lt]
30
25
20
0
2 4 6 8 10 12 14 16 18
a) Tiempo [Horas]
202
x 10
-3 Entrada para el Filtro - Oxígeno Disuelto
8
0
0 2 4 6 8 10 12 14 16 18 20
b) Tiempo [Horas]
10
0
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
Figura 5.43. Estima de Biomasa para la fermentación B usando como modelo de
incertidumbre el construido a partir de la fermentación B
203
15
Células[g/l]
10
Fermentación Real B
5 Filtro SIR
0
0 2 4 6 8 10 12 14 16 18
Tiempo [Horas]
Figura 5.44. Estima de Biomasa para la fermentación B usando como modelo de
incertidumbre el construido a partir de la fermentación C
2.5
Células[g/l]
2
Fermentación Real C
1.5 Filtro SIR
0.5
0
0 1 2 3 4 5 6 7 8 9
Tiempo [Horas]
Figura 5.45. Estima de Biomasa para la fermentación C usando como modelo de
incertidumbre el construido a partir de la fermentación B
Estimación de Biomasa Modelos Arima C
3.5
2.5
Células[g/l]
0.5
0
0 1 2 3 4 5 6 7 8 9
Tiempo [Horas]
Figura 5.46. Estima de Biomasa para la fermentación C usando como modelo de
incertidumbre el construido a partir de la fermentación C
205
Capítulo 6
CONTROLADORES BASADOS EN MÉTODOS
NUMÉRICOS PARA BIOREACTOR DE
FERMENTACIÓN ALCOHÓLICA EN
CONTINUO
Motivación
Un bioproceso típico hace referencia a reactores por lotes, semi lotes o continuos en
los cuales la materia prima se transforma en los productos biológicos deseados. En
muchas aplicaciones, los preferidos son los bio reactores en continuo, debido a su
mayor posibilidad de manejo en la operación y una alta productividad (Lee, 1992). Un
bioreactor tipo tanque continuamente agitado (continuous stirred tank bioreactor
CSTB) conocido también como fermentador en continuo, se muestra en la Figura 6.1.
Se dice en continuo ya que el medio de cultivo se suministra de forma continua al
reactor para sostener el crecimiento de la población microbiana. El medio sintético
contiene los sustratos metabolizados por las células durante la fase de crecimiento, así
como también otras componentes como minerales y sales que se necesitan para poder
reproducir un ambiente natural para los microorganismos.
del tiempo de residencia del reactor. La tasa de dilución se controla con un lazo simple
de realimentación que manipula la tasa media de flujo. La corriente efluente contiene
sustrato que no reaccionó, trazas de producto y biomasa, que es una mezcla compleja
de células y varios metabolitos producidos por las células. Los productos deseados
pueden ser las células mismas, uno o más metabolitos o alguna combinación de
células y metabolitos. Los productos son separados de los otros componentes vía
series de recuperación, separación y purificación. En adición a la corriente líquida de
producto, también como bioproductos de las reacciones bioquímicas se pueden
producir gases como el dióxido de carbono.
Air
F2,X2,S,D F2,X,S
HLFC
F1,S1,Ds Stream
Water TT
pHT
F3,X3,S3
PC
Separator
F4,X4,S4
F6,X4,S4,Dr
F5,X4,S4
Una operación exitosa de un bioreactor en continuo requiere mucho más que una
estrategia simple para dar los nutrientes necesarios y extraer los productos deseados.
Es esencial que haya una preparación cuidadosa del medio de cultivo ya que los
microorganismos se ven afectados fuertemente por los cambios en el ambiente de
cultivo. El microorganismo debe inocularse en el reactor al iniciar el crecimiento
celular. La inoculación típica se alcanza vía un procedimiento de múltiples pasos en el
que las células crecen en unidades externas y se transfieren a bioreactores de
volúmenes incrementales hasta que el bioreactor de producción es inoculado (Lee,
1992). El procedimiento mencionado es necesario para alcanzar una población lo
suficientemente grande como para sostener el crecimiento. Un requerimiento crítico es
mantener la esterilidad del medio y todo el equipo de procesamiento. Inclusive una
cantidad pequeña de contaminantes puede llevar a una pérdida completa de
productividad y un paro total del bioreactor. Como resultado de esas complejidades, la
operación efectiva de un bioreactor es un problema que representa un reto importante.
207
En este capítulo se propone usar métodos numéricos no solo para simular la evolución
y el comportamiento dinámico de la fermentación alcohólica a partir de Zymomonas
mobilis (Capítulo 5) sino también para encontrar las acciones de control que permitan
llevar las variables de estado del bioreactor desde el estado actual al estado deseado
con el fin de maximizar su productividad. Para el diseño de los controladores se hace
uso del modelo presentado en el capítulo anterior usando la tasa de dilución, la
concentración de sustrato de entrada y la recirculación de microorganismos para el
control de la concentración de sustrato, producto y posteriormente el de biomasa, bajo
condiciones de observabilidad total de las variables primarias del bioproceso. La
contribución de esta parte del trabajo se basa en que no es necesario hacer cálculos
complejos para obtener las acciones de control del bioreactor, tampoco se usan
linealizaciones del modelo dinámico; lo que hace que esta solución sea fácil de
implementar para afrontar el reto de control de este bioreactor en continuo (Quintero,
Scaglia et al., 2008b, 2008e).
La metodología usada fue desarrollada por (Scaglia, 2006) y aplicada sobre robots
móviles para seguimiento de trayectorias (Scaglia et al., 2005)(Scaglia, Quintero, et
al., 2006, 2007, 2008a, 2008b).
x1 f1 ( x1 , x2 , t )u
(6.1)
x2 f 2 ( x1 , x2 , t )u
210
donde x1, x2 representan los estados del sistema, u la señal de control, y t , el tiempo.
Los valores de x1 y x2 en instantes discretos de tiempo t = nT0, donde T0 es el período
de muestreo y n 0,1, 2,3, se simbolizará como x1n y x2n.
( n 1)To
Existen distintos métodos de integración que permiten obtener el valor de x1n+1. Por
ejemplo, el método de Euler (6.3) y el de Runge Kutta de 2do orden (6.4).
T0
x2n 1 x2n f 2 ( x1n , x2n , tn )un f 2 ( x1n 1 , x2n 1 , tn 1 )un 1 (6.6)
2
f1n f1 ( x1n , x2n , tn ) , f1n 1 f1 ( x1n 1 , x2n 1 , tn 1 ) , f2n f2 ( x1n , x2n , tn ) , f2n 1 f2 ( x1n 1 , x2n 1 , tn 1 ) .
2
( x1 x1 )
f1n f1n 1 un T0 n 1 n
(6.7)
f 2n f 2n 1 u n 1 2
( x2 x2n )
T0 n 1
211
Donde la Ecuación (6.7) representa un sistema de dos ecuaciones con dos incógnitas,
las acciones de control un y un+1, si se conocen las trayectorias deseadas por los
estados x1d y x2d, se puede reemplazar x1n+1 y x2n+1 por x1dn+1 y x2dn+1 y de esta manera
calcular un y un+1. Las trayectorias a considerar son continuas por partes y con
derivadas continuas por partes. Del sistema de Ecuaciones (6.7) se puede ver que
debido a la aproximación utilizada, es posible calcular en el instante nT0, el valor de la
señal de control un y un+1, lo cual representa la señal de control en el próximo instante
de muestreo. Esto representa una ventaja importante por dos motivos, primero para
sistemas complejos, las Ecuaciones (6.7) se pueden resolver utilizando métodos
iterativos para la solución de sistemas de ecuaciones lineales, los cuales solo necesitan
un valor inicial para comenzar la iteración que puede ser precisamente la estima
calculada en el instante de muestreo anterior. La segunda es que se pueden combinar
un y un+1 y así eliminar algunas oscilaciones indeseables en la señal de control y
también en los estados.
Asimismo existe una ventaja adicional al plantear el cálculo de las señal de control
como la solución de un sistema de ecuaciones, ya que ahora es más sencillo sacar
conclusiones acerca de los estados que puede alcanzar el sistema, por ejemplo, si la
matriz de coeficientes
f1n f1n 1
A (6.8)
f 2n f 2n 1
Tiene rango 2, entonces la Ecuación (6.7) tendrá solución lo cual indica que existe una
secuencia de control que llevará al sistema al estado deseado. Si el rango de la matriz
T
de coeficientes A es 1, en este caso quiere decir que el vector f1n 1 f 2n 1 se puede
T
obtener como una combinación lineal del vector f1n f 2n y entonces en este caso
se puede utilizar las Ecuaciones (6.3) y (6.5),
x1dn 1 x1n
f1n T0
un (6.9)
f 2n x2 dn 1 x2n
T0
lo cual representa un sistema de dos ecuaciones y una incógnita, cuya solución óptima
según mínimos cuadrados es, (Strang, 1980),
212
Para que el sistema de Ecuaciones (6.9) tenga solución se debe cumplir que,
6.2.1 Controlador 1
Para este propósito, serán usadas las ecuaciones de estado definidas por las
Ecuaciones (5.1) - (5.10). Primero se diseñará un controlador para las variables
principales sustrato y producto, usando la tasa de dilución, la concentración de
sustrato de entrada y eliminando la influencia de la recirculación de microorganismos
propuesta por (Echeverry et al., 2003, 2004).
Sean las ecuaciones que definen el comportamiento dinámico del bioreactor, bajo las
consideraciones anteriores:
dX (6.14)
DS X
dt
dS 1
Qp X DsSin Ds S
(6.15)
dt Yp / s
dP (6.16)
Qp X Ds P
dt
dZ (6.17)
I Z
dt
dI (6.18)
Qp X DsP I
dt
Xn 1 Xn (6.19)
n Dsn X n
T0
Sn 1 Sn 1
Qpn X n Dsn Sinn Dsn Sn
(6.20)
T0 Yp / s
Pn 1 Pn
Q pn X n Dsn Pn (6.21)
T0
Zn 1 Zn (6.22)
In Zn
T0
In 1 In (6.23)
Qpn X n Dsn Pn In
T0
Xn 1 X n T0 n Dsn X n (6.24)
Sn Sn T0
1
Q pn X n Dsn Sinn Dsn S n
(6.25)
1
Yp / s
Zn 1 Z n T0 In Zn (6.27)
In 1 I n T0 Qpn X n Dsn Pn In (6.28)
X n T0 X n n
Xn T0 X n 0 1
1 Sn Q pn X nT0
Sn 1 T0 S n To Yp / s (6.29)
Dsn
Pn 1 T0 Pn 0 Pn T0Q pn X n
Dsn Sinn
Zn 1 0 0 Z n T0 In Zn
In 1 T0 Pn 0
I n T0 Qpn X n In
X n T0 X n n
T0 X n 0 Xn 1
1 Sn Qpn X nT0
T0ToSn To Sn 1
Yp / s (6.30)
Dsn
T0 Pn 0 Pn 1 Pn T0Q pn X n
Dsn Sinn
0 0 Zn 1 Z n T0 In Zn
T0 Pn 0 In 1
I n T0 Q pn X n In
T0 X n 0 1
Sn 1 Sn Q p X nT0
T0 Sn To Yp / s n (6.31)
Ds
n
T0 Pn 0 Pn 1 Pn T0Q p X n
Ds Sin n
0 0 n n
Z n 1 Z n T0 In Zn
T0 Pn 0 U
I n 1 I n T0 Qp Xn In
A n
Xn 1 X n T0 X n n
T0 X n 0 1
Sn 1 Sn Qpn X nT0
T0 S n To Yp / s (6.32)
Dsn
T0 Pn 0 Pn 1 Pn T0 Q pn X n
Dsn Sinn
T0 Pn 0 In I n T0 Q pn X n In
1
0 0
Zn 1 Z n T0 In Zn
Xn 1 X n T0 X n n
T0 X n 0 1
T0 S n To Dsn Sn 1 Sn Qpn X nT0 (6.33)
Yp / s
T0 Pn 0 Dsn Sinn
Pn 1 Pn T0 Q pn X n
T0 Pn 0
In 1 I n T0 Q pn X n In
Xn 1 X n T0 X n n
T0 X n 0 1
Dsn Sn 1 Sn Q pn X nT0 (6.34)
T0 Sn T0 Yp / s
pinv
Dsn Sinn T0 Pn 0
Pn 1 Pn T0Q pn X n
T0 Pn 0
A I n 1 I n T0 Q pn X n In
b
217
Las trayectorias deseadas fueron diseñadas con base en los datos reales del proceso de
fermentación reportado por Raposso (Raposso et al, 2005), que alcanzan los estados
estables a partir de 3 condiciones iniciales reales diferentes.
Los resultados de simulación se presentan en las Figuras 6.2, 6.4 y 6.6 en las cuales se
tiene como trayectorias deseadas las seguidas por los datos reales del trabajo de
Raposso et al, 2005. Las acciones de control usadas están dentro de las condiciones
físicas, definidas por los flujos y la velocidad de actuadores reales. Las acciones de
control están limitadas a valores entre 0.03< Ds<0.2 y 10<Sin<250. Éstas se presentan
en las Figuras 6.3, 6.5 y 6.7 respectivamente.
Parámetros 1 2 3 4 5 6
Cond. Inic. (Daugulis (Raposso (Pinheiro (Raposso (Raposso (Raposso
Osc) Osc) Osc) EE 1) EE 2) EE3)
X(1) (g/l) 3 1.5 1.25 0.05 0.07 0.4
S(1) (g/l) 60 57 60 145 160 150
P(1) (g/l) 70 57 60 0.01 0.01 0.01
Sin(1) (g/l) 187 200 200 187 187 187
Ds(1) (1/h) 0.0879 0.0675 0.067 0.06 0.0675 0.07
Tabla 6.1 Resumen de datos de las simulaciones para prueba de controladores
Puede verse de la Figura 6.2 que el sistema (en línea sólida) sigue de manera aceptable
la trayectoria deseada (en línea punteada); en el tiempo 60 horas el sistema fue
218
Concentración de Producto
Concentración de Etanol [g/L]
60
50
40
30
Modelo
Trayectoria Deseada
20
Datos Experimentales
10
0
0 50 100 150
Tiempo [Horas]
a. Desempeño del controlador de producto, asumiendo observabilidad completa
siguiendo la trayectoria definida por datos experimentales.
219
Concentración de Sustrato
Concentración de Sustrato [g/L]
200
150
Modelo
100 Trayectoria deseada
Datos Experimentales
50
0
0 50 100 150
Tiempo [Horas]
b. Desempeño del controlador de sustrato, asumiendo observabilidad
completa siguiendo la trayectoria definida por datos experimentales.
Figura 6.2 Desempeño del controlador de sustrato y producto, para fermentación 1
150
100
50
0
0 50 100 150
Tiempo [Horas]
a. Acción de Control Sustrato de alimento de la ecuación 6.34 para la
fermentación 1
220
Acción de Control Ds
0.12
Tasa de Dilución [1/h]
0.1
0.08
0.06
0.04
0.02
0 50 100 150
Tiempo [Horas]
b. Acción de control tasa de dilución, de la ecuación 6.34 para la fermentación 1
Figura 6.3 Acciones de control Ds y Sin para fermentación 1
200
Modelo
150
Trayectoria deseada
Datos Experimentales
100
50
0
0 20 40 60 80 100 120
Tiempo [Horas]
a. Concentración de sustrato de alimento dentro del reactor (línea sólida) y
trayectoria deseada (línea punteada) en la fermentación 2
221
Concentración de Producto
Concentración de Etanol [g/L]
70
60
50
40 Modelo
Trayectoria Deseada
30 Datos Experimentales
20
10
0
0 20 40 60 80 100 120
Tiempo [Horas]
b. Concentración de etanol (línea sólida) y trayectoria deseada (línea punteada)
dentro del reactor para la fermentación 2
Figura 6.4 Desempeño del controlador de Sustrato y producto, para fermentación 2
De manera análoga se tiene en la Figura 6.5 las acciones de control que se usaron en la
fermentación 2 para el control de las concentraciones de sustrato y producto de las
Figura 6.4 a y b. Puede notarse que la mayoría del tiempo de la fermentación las
acciones de control estuvieron saturadas, precisamente en la fase batch.
Específicamente se encontró que los parámetros del controlador estaban aumentando
las magnitudes de las acciones de control y estaban perdiendo la suavidad necesaria
para el seguimiento de las trayectorias, lo que sugirió una reducción de las constantes.
Acción de Control Sustrato de Entrada Sin
200
Concentración de Sin [g/L]
150
100
50
0
0 20 40 60 80 100 120
Tiempo [Horas]
a. Acción de control de sustrato de alimento para la fermentación 2.
222
Acción de Control Ds
Tasa de Dilución de sustrato [1/h]
0.12
0.1
0.08
0.06
0.04
0.02
0 20 40 60 80 100 120
Tiempo [Horas]
b. Acción de control de tasa de dilución para la fermentación 2
Figura 6.5 Acciones de control Ds y Sin para fermentación 2
En las Figuras 6.6 y 6.7 se presentan los resultados de la fermentación 3. Puede verse
que la concetración de sustrato de alimento, la concentración de producto siguen lñas
trayectorias deseadas y las acciones de control son mas suaves en el transcurso de la
fermentación, mejorando asi los resltados anterioemente descritos.
Concentración de Sustrato
Concentración de Glucosa [g/L]
200
150 Modelo
Trayectoria Deseada
Datos Experimentales
100
50
0
0 20 40 60 80 100 120
Tiempo [Horas]
a. Concentración de sustrato de alimento para la fermentación 3
223
Concentración de Producto
Concentración de etanol [g/L]
70
60
50
40
Modelo
30 Trayectoria Deseada
Datos Experimentales
20
10
0
0 20 40 60 80 100 120
Tiempo [Horas]
b. Concentración de producto dentro del reactor para la fermentación 3
Figura 6.6 Desempeño del controlador de Sustrato y producto, para fermentación 3
Acción de Control Sustrato de Entrada Sin
200
Concentración de Sin [g/L]
150
100
50
0
0 20 40 60 80 100 120
Tiempo [Horas]
a. Acción de control de sustrato de alimento para la fermentación 3
Acción de Control Ds
Tasa de Dilución de sustrato [g/L]
0.12
0.1
0.08
0.06
0.04
0.02
0 20 40 60 80 100 120
Tiempo [Horas]
b. Acción de control de la tasa de dilución para la fermentación 3
Figura 6.7 Acciones de control Ds y Sin para fermentación 3
224
6.2.2 Controlador 2
dX (6.35)
X RD / 4 DS R 1 X
dt
Xn 1 Xn T X n T0 X n RD / 4 Dsn R 1
n 0
(6.36)
Donde el término D es la tasa de dilución total, esto significa que es la suma de la tasa
de dilución de sustrato Ds y la tasa de dilución Dr asociada a la recirculación de
microorganismos R. Es importante recalcar que los términos relaiconados con la
recirculación R permanecerán constantes porque el propósito es analizar su influencia
sobre el valor de las acciones de control Ds y Sin.
Xn 1 Xn n 0 T Xn T0 X n RD / 4 T0 X n Dsn R 1
Xn 1 Xn n 0 T Xn 2
T0 X n Dsn R T0 X n R D / 4 T0 X n Dsn T0 X n RD / 4 (6.37)
2 2
Xn 1 Xn n 0 T X n T0 X n R Dr / 4 T0 X n R Dsn / 4 T0 X n Dsn R T0 X n Dsn T0 X n RDsn / 4 T0 X n RDr / 4
2
Xn 1 Xn n 0 T Xn Dsn T0 X n R / 4 T0 X n R T0 X n T0 X n R / 4 Dr T0 X n R 2 / 4 T0 X n R / 4
3
Xn 1 Xn n ToX n Dsn T0 X n R 2 / 4 T0 X n R T0 X n Dr T0 X n R 2 / 4 T0 X n R / 4
4
Entonces
3
Dsn T0 X n R 2 / 4 T0 X n R T0 X n Xn 1 Xn T Xn
n 0 Dr T0 X n R 2 / 4 T0 X n R / 4
4
(6.38)
Xn 1 X n 1 T0 n Dr T0 X n R 2 / 4 T0 X n R / 4
2 3
T0 X n R / 4 ToX n R T0 X n 0
4 1
Dsn Sn 1 Sn Q pn X nT0
T0 S n To Yp / s
Dsn Sinn
T0 Pn 0 Pn 1 Pn T0Q pn X n
T0 Pn 0
In 1 I n T0 Q pn X n In
(6.39)
Xn 1 X n 1 T0 n Dr T0 X n R 2 / 4 T0 X n R / 4
3
T0 X n R 2 / 4 ToX n R T0 X n 0
4 1
Dsn Sn 1 Sn Q pn X nT0
pinv T0 S n To Yp / s
Dsn Sinn
T0 Pn 0 Pn 1 Pn T0 Q pn X n
T0 Pn 0
In 1 I n T0 Q pn X n In
(6.40)
El desempeño del controlador 2 con una recirculación fija del 10%, puede verse de la
Figura 6.8. Como caso de estudio y al igual que en el controlador 1, las trayectorias
deseadas fueron generadas mediante interpolación y suavizado de las trayectorias
definidas con base en los datos experimentales de (Raposso et al., 2005) para Sustrato
y Producto. Solo se presentará los resultados de una fermentación contrario a las tres
presentadas en la sección anterior, como ilustración, dado que este controlador
226
60
Concentración de Producto
Concentración de Etanol [g/L]
50
40
30 Modelo
Trayectoria Deseada
20 Datos Experimentales
10
0
0 20 40 60 80 100 120 140 160
Tiempo [Horas]
a. Concentración de Producto en la fermentación controlada (línea solida)
trayectoria deseada (línea punteada)
Concentración de Sustrato
Concentración de Glucosa [g/L]
200
150
Modelo
100 Trayectoria Deseada
Datos Experimentales
50
0
0 20 40 60 80 100 120 140 160
Tiempo [Horas]
b. Concentración de sustrato en la fermentación controlada (línea solida)
trayectoria deseada (línea punteada)
227
200
Acción de Control Sustrato de Entrada Sin
Concentración de Sin [g/L]
150
100
50
0
0 20 40 60 80 100 120 140 160
Tiempo [Horas]
c. Acción de control concentración de sustrato de alimento
200
Acción de Control Sustrato de Entrada Sin
Concentración de Sin [g/L]
150
100
50
0
0 20 40 60 80 100 120 140 160
Tiempo [Horas]
d. Acción de control de tasa de dilución.
Figura. 6.8. Desempeño del controlador 2 ante perturbaciones.
Debe recalcarse que la diferencia del controlador 2 de esta sección y la del controlador
1 de la sección 6.2.1 radica esencialmente en la inclusión del término asociado a la
recirculación de microorganismos al reactor, lo que hace que halla modificaciones en
las dinámicas de las concentraciones de microorganismos que afectan
considerablemente las concentraciones de sustrato y de producto dentro del reactor.
En las Figura 6.8 c. y d. se puede ver claramente que las acciones de control
responden dinámicamente a los cambios presentados dentro del reactor y que dado
que dichas variables manipuladas tienen un efecto directo y más rápido sobre la
concentración de sustrato, es ésta la que se mantiene mejor controlada durante toda la
fermentación.
De los resultados de las Figura 6.9, puede concluirse que la metodología propuesta es
simple en el sentido de la selección de los parámetros del controlador para alcanzar un
228
buen desempeño del sistema en lazo cerrado. Esta metodología puede ser aplicada a
otros tipos de sistemas. La precisión requerida del método numérico propuesto para la
aproximación del modelo es menor que la necesaria para simular el comportamiento
del sistema. Esto se debe a que se consideró que la información de los estados
realimentados está disponible en cada instante de muestreo y es corregida cualquier
diferencia debido a errores acumulativos (por ejemplo, errores de redondeo). Como
consecuencia, la aproximación usada para encontrar la mejor forma de ir de un estado
a otro, puede variar de acuerdo a la disponibilidad de la información y del modelo del
sistema. El diseño del controlador puede plantearse como la minimización de un
índice cuadrático, lo que es un problema simple y permite considerar otras
propiedades de las trayectorias.
Concentración de Sustrato
Concentración de Glucosa [g/L]
200
150 Modelo
Trayectoria deseada
Datos Experimentales
100
50
0
0 20 40 60 80 100 120
Tiempo [Horas]
a. Concentración de sustrato en la fermentación controlada (línea solida)
trayectoria deseada (línea punteada)
Concentración de Producto
Concentración de etanol [g/L]
70
60
50
Modelo
40 Trayectoria Deseada
Datos Experimentales
30
20
10
0
0 20 40 60 80 100 120
Tiempo [Horas]
b. Concentración de Producto en la fermentación controlada (línea solida)
trayectoria deseada (línea punteada)
229
0.12
Acción de Control Ds
Tasa de Dilución [1/h]
0.1
0.08
0.06
0.04
0.02
0 20 40 60 80 100 120
Tiempo [Horas]
c. Acción de control de la tasa de dilución
200
Acción de Control Sustrato de Entrada Sin
Concentración de Sin [g/L]
150
100
50
0
0 20 40 60 80 100 120
Tiempo [Horas]
6.2.3 Controlador 3
dX RD
X DS R 1 X
dt 4 (6.41)
dS 1
Qp X DSin DS
dt Yp / s
(6.42)
dP
Qp X DP
dt (6.43)
dZ
I Z
dt (6.44)
dI
Qp X DP I
dt (6.45)
De la misma manera que los controladores anteriores, los valores de las salidas en el
tiempo discreto t nT0 , donde T0 es el periodo de muestreo y n 0,1, 2,3, , serán
los valores de X n , S n y Pn . Si se desea conocer el valor de X n 1 , Sn 1 y
Pn 1 conociendo los valores de X n , S n y Pn , el conjunto de ecuaciones previamente
descritas deben integrarse en el intervalo de tiempo nT0 t (n 1)T0 . En
consecuencia, la aproximación de Euler para calcular los valores de
X n 1 , Sn 1, Pn 1 , Z n 1 , I n 1 es la siguiente:
Xn 1 Xn Rn Dn
n Xn Dsn Rn 1 X n
T0 4
(6.46)
231
Sn 1 Sn 1
Qpn X n Dn Sinn Dn Sn
T0 Yp / s
(6.47)
Pn 1 Pn
Qpn X n Dn Pn
T0
(6.48)
Zn 1 Zn
In Zn
T0
(6.49)
In 1 In
Qpn X n Dn Pn In
T0
(6.50)
T0 X n Rn 2 Dn T0 X n Rn Dn
Xn 1 X n 1 T0 n T0 X n Rn Dsn T0 X n Dsn
4 4 (6.51)
1
Sn 1 Sn T0 Q pn X n Dn Sinn Dn S n
Yp / s
(6.52)
Pn 1 Pn T0 Qpn X n Dn Pn
(6.53)
232
Zn 1 Z n T0 In Zn
(6.54)
In 1 I n T0 Qpn X n Dn Pn In
(6.55)
T0 X n Rn 2 Dn T0 X n Rn Dn
X n 1 T0 n T0 X n Rn Dsn T0 X n Dsn
4 4
Xn 1
1
Sn T0 Q pn X n T0 Dn Sinn T0 Dn S n
Sn 1
Yp / s
Pn 1 Pn T0Q pn X n T0 Dn Pn
Zn 1 Z n T0 In Zn
In 1
In 1 T0 T0 Q pn X n T0 Dn Pn
(6.56)
Rn 2 Dn X n 1 T0 n
T0 X n T0 X n
Xn 1 T0 X n T0 X n 0 0
4 4 Rn Dsn 1
Sn S n T0 Q pn X n
1 Yp / s
0 0 0 0 T0 T0 S n Rn Dn
Pn 1
0 0 0 0 0 T0 Pn Dsn Pn T0 Q pn X n
Zn 1
0 0 0 0 0 0 Dn Sinn Z n T0 In Zn
In 1 0 0 0 0 0 T0 Pn
Dn In 1 T0 T0 Q pn X n
(6.57)
(6.58)
Re ordenando las filas 4 y 5
233
Xn 1 X n 1 T0 n
T0 X n T0 X n Rn 2 Dn 1
T0 X n T0 X n 0 0 Sn 1 S n T0 Q pn X n
4 4 Rn Dsn Yp / s
0 0 0 0 T0 T0 S n Rn Dn
Pn 1 Pn T0 Q pn X n
0 0 0 0 0 T0 Pn Dsn
0 0 0 0 0 T0 Pn Dn Sinn In 1 In 1 T0 T0 Q pn X n
0 0 0 0 0 0 Dn Zn Z n T0 In Zn
1
(6.59)
Reduciendo la matriz A mediante la eliminación de la variable de estado Z, se obtiene
la siguiente expresión
Rn 2 Dn Xn 1 X n 1 T0 n
T0 X n T0 X n
T0 X n T0 X n 0 0 Rn Dsn 1
4 4 Sn S n T0 Q pn X n
Rn Dn 1
Yp / s
0 0 0 0 T0 T0 S n
Dsn
0 0 0 0 0 T0 Pn Pn Pn T0 Q pn X n
1
0 0 0 0 0 T0 Pn Dn Sinn
Dn In 1 In 1 T0 T0 Q pn X n
(6.60)
Estas condiciones son similares a las de la sección 6.2.1. la solución óptima según
mínimos cuadrados esta dada por
Rn 2 Dn Xn 1 X n 1 T0 n
T0 X n T0 X n
Rn Dsn T0 X n T0 X n 0 0 1
4 4 Sn 1 Sn T0 Q pn X n
Rn Dn Yp / s
0 0 0 0 T0 T0 Sn
Dsn
0 0 0 0 0 T0 Pn Pn 1 Pn T0Q pn X n
Dn Sinn 0 0 0 0 0 T0 Pn
Dn In 1 In 1 T0 T0 Q pn X n
(6.61)
posible hacer un análisis dinámico que permita ajustar adecuadamente las acciones de
control y sus efectos sobre el desempeño del bioreactor y la evolución temporal del
mismo. De manera similar a los casos anteriores será necesario probar la capacidad
del controlador de responder ante perturbaciones considerables en las corrientes de
entrada del fermentador.
El escenario de prueba será una simulación de 150 horas de duración aproximada, que
busca reproducir las condiciones de los datos reales de la fermentación reportada por
(Raposso et al., 2005) usando sus condiciones iniciales y los parámetros calculados
por ellos.
Se desea que el controlador siga una trayectoria natural, pero que no contenga picos
indeseados ni oscilaciones que puedan afectar la productividad del microorganismo, ni
lo sometan a condiciones de estrés. Por ellos las trayectorias deseadas son más suaves
que las seguidas en el reporte de Raposso. Estas condiciones de simulación nos
permiten analizar la factibilidad de uso de un controlador para su implementación en
línea, consecuente con la concepción de usar la estrategia de alimentación como
herramienta para la maximización de la productividad de etanol.
Los resultados que se muestran a continuación en las Figura 6.10 son satisfactorios ya
que el controlador tiene un buen desempeño y las variables biomasa, sustrato y
producto siguen las trayectorias dentro de los límites permitidos como aceptables (5%
alrededor del valor deseado). Las acciones de control usadas que se ven en las Figura
6.11, se limitan a los valores condicionados por los actuadores y la capacidad física de
la bacteria, esto es 0.03< Ds<0.2, 10<Sin<250 y para el término de recirculación las
condiciones de 0<R<1.
235
2
Biomasa [g/L]
Modelo
1.5 Trayectoria Deseada
Datos Experimentales
1
0.5
0 20 40 60 80 100 120
Tiempo [Horas]
200
Glucosa[g/L]
Modelo
Trayectoria Deseada
100 Datos Experimentales
0
0 50 Tiempo [Horas] 100 150
80
Etanol [g/L]
60
40
Modelo
20 Trayectoria Deseada
Datos Experimentales
0
0 50 Tiempo [Horas] 100 150
0.5
0
0 50 Tiempo [Horas] 100 150
Concentración de Sin [g/L]
100
0
0 50 Tiempo [Horas] 100 150
0.1
0.05
0
0 50 Tiempo [Horas] 100 150
Acción de Control Tasa de Dilución Total D
Tasa de dilución total [1/h]
0.2
0.1
0
0 50 100 150
Tiempo [Horas]
Figura 6.11. Acciones de control para seguimiento de trayectoria del controlador
2
Biomasa [g/L]
Modelo
1.5 Trayectoria Deseada
1 Datos Experimentales
0.5
0 50 Tiempo [Horas] 100 150
200
Glucosa [g/L]
Modelo
Trayectoria Deseada
100 Datos Experimentales
0
0 50 Tiempo [Horas] 100 150
Etanol [g/L]
50
Modelo
Trayectoria Deseada
Datos Experimentales
0
0 50 Tiempo [Horas] 100 150
a. Variables de control Biomasa, Susutrato y Producto
Acción de Control de Recirculación de Biomasa
1
Fracción
0.5
0
0 50 Tiempo [Horas] 100 150
Acción de Control Sustrato de Entrada Sin
Concentración de Sin [g/L]
200
100
0
0 50 Tiempo [Horas] 100 150
Acción de Control Tasa de Dilución de Sustrato Ds
0.1
Tasa de Dilución [1/h]
0.05
0 Tiempo [Horas]
0 50 100 150
Acción de Control Tasa de Dilución Total D
0.2
Tasa de Dilución [1/h]
0.1
0
0 50 Tiempo [Horas] 100 150
El uso de las trayectorias basadas en datos reales mostró que el controlador es factible
y puede ser implementado con acciones de control acotadas para las especificaciones
del proceso real, el término de recirculación que fue agregado en el controlador 3
como acción de control representa la probabilidad de usar la recirculación de biomasa
como variable de control para la fermentación en continuo, y que puede mejorar el
comportamiento dinámico. De los resultados obtenidos en simulación puede verse que
el error entre la trayectoria deseada y la seguida por el proceso es pequeño, dentro de
los límites tolerables alrededor del valor deseado.
239
MODELO DE
FERMENTADOR
ALCOHÓLICO
ut yt
ESTIMADOR DE
ESTADOS FILTRO DE
PARTÍCULAS
xt
CONTROLADOR
BASADO EN MÉTODOS
NUMÉRICOS
Este caso nos ubica en el uso del estimador de biomasa en lazo cerrado con los
controladores de las secciones 6.2.1 y 6.2.2 ya que la tasa de recirculación de
microorganismos puede permanecer fija, inclusive en un valor nulo. En la Figura 6.14
puede verse la simulación del controlador de sustrato y producto, en lazo cerrado con
el estimador de estados, se grafican únicamente la trayectoria real y la trayectoria
deseada. En esta simulación se usa una recirculación de microorganismos del 10%
160 80
Prod Model
ProdDesired
Prod Exp Data
70
132
60
Glucose Concentration [g/L]
40
76
Sust Model 30
SustDesired
Sust Exp Data
20
48
10
20 0
0 30 50 60 90 100 120 150
Time [Hours]
La Figura 6.15 presenta el mismo controlador de sustrato y de producto que usa como
deseada, la trayectoria calculada al igual que en la sección 6.2 basada en datos reales.
A través del análisis de estas simulaciones, puede concluirse que el desempeño en lazo
cerrado de ambas soluciones propuestas es aceptable y cumple con las expectativas de
control y de estimación. Es importante resaltar que ambas aproximaciones se realizan
en el tiempo de muestreo simulado de 6 minutos, tiempo suficiente para realizar no
solo las acciones de control sino el cálculo del estimador recursivo inclusive con el
número máximo de partículas usado en este trabajo.
241
Se observó que el uso del estimador basado en filtros de partículas para las variables
de biomasa e inhibición es satisfactorio y provee información que el controlador usa
para llevar el sistema al estado deseado. Este primer resultado supone que ambas
herramientas son implementables en el proceso de fermentación real. El uso de la
herramienta de estimación permite resolver el problema de falta de información de la
concentración de biomasa en línea y las demás variables estimadas permiten
desarrollar estrategias de mayor jerarquía en caso de una aplicación a escala industrial.
Computacionalmente esta metodología puede ser implementada en una computadora
personal de características comunes con acceso a puertos de comunicación con
dispositivos electrónicos de medida.
Concentración de Glucosa [g/L]
60
148
50
116 40
Modelo de Producto
30
Prod Deseado
84 Modelo de Sustrato
Prod Dato Experimental 20
Sust Deseado
Sust Dato Experimental 10
52
0
20 -10
0 30 50 60 90 100 120 150
Tiempo [Horas]
Retomando el control de la sección 6.2.3, debe recordarse que las acciones de control
R, Ds y Sin están acopladas. En este caso se usa la recirculación de microorganismos
como variable manipulada para hacer que la concentración de biomasa siga una
trayectoria previamente definida. En esta sección en particular se verá con mayor
relevancia, la utilidad del estimador de biomasa, ya que no solo servirá para calcular
acciones de control para el sustrato y el producto, sino también que esta medida está
involucrada directamente en la acciones para controlar la cantidad de
microorganismos de la colonia, lo que exige un cálculo más preciso y un control que
pueda ser menos susceptible a la pérdida de precisión de la medida o los errores de
redondeo, característica con que cuenta el controlador basado en métodos numéricos
(Scaglia, 2006)(Scaglia, Quintero et al., 2006, 2007, 2008a)
En la Figura 6.16 se presentan las trayectorias seguidas por el bioproceso para una
simulación de 150 horas de duración, reproduciendo los datos reales de (Raposso et
al., 2005) bajo las consideraciones de la disponibilidad de información únicamente de
las concentraciones de sustrato y producto, extraídos de los datos reales.
1.8
1.6
1.4
1.2 Fermentación
1 Trajectoria Deseada
0.8
0.6
0.4
0.2
Concentración de Sustrato
Concentración de Glucosa [g/L]
300
250
200
Fermentación
150
Trajectoria Deseada
100
50
0
0 50 Tiempo [Horas] 100 150
80
Fermentación
70 Trayectoria deseada
60
50
40
30
20
10
0
0 50 Tiempo [Horas] 100 150
Capítulo 7
CONCLUSIONES
2009
Revistas Indexadas:
En primera Revisión
2008
Revistas Indexadas
Congresos:
2007
Congresos:
G. Scaglia,G., V. Mut, O Quintero et al,. “Tracking Control of a Mobile
Robot using LinearInterpolation” IMAACA, 2007 vol. 1, pp. 11-15, ISBN:
978-2-952071277.
O. Quintero, A. Amicarelli, F. di Sciascio, Estimador de Estados en
Fermentación Alcohólica en Continuo de Z.m Mediante Filtrado Bayesiano
Recursivo .XII Reunión de Trabajo en Procesamiento de la Información y
Control, 2007, Río Gallegos, Argentina
2006
Congresos:
G. Scaglia, O. Quintero, V Mut, F. di Sciascio, “Control de Trayectoria de
Robots Moviles Usando Álgebra lineal”. XII Congreso Latinoamericano de
Control Automático CLCA 2006, Bahia, Brazil.
G. Scaglia, O. Quintero, V Mut, F. di Sciascio, “Control de Trayectoria de
Robots Moviles Usando el método de integración trapezoidal”. XIX
Congreso Argentino de Control Automático AADECA 2004, Buenos Aires,
Argentina, 2006.
2005
O. Quintero, F. di Sciascio, “Observadores De Estados De Un Proceso De
Fermentación Alcohólica En Continuo Mediante Filtros De Kalman
Extendidos” XI Reunión de Trabajo en Procesamiento de la Información y
Control, 21 al 23 de septiembre de 2005, Río IV, Argentina.
249
Bibliografía
Adilson J., Rubens M., “Soft sensors development for on-line bio-reactor state
estimation.” Computers and Chemical Engineering 24 , 1099-1103, (2000)
Ahlstrom, M. and Calais, M. “Bayesian terrain navigation with Monte Carlo Method”.
Master's Thesis LiTH-ISY-EX-3051, Department of Electrical Engineering.
Linkoping University, Linkoping, Sweden, 2000. In Swedish.
Amicarelli A., di Sciascio F. and Álvarez H., “Including Dissolved Oxygen Dynamics
to the bacillus thuringiensis δ-endotoxins production process”. Manuscript submitted
to The Canadian Journal of Chemical Engineering, (2007).
Andersen, D. C.; Swartz, J.; Ryll, T.; Lin, N. e Snedecor, B. (2001), Metabolic
oscillations in an E. coli fermentation Dana C. Andersen, Biotechnology and
Bioengineering,, v. 72, n. 2, p. 212 – 218.
Andrieu, C. and Godsill, S. J. (2000). A particle filter for model based audio source
separation. In Proceedings of the International Workshop on Independent Component
Analysis and Blind Signal Separation (ICA 2000), Helsinki, Finland.
Andrieu, C. and Doucet A. (2000) Stochastic algorithms for marginal MAP retrieval
of sinusoids in non-Gaussian noise, Proc. Work. IEEE SSAP.
250
Andrieu, C. and Doucet, A. “Particle filtering for partially observed Gaussian state
space models”. Journal of the Royal Statistical Society, 64(4):827–836, 2002.
Arnold, L. and Xu Kedai (1995). Normal forms for random differential equations.
Journal of Differential Equations 116, 484–503.
Arnold, L. and P. Imkeller (1998). Normal forms for stochastic differential equations.
Probabability Theory and Related Fields 110, 559–588.
Arulampalam M. S., Maskell S., Gordon N., and Clapp T., “A Tutorial on Particle
Filters for Online Nonlinear/Non-Gaussian Bayesian tracking,” IEEE Trans. Signal
Processing, vol. 50, no. 2, pp. 174–188, Feb. 2002.
Åström, K.J. Introduction to Stochastic Control Theory. Academic Press. New York.
(1970)
Atehortúa P., Álvarez H., Orduz S. “Comments on: A Sporulation Kinetic Model for
Batch Growth of B. thuringiensis”, The Canadian Journal of Chemical Engineering,
V. 84, No. 3. (2006).
Atehortua P., Alvarez H.D., S. Orduz. Modeling of growth and sporulation of Bacillus
thuringiensis in an intermittent fed batch culture with total cell retention. Bioprocess
Biosyst Eng. 2007 Jul 22. (2007)
Ba, S.O., and Odobez, J.M., “A Rao-Blackwellized Mixed state particle filter for head
pose tracking in meetings”, Workshop MMMMP, ICMI Trento, October 2005.
Bar-Shalom, Y. and Li, X.-R. “Estimation and Tracking: Principles, Techniques, and
Software”. Artech House, Norwood, MA, USA. 1993.
Bartholomaus, R., “Continuation methods for robust controller design for plants with
operating-point-dependent behaviour” Contr. Eng. Practice, 10, 941. 2002.
Benes, V. “New exact finite dimensional with large Lie algebras”Syst. Control Lett,
vol 5, pp 217-221, 1985.
Berger, J.O. Statistical Decision Theory and Bayesian Analysis. Springer- Verlag,
2nd edition, 1985
Beuse, M.; Kopmann, A.; Diekmann, H.; Thoma, M. (1999), Oxygen, pH value, and
carbon source induced changes of the mode of oscillation in synchronous continuous
culture of Saccharomyces cerevisiae, Biotechnology and Bioengineering, v. 63, n. 4,
p. 410 – 417.
Bishop, G. and Welch, G. “An Introduction to the Kalman Filter”. University of North
Carolina at Chapel Hill Department of Computer Science, SIGGRAPH 2001, Course
8.
Boillereaux L., Flaus J., “A New Approach for Designing Model-Based Indirect
Sensors” IEEE Transaction on Control Systems Technology, Vol. 8, N. 4, July
(2000).
Box, G. and Tiao, G. “Bayesian Inference in Statistical Analysis”. Ed. John Wiley,
1992.
Bravo, S, Mahn, A and Shene, C. Effect of feeding strategy on Zymomonas mobilis
CP4 fed-batch fermentations and matemática modelling of the system.. Appl
Microbiol Biotechnol (2000) 54: 487-493
Briers, M., Doucet, A, and Maskell, S. (2004). “Smoothing algorithms for state-space
models,” Cambridge University Engineering Department Technical Report, CUED/F-
INFENG/TR.498.
Briers, M., Doucet, A., and Singh, S. (2005). “Sequential auxiliary particle belief
propagation,” International Conference on Information Fusion.
Bringer, S.; Hartner, T.; Poralla, K. e Hermann, S. (1985), Influence of ethanol on the
hopanoid content and the fatty acid pattern in batch and continuous cultures of
Zymomonas mobilis. Archives of Microbiology, v. 140, p. 312–316.
Briers M., Doucet, A., “Sequential Auxiliary Particle belief Propagation” Cambridge
University Engineersing department 2005.
Bryson E. and KI. Frazier. “Smoothing for linear and non- I . linear dynamic
systems,” Aeronautical SGstems Div Patterson AFB, Ohio, Tech. Rept. ASD-TDR-
63-119, February 196.3.
Bryson, A.E. Jr., and Henrikson, L.J. “Estimation using sampled data containing
sequentially correlated nose” in Journal of Spacecraft and Rockets, June 1968
Bucy R.S. “Non Linear Filtering Theory”. IEEE Transactions on Automatic Control
1965.
Bucy, R.S. and P.D. Joseph (1968). Filtering for Stochastic Process with Applications
to Guidance. Wiley-Interscience.
Cai, Z., Le Gland F., and Zhang H., “An adaptive local grid refinement method for
nonlinear filtering”. Internal Pubication. INRIA, Tech. Rep. 954.
Carpenter J., Clifford P., and Fearnhead P., “Improved Particle Filter for Nonlinear
Problems,” IEE Proc. F Radar, Sonar Navigation., vol. 146, no. 1, pp. 2–7, 1999.
Casella G., “Statistical inference and Monte Carlo algorithms” Test, vol. 5, pp. 249–
344, 1997.
253
Chang C., and Ansari, R., “Kernel particle filter: Iterative Sampling for efficient visual
tracking” IEEE 2003.
Chen T, Morris J and Martin E. “Particle Filters for the Estimation of a State Space
Model” Centre for Process Analytics and Control Technology School of Chemical
Engineering and Advanced Materials University of Newcastle, Newcastle upon Tyne,
NE1 7RU, UK 2003.
Chen, R., Wang, X., and Liu, J. S. (2000). Adaptive joint detection in flat-fading
channels via mixture Kalman filtering. IEEE Transactions on Information Theory,
46(6):2079– 2094.
Chen, R. and Liu, J. S. “Mixture Kalman filters”. Journal of the Royal Statistical
Society, 62(3):493–508. 2000.
Chen, Z “Bayesian Filtering: From Kalman Filters to Particle Filters, and Beyond”,
Technical report, the Communications Research Laboratory, McMaster University,
Hamilton, Ontario, Canada 2003
Chi, C. T.; Howell, J. A. e Pawlowsk, U., (1974), Regions of multiple stable steady
states of a biological reactor with wall growth, Chemical Engineering Science, v. 29,
p. 207-211.
Chen, Z and Haykin, S. “On different facets of regularization theory”, Neural comput
vol 14, No 12 pp 279-2846, 2002.
Clancy EA, Hogan N. Single site electromyograph amplitude estimation. IEEE Trans
Biomed Eng 41: 159–167, 1994.
Clancy EA, Hogan N. Probability density of the surface electromyogram and its
relation to amplitude
detectors. IEEE Trans Biomed Eng 46: 730–739, 1999.[CrossRef][ISI][Medline]
Cox, H. “On the estimation of state variables and parameters for noisy dynamic
systems,” IEEE Trans. on Automatic Control, vol. AC-9, pp. 5-12, January 196:.
Daum, F. “Exact finite dimensional nonlinear filter” IEEE Trans, Automatic, Control,
vol 31 no 7, pp 616 – 622, 1986.
Daum, F. “New exact nonlinear filters”Bayesian analysis of time series and dynamic
models, J.C Spall, Ed New York: Mariel Dekker, 1988. Pp 199- 226.
Daum F, Nonlinear Filters: beyond The Kalman Filter, IEEE A& E Systems
Magazine, 2005
Dien, B. S.; Cotta, M. A. e Jeffries, T. W. (2003), Bacteria engineered for fuel ethanol
production: current status, Applied Microbiology and Biotechnology, v. 63, n. 3, p.
258 – 266
256
Doucet A. “Monte Carlo Methods for Bayesian Estimation of Hidden Markov Models.
Application to Radiation Signals” Ph.D. Thesis (in French with chapters 4 and 5 in
English) Univ. Paris-Sud, Orsay, 1997.
Doucet A., Godsill S., and Andrieu C. “On Sequential Monte Carlo Sampling
Methods for Bayesian Filtering”, 1998.
Doucet, A., Godsill, S. J., and Andrieu, C. (2000a). On sequential Monte Carlo
sampling methods for Bayesian filtering. Statistics and Computing, 10(3):197–208.
Doucet, A., Godsill, S. J., and West, M. (2000b). Monte Carlo filtering and smoothing
with application to time-varying spectral estimation. In Proceedings of the IEEE
International Conference on Acoustics, Speech, and Signal Processing, volume 2,
pages 701–704, Istanbul, Turkey.
Doucet A., Gordon N., and Krishnamurthy V., “Particle filters for state estimation of
jump Markov linear systems,” IEEE Trans. Signal Processing, vol. 49, pp. 613–624,
Mar. 2001.
257
Doucet, A., Briers, M. “Efficient Block Sampling strategies for sequential Monte
Carlo Methods” Cambridge University Engineersing Department 2006.
Dondo, R. y Marques, D., “Mass and energy balances as state-space models for
aerobic batch fermentations.,” Lat. Am. Appl. Res.. V. 32, no.2 p.195-204,(2002).
Eidehall, A., Schön, T. B., and Gustafsson, F. (2005). The marginalized particle filter
for automotive tracking applications. In Proceedings of the IEEE Intelligent Vehicle
Symposium, pages 369–374, Las Vegas, USA.
Fang B., “Kalman-Bucy Filter for Optimum Radio- Inertial Navigation” Short Paper
IEEE Transactions on Automatic control, August 1967.
258
Fearnhead, Paul. Sequential Monte Carlo methods in Filter theory, Merton College
University of Oxford. A thesis submitted in partial fulfillment of the requirements for
the degree of DPhil at the University of Oxford Trinity Term 1998.
Fearnhead, Paul. “Particle filters for mixture models with unknown number of
components,” paper preprint, 2001. Available on line
http://www.maths.lancs.ac.uk/_fearnhea/.
Fujisaki M., Kallianpur G., and Kunita H. “Stochastic differential equations for the
nonlinear filtering problem,” Osaka J. Math., vol. 9, pp. 19–40, 1972.
Gadeyne K., and Lefebvre T., “The application of probabilistic techniques for the
state/parameter estimation of (dynamical) systems and pattern recognition problems”
Katholieke Universitet Leuven, 2004.
Gentil, I., Remillard, B. and DelMoral, P. (2005). Filtering of images for detecting
multiple targets trajectories. In Statistical Modeling and Analysis for Complex Data
Problems. Springer. To Appear.
Geweke, J. y Tanizaki, H. “On Markov chain Monte Carlo methods for nonlinear and
non-gaussian state-space models,” Commun. Stat. Simul. C, vol. 28, pp. 867–894,
1999.
259
Gillespie Dan T.. The chemical langevin equation. J. Chem. Phys.,113:297–306, 2000.
Glynn P., and Szechtman, R., Some perspectives on the method opf control variates.
In Monte Carlo, cuasi Monte Carlo Methods. Springler Verlag 2002
Gordon, N., Salmond, D., and Smith, A. “Novel Approach to Non Linear-Non
Gaussian Bayesian State Estimation” IEE Proceedings F, Abril 1993.
Guergachi, A, Patry, G “Using statistical learning theory for modeling the uncertainty
in Bussiness and engineering systems: A qualitative Introduction”, Systems, Man, and
Cybernetics, 2001 IEEE International Conference onVolume 1, 7-10 Oct. 2001
Page(s):423 - 428 vol.1 Digital Object Identifier 10.1109/ICSMC.2001.969849
Haddad W., Challaboina V., Nersesov S., “Time Reserval Symmetry, Poincare
Recurrence, Irreversibility and Entropic Arrow of time: From Mechanics to System
Thermodynamics”. Proc. 44th IEEE CDC-ECC, España, pp. 5995-6002 (2005).
Haykin, S., Yee, P., and Derbez, E., “Optimum Nonlinear Filtering” IEEE
Transactions on Signal processing, Vol 45, No 11, November 1977.
Heffes, H. “The effect of erroneous models on the Kalman Filter response”. IEEE
Transactions on Automatic Control, July 1966.
Higdon D. M., “Auxiliary variable methods for Markov chain Monte Carlo with
applications,” J. Amer. Statist. Assoc., vol. 93, pp. 585–595, 1998.
Ho, Y.C y Lee, R.C. “A Bayesian Approach to problems in stochastic Estimation and
Control” IEEE Transactions on Automatic Control, October 1964.
Jazwinski, A.H. “Stochastic Processes and Filtering Theory”. Academic Press, 1970.
Johansen A., Singh, S., Doucet, A., Vo, B-N., “Convergence of the SMC
Implementation on the PHD Filter” Signal Procesing Laboratory Department of
Engineering University of Cambridge, 2005.
Julier, S. J., Uhlmann, J. K., and Durrant-Whyte, H. F. “A new approach for filtering
nonlinear systems”. In Proceedings of the 1995 American Control Conference (pp.
1628–1632), 1995.
Julier, S. J., and Uhlmann, J. K. “A new extension of the Kalman filter to nonlinear
systems”. In Proceedings of AeroSense: The 11th International Symposium on
Aerospace/Defense Sensing, Simulation and Controls (pp. 182–193). Orlando, FL,
USA: SPIE. 1997.
Julier, S. J., Uhlmann, J., and Durrant-Whyte, H. F. “A new method for the nonlinear
transformation of means and covariances in filters and estimators”. IEEE Transactions
on Automatic Control, 45(3), 477–482. 2000.
Kailath T. and Frost P. “An Innovations Approach to least squares estimation- Part III:
Non Linear Estimation In White Gaussian Noise” IEEE Transactions on Automatic
Control, 1971.
Kalman, R.E and Bucy, R.S. “New Results in Linear Filtering and Prediction Theory”
Transactions of the ASME, Journal of basic engineering, March 1961.
262
Kamisnski, P. G. Bryson, Jr A.E. and Schmidt, S.F. “Discrete square root filtering: A
survey of current techniques”. IEEE Transactions on Automatic Control, December
1971.
Karlson R. “Particle Filters for Positioning and tracking applications” PhD thesis,
Linkoping Studies in Science and Technology. Dissertations. No. 924, Department of
Electrical Engineering Linkoping University, Linkoping, Sweden 2005.
Karlsson, R., Schön, T., and Gustafsson, F. (2004). Complexity analysis of the
marginalized particle filter. Technical Report LiTH-ISY-R-2611, Department of
Electrical Engineering, Linköping University.
Klaas M., de Freitas N. and Doucet A. “Toward Practical N2 Monte Carlo: the
Marginal Particle Filter”, 2005.
263
Koop R.E. and Orford R.J. “Linear Regresion Applied to system identification for
adaptive control systems” AIAA Journal, October 1963.
Kotecha, JH. and Djuric´ P. “Gaussian sum particle filtering” IEEE transactions on
signal processing, vol. 51, no. 10, october 2003
Kushner H., “Non Linear Filtering: The exact Dynamical Ecuations Satisfied by the
Conditional Mode” IEEE Transactions on Automatic Control, Vol. AC-12 No. 3, June
1967.
Kushner, H.J. and P.G. Dupuis (1992). NumericalMethods for Stochastic Control
Problemsin Continuous-Time. Springer-Verlag.
Leal Ascencio R. "Artificial Neural Networks as a Biomass Virtual Sensor for a Batch
Process., "Proceedings of the IEEE International Symposium on Intelligent Control,
September, Mexico City, (2001).
Leal R., Butler P, Lane P., Payne P.,”Data fusion and artificial neural networks for
biomass estimation,” Science Measurement and Technology. IEEE Procedings V. 144,
Issue 2, pp.69-72(1997)
Li, P., Goodall, R., and Kadirkamanathan, V. (2003). Parameter estimation of railway
vehicle dynamic model using Rao-Blackwellised particle filter. In Proceedings of the
European Control Conference, Cambridge, UK.
Lin, Y., Tanaka, S. Ethanol fermentation from biomass resources: current state and
prospects Appl Microbiol Biotechnol (2006) 69: 627–642 DOI 10.1007/s00253-005-
0229-x 2006
Liu, J. and Chen “Sequential Monte Carlo methods for dynamical systems” J. Amer.
Statistic. Assoc., vol. 93, pp. 1032–1044, 1998.
Ljung, L. System Identification Theory for the user. Prentice hall, Englewood cliffs,
Nj, 1987.
Lynd LR (1996) Overview and evaluation of fuel ethanol production from cellulosic
biomass: technology, economics, the environment, and policy. Annu Rev Energy
Environ 21:403–465
Maher, B., and Zeng, R., “Experimental Results and Model Reference Adaptive
Estimation and Control for a Fermentation Process”, Control Engineering Practice,
Vol 3, (1995).
Mauch, K., S. Arnold, and M. Reuss, “Dynamic Sensitivity Analysis for Metabolic
Systems,” Chem. Eng. Sci., 52, 2589–2598 (1997).
266
Mehra, R. K. “On the identification of variances and adaptive kalman filtering” IEEE
Transactions on Automatic Control, April 1970.
Menzel, K.; Zeng, A.; Biebl, H. e Deckwer, W. (2000), Kinetic, dynamic, and
pathway studies of glycerol metabolism by lebsiella pneumoniae in anaerobic
continuous culture: I. The phenomena and characterization of oscillation and
hysteresis, Biotechnology and Bioengineering, v. 52, n. 5, p. 549 – 560.
Meng X.-L. y van Dyk D. A., “Seeking efficient data augmentation schemes via
conditional and marginal augmentation” Biometrika, vol. 86, pp. 301–320, 1999.
Mitter, Sanjoy. Filtering and Stochastic Control: A historical perspective, MIT report
1996
Murthy K.P.N, “Monte Carlo: Basics. Theoretical Studies Section” Indira Gandhi
Centre for Atomic research, Indian Society for Radiation Physics, 2000.
Neˇsi´c, D. and A.R. Teel (2001). Perspectives in Robust Control. Chap. 14: Sampled-
Data Control of Nonlinear Systems: an Overview of Recent Results. Springer.
Neˇsi´c, D., A.R. Teel and P.V. Kokotovi´c (1999). Sufficient conditions for
stabilization of sampled-data nonlinear systems via discrete-time approximations.
Systems and Control Letters 38, 259–270.
Nielsen, J. and J. Villadsen, “Modelling of Microbial Kinetics,” Chem. Eng. Sci., 47,
4225–4270 (1992).
Nordlund, P.-J. (2002). Sequential Monte Carlo Filters and Integrated Navigation.
Licentiate Thesis No 945, Department of Electrical Engineering, Linköping
University, Sweden.
Ocone D. and Pardoux E., “A lie algebra criterion for the nonexistence of finite
dimensionally computable filters,” Lecture Notes Math., vol. 1390, pp. 197–204, 1988
Oliveira I., Rapôso A.C., Souto-Maior A.M., Lopes C.E. “Oscillatory behavior of
zymomonas mobilis in continuous stirred tank bioreactor” 2nd Mercosur Congress on
Chemical Engineering, Rio de Janeiro Brasil, (2005).
268
Onder, M., et al, Neural Network Assisted Nonlinear Controller for a Bioreactor.
Bogazici Univessity Estanbul, 1998.
Ozaki, T., Iino, M., “An Innovation approach to Non gausian time series analysis”. J.
Appl. Probab. 38A, 78-92, 2001.
Pan, Z. (2002). Canonical forms for stochastic nonlinear systems. Automatica 38,
1163–1170.
Parker, C.; Peekhaus, N.; Zhang, X. e Conway, T. (1997), Kinetics of sugar transport
and phosphorylation influence glucose and fructose cometabolism by Zymomonas
mobilis. Applied and Environmental Microbiology, v. 63, n. 9, p. 3519–3525.
Pitt, M. and Stephard N. “Filtering Via Simulation: Auxiliary particle Filter” Imperial
College Oxford Report, 1997.
Pitt, M and Shephard, N. “A fixed lag auxillary particle filter with deterministic
sampling rules,” unpublished paper, 1998.
269
Pons, M.-N., Bioprocess Monitoring and Control. Hanser Publishers, New York
(1992).
Poyiadjis, G., Doucet, A., Singh, S., “Particle Methods for Optimal filter derivative:
Application to parameter estimation” IEEE 2005.
Quintero, O., di Sciascio, F., and Álvarez, H., (2004). “Estructura de dos grados de
libertad y sensor virtual para controlar un bioproceso,” XIX Congreso Argentino de
Control Automático AADECA, Bs. As., Argentina.
Quintero, O., and di Sciascio, F. (2005). “Aplicación del filtro de Kalman extendido a
un proceso de fermentación alcohólica en continuo,” XI Reunión de Trabajo en
Procesamiento de la Información y Control RPIC.
Quintero, O., Scaglia. G.,. di Sciascio, F. (2008b) “Bio process control strategy based
on numerical methods and linear algebra: second approach” in IASTED conference on
“Modelling, Identification and Control” (MIC 2008).
Quintero, O., Scaglia. G.,. di Sciascio, F. (2008c) “Numerical Methods Based Strategy
and Particle Filter State Estimation For Bio Process Control” in IEEE- ICIT 2008,
April 2008 China.
Quintero, O., Scaglia. G., di Sciascio, F. (2008e) “Numerical method based controllers
applied to a biotechnological process”. Submitted on 26th November 2007, Revised
18th february, to Journal of Industrial Microbiology & Biotechnology.
270
Rallo R., Ferre-Gine J., Arenas A., Giralt F.,"Neural virtual sensor for the inferential
prediction of product quality from process variables” Computers and Chemical
Engineering 26,1735/1754, (2002).
Raposso Camêlo, A., Pinto, J., Pinheiro, I., “Parameter Estimation od Continuous
fermentations using Particle Swarm”. 2nd Mercosur Congress on Chemical
Engineering, Rio de Janeiro Brasil (2005).
Rogers, P, Jeon, Y, Lee, K, Lawford, H. Zymomonas mobilis for Fuel Etanol and
Higher Value Products Adv Biochem Engin/Biotechnol (2007) 108: 263–288 DOI
10.1007/10_2007_060 © Springer-Verlag Berlin Heidelberg Published online: 24
May 2007
Romanenko A., Castro J., “The unscented filter as an alternative to the EKF for
nonlinear state estimation: a simulation case study”. Computers and Chemical
engineering 2003.
Rossi Vivien and Vila Jean-Pierre Filtering discrete time nonlinear system with
unknown parameters: a non parametric approach Preprint submitted to Elsevier
Science 3rd February 2005
Rubin, Multiple Imputation for Nonresponse in Surveys, New York: Wiley, 1987.
271
Russell S. A., Robertson D. G., Lee J. H. and. Ogunnaike B. A. “Model -based quality
monitoring of batch and semi-batch processes,” Journal of Process Control, V. 10,
Issue 4, Pages 317-332, (2000)
Scaglia,G., Mut .V, Rosales A., Quintero O. “Tracking Control of a Mobile Robot
using LinearInterpolation“ IMAACA,07 vol. 1, pp. 11-15, ISBN: 978-2-9520712-7-7
Schon, Thomas, On computational Methods for Non Linear Estimation – Lic Thesis,
Division of Automatic Control and Communication Systems Department of Electrical
Engineering Linkoping University, SE581 83 Linkoping, Sweden, 2003.
272
Schön, T., Gustafsson, F., and Hansson, A. (2003b). A note on state estimation as a
convex optimization problem. In Proceedings of the IEEE International Conference on
Acoustics, Speech, and Signal Processing, volume 6, pages 61–64, Hong Kong.
Schön, T. B., Eidehall, A., and Gustafsson, F. “Lane departure detection for improved
road geometry estimation”. Technical Report LiTH-ISY-R-2714, Department of
Electrical Engineering, Linköping University, Sweden. Submitted to the IEEE
Intelligent Vehicle Symposium, Tokyo, Japan. 2005.
Schön, T. B., Karlsson, R., and Gustafsson, F. (2006a). The marginalized particle
filter in practice. In Proceedings of IEEE Aerospace Conference, Big Sky, MT, USA.
Invited paper, accepted for publication.
Schön, T. B., Wills, A., and Ninness, B. (2006b). Maximum likelihood nonlinear
system estimation. In Proceedings of the 14th IFAC Symposium on System
Identification, Newcastle, Australia. Accepted for publication.
Schön, T., Gustafsson, F., and Nordlund, P.-J. (2005). Marginalized particle filters for
mixed linear/nonlinear state-space models. IEEE Transactions on Signal Processing,
53(7):2279–2289.
Sing, S., Kantas, N., Vo, B-N., Doucet, A., “Simulation Based Optimal sensor
scheduling with application to observer trajectory planning” Submitted to Automatica,
29 March 2005.
Stoica P and Jansson M. MIMO system Identification state space and subspace
approximations versus transfer function and instrumental variables. Subbitted to IEEE
trans. Sig. Proc - Royal Institute of Technology, 2000.
Strassle, C., B. Sonnleitner, and A. Fiechter, “A Predictive Model for the Spontaneous
Synchronization of Saccharomyces cerevisiae Grown in Continuous Culture. II.
Experimental Verification,” J. Biotech., 9, 191–208 (1989).
Sun M. and Glowinski R., “Pathwise approximation and simulation for the Zakai
filtering equation through operator splitting” Calcolo., vol. 30, pp. 219–239, 1993.
Tanizaki. H. “Nonlinear and non normal filter using importance sampling: antithetic
monte carlo integration” 1991.
Tanizaki, H. and Mariano R., “Non-Linear and Non-Normal Filter Based on Monte-
Carlo Technique” 1995a, 1995b
Tanizaki, H. and Mariano R., “Nonlinear filters based on taylor series expansions”
1994.
Tano, M., Buzato, J., Celliogoi, M.A., “Sugar Cane Juice Fermentation by
Zymomonas Mobilis CP4 Subjected to Inhibition by ethanol and High Initial
Concentratoin of Substrate”. Brasilian Arch. Biol. Tech, Vol43, April 2000.
Thrun, S. Fox “Particle filters for mobile robot localization” in Sequential Monte
Carlo Methods in Practice, A. Doucet, N. de Freitas, and N. Gordon, Eds. New York:
Springer Verlag, 2001.
Twardoska K., Marnik, T., Paslawska, M., “Approximations of the Zakai equation in a
nonlinear filtering problem with delay” Int. J. Appl. Math. Comput. Sci, 2003 Vol 13,
No 2.
Vo, B., Singh, S. and Doucet, A. (2005). Sequential Monte Carlo methods for multi-
target _ltering with random _nite sets. IEEE Transactions on Aerospace and
Electronic Systems. To Appear.
http://www2.ee.unimelb.edu.au/staff/bv/vo/phdfilterAESrevfinal2cola.pdf.
275
Wan, E van der Merwe, A., Rudolph A T. “Dual Estimation and the Unscented
Transformation” Oregon Graduate Institute of Science & Technology Department of
Electrical and Computer Engineering.
Wang, X., R. Chen, R., and Guo, D. “ Delayed-pilot sampling for mixture Kalman
filter with application in fading channels”. IEEE Transactions on Signal Processing,
50(2):241–254. 2002.
Wang, F., Bahri, P., Lee, P., Cameron, I., A multiple Model, State feedback strategy
for robust Control of Non linear Process. European Symposium on Computer arded
Process engineering-15, Elsevier, 2005 .
Yan Lin . Shuzo Tanaka. Ethanol fermentation from biomass resources: current state
and prospects. Appl Microbiol Biotechnol (2006) 69: 627–642 DOI 10.1007/s00253-
005-0229-x.
Yuz, J.I. (2005). Sampled-data models for linear and nonlinear systems. PhD thesis.
School of Electrical Eng. and Computer Science, The University of Newcastle.
Australia.
Yuz, J.I. and G.C. Goodwin (2005a). Generalised ¯lters and stochastic sampling zeros.
In: Joint CDC-ECC'05. Seville, Spain.
Yuz, J.I. and G.C. Goodwin (2005b). On sampleddata models for nonlinear systems.
IEEE Trans. on Automatic Control (to appear).
Yuz, J.I., G.C. Goodwin and H. Garnier (2004). Generalised hold functions for fast
sampling rates. In: Proc. of the 43rd IEEE Conference on Decision and Control.
Nassau, Bahamas.
Yuz, J and Goldwin G. Sampled data models for stochastic non linear systems. 2006
Zhang J., Liu, J., “A new sequential importance sampling method and its applications
to the two dimensional hydrophobic-hydrophilic model”. Journal of chemical physics
volume 117, number 7, august 2002.
Zhou K. with Doyle J. and Glover K. Robust and optimal control, Prentice Hall, NJ
APÉNDICE A
PRELIMINARES MATEMÁTICOS
Sigma Álgebras
álgebra :
Sea S un conjunto y F una familia de subconjuntos de S. F es una sigma álgebra si
F
A F Ac F
A1 , A2 ,... F U i 1 Ai F
Espacios de Probabilidad
Función de Probabilidad
dP( x)
p( x)
d
Nx
Con x X y p ( x)
278
Proceso estocástico
T ( t ,w ) x( t , w )
Tal que para tiempo fijo t = t0, x(t0,.) es una variable aleatoria.
Además, para cada w fijo, w = w0, obtenemos una realización (path) específica que es
función del tiempo.
Filtración
Proceso de Markov
Martingala
279
1. 0 0
2. t1 , t2 t1 ,..., tk tk 1 son independientes para todo
t1 t2 ... tk 1 tk
3. t s N 0, t s 0 s t
4. La trayectoria simple definida por t t; w es continua para todo w
T
Un vector n dimensional de procesos t 1 t ... n t , donde casa proceso
escalar t es un movimiento Browniano estándar e independiente se llama un
1
Ruido blanco:
Se define como ruido blanco al proceso estocástico ( t ) que verifica las siguientes
propiedades
( t ) es gaussiano para todo t, Tiene media nula, E{ ( t )} 0
Función de covarianza:
T
r( t , s ) E{ ( t ) ( s )} (t s ) t ,s : s t
T
r( ) E{ ( t ) ( t )} ( ), t s
t
W(t ) ( )d donde W(t) es un proceso de Wiener.
0
Procesos estacionarios:
Ergodicidad:
Integral de Itô
1. t , w f t , w esB 0, Fmedible
2. Existe una filtración Ft , de modo que t sea una martingala con respecto a Ft y
f t , w este Ft adaptado.
T
2
3. E f t , w dt
s
dx t f x, t dt g x, t d t
n n n d
donde f : 0, es la función de base o drift, g : 0, es la
matriz de dispersión y t es el movimiento Browniano con matriz de difusión
Qc t . La matriz definida por g x, t Qc t g T x, t es la matriz de difusión de la
ecuacion diferencial estocástica.
La ecuacion diferencial estocástica puede ser expresada de la siguiente manera:
t t
x t x s f x, t dt g x, t d t
s s
Integral de Stratonovich
dx t f x, t dt g x, t d t
Formula de Itô
dx t f x, t dt g x, t d t
Sea g una Ecuacion doblemente diferenciable. Entonces las componentes del proceso
estocástico
y(t ) g x(t ), t , satisfacen la ecuaciones diferenciale estocástica
2
gk gk 1 gk
dyk dt dxi dxi dx j
t i xi 2 ij xi x j
donde los términos dxi dx j se calculan de acuerdo a las reglas
dtd 0; d dt 0; d d T Qc (t )dt
en caso que la ecuacion fuera en el sentido de Stratonovich, el diferencial seria
gk gk
dyk dt dxi
t i xi
que es igual al resultado del calculo ordinario.
Distribución empírica
Np
1
Pˆ ( x) ( x xi )
Np i 1
1 si x xi
x X es discreta, la función ( x xi ) ,
0 si x xi
dPˆ ( x) pˆ ( x)dx 1
X X
Teorema de Girsanov
n
Asuma que t :0 t T es un proceso F medible adaptado a la filtración
n
natural Ft F de un movimiento Browniano n dimensional t :0 t T
con respecto a la medida P.
Si
t
2
E exp t dt
0
entonces
t t
T 1 2
Z t exp td t t dt
0
20
satisface la ecuacion
t
T
Z t 1 Z t td t
0
En este tipo de estimación, desde el punto de vista de Fisher, también llamado como
paramétrico, los parámetros del modelo o los estados, no son considerados como
variables aleatorias, pero como se menciono anteriormente se consideran fijas pero no
conocidas y el valor del parámetro afecta explícitamente las propiedades estadísticas
de las observaciones de manera conocida. La conexión entre las medidas y los
parámetros está definida por la función de densidad de el vector de observación y,
parametrizado por el vector de parámetros x, de la forma p(y/ x), la cual es función de
los parámetros después de haber hecho las medidas. Esta función de vecindad se
escribe de la siguiente forma combinando las vecindades para diferentes instantes de
tiempo, asumiendo independencia de eventos:
l ( x y) p( y x)
Nótese que para enfatizar, dicha vecindad es función del primer argumento y luego, se
consideran las medidas u observaciones y. Luego de considerar dichas medidas, el
vector de parámetros x corresponde a la mejor densidad elegida a partir de las medidas
observadas.
Fisher plantea que una buena estimación del vector de parámetros o del vector de
estados, esta dado por el argumento que maximiza la función de vecindad; la teoría se
basa en que asintóticamente, la estimación de máxima verosimilitud converge casi
seguramente al valor verdadero, bajo las condiciones generales:
x ML arg max l ( x y )
x
Lo que da pie a la teoría de estimación de máxima verosimilitud o de vecindad. Para
mayor detalle ver (E.L. Lehmann. 1991).
Puede verse que el posteriori debe ser proporcional a la vecindad y en eso, las
aproximaciones Bayesiana y Fisheriana coinciden.
Algunos argumentan que la necesidad de conocimiento previo para definir p(x), hace
que el punto de vista Bayesiano sea menos atractivo, lo que favorecería la
aproximación paramétrica donde no es necesario tanto conocimiento previo. Sin
287
embargo, si este conocimiento previo existe, es difícil de describir, por lo menos como
función de densidad de probabilidad; un razonamiento similar afirma que dado que la
elección de la función previa es algo subjetiva, puede afectar la inferencia en forma
destructiva. Estos argumentos se encuentran a favor de la estadística paramétrica,
frente a la Bayesiana, mas allá de esto, el previo y la vecindad define el problema de
inferencia como se presentó anteriormente y definir una función de vecindad correcta,
refleja el mismo problema que definir una buena función de densidad previa.
T
J (u ) E x 't Qxt u 't Rut dt x 't Mxt , Q 0, R 0, simetricas
0
288
Cuando las perturbaciones estocásticas están presentes hay una diferencia fundamental
entre el lazo de control abierto (donde el control no es función de las observaciones) y
el lazo de realimentación (donde la acción de control es función de las observaciones).
La respuesta para resolver este problema estocástico totalmente observable
vt 0, t 0, T H I
Considerando a
En [Mitter, 1996] puede verse que la acción de control óptima esta dada por
u o (t ) R 1 B ' P (t ) xt
T
J (uo ) tr G ' P( s)G ds m0' P(0)m0 tr 0 P(0)
0
E ( x0 ) m0
cov( x0 ) 0
Pu ( xt yt ,0 s t ) ,
Con v siendo un movimiento Browniano adaptado a Fyt, el sigma álgebra generado por
289
ys 0, 0 s t .
T T
J ( n) E xˆ t' Qt xˆt u 't Rut dt xˆ t' Mxˆt tr P (t )Q dt tr P (T )M
0 0
Con
u o (t ) R 1 B ' (t ) xˆt