martes, 24 de enero de 2012

ARX - Algoritmo propio para identificación del sistema

El algoritmo provisto por Matlab demostró funcionar a la perfección con las series de ejemplo pero al intentar utilizarlo con volúmenes de datos más grandes los tiempos de cómputo crecen exponencialmente. Como agregado, el sistema final debe programarse en C/C++/CUDA donde la utilización del toolkit de Matlab es poco práctica.
En consecuencia se decide implementar la identificación del modelo ARX. A tal fin se analizan ciertas publicaciones de Lennart Ljung (autor del toolkit de Matlab) y se consigue un algoritmo que entrega resultados similares en menor tiempo.
Para comprobar la bondad del algoritmo creado se realizan comparaciones contra los modelos obtenidos con el toolkit de Matlab para las series 'Motorized Camera' y 'Missile' previamente estudiadas.
Se utiliza la función compare del toolkit para verificar el ajuste con muy buenos resultados. Las siguientes imágenes muestran la comparación gráfica entre los resultados obtenidos por Matlab (llamado model) y por el algoritmo basado en los escritos de Ljung (llamado modelLjung):
Serie 'Motorized Camera'. Comparación entre algoritmos de identificación.





Serie 'Missile'. Comparación entre algoritmos de identificación.

lunes, 16 de enero de 2012

ARX - Identificación del Sistema

El modelo físico a estudiar debe estar representado en el software de control de forma tal que sea posible predecir los futuros estados y determinar cual es la mejor actuación para poder conseguir el resultado deseado.
Existen varias maneras de modelar numéricamente una experiencia física. En este caso optamos por un modelo de caja negra llamado Auto Regressive with eXogenous input (ARX). En este modelo, es necesario definir grandes matrices de coeficientes que permitan, dado un vector de estado con información de pasos previos, obtener una predicción a futuro del sistema. Como todo modelo de caja cerrada, los coeficientes necesarios son definidos 'por inspección' o con métodos de identificación en base a mediciones previas.

La ecuación que define el modelo ARX más simple, una entrada y una salida (SISO), es la siguiente:


donde:
y: variable de salida del sistema
u: variable de entrada del sistema
a y b: coeficientes asociados

En caso de tener un sistema con entradas y salidas múltiples (MIMO), las variables u e y se transforman en vectores (un componente de u por cada entrada, un componente de y por cada salida). Sin cambios en la ecuación, se la puede extender a la forma matricial para representar ese esquema.
Por último, los coeficientes a y b se pueden agrupar utilizando una notación reducida:


Donde se definen las matrices de polinomios A y B según:



Siendo Z el operador de corrimiento:






En todos los casos se debe entender la existencia de un término de error en cualquiera de los lados de la igualdad que corrige las incertezas que provoca el modelo a caja cerrada. Con esto en mente, se puede entender el modelo con el siguiente esquema:




Matlab posee un toolkit de identificación de sistemas que entrega el modelo ARX computado a partir de una serie de muestras. Para probar el toolkit utilizamos series predefinidas en Matlab llamadas 'Motorized Camera' y 'Missile' que incluyen datos reales de varios sensores.
Se puede obtener más información sobre el algoritmo y las funciones utilizadas desde aquí.

% REAL DATA 1: MISSILE
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'missiledata.mat'));
% o bien: 
% load motorizedcamera
data = iddata(y, u);

p = 2;
outputs = size(y,2);
inputs = size(u,2);
model = arx(data, [p*ones(outputs, outputs), p*ones(outputs, inputs), ones(outputs, inputs)]);

viernes, 2 de diciembre de 2011

Video de la Simulación con Actuadores Programados

Gracias al código Fortran agregado a los loops de simulación del CodeSaturne, pudimos agregar actuadores de plasma a los resultados anteriores.
De esta forma simulamos la actuación en dos sectores del cilindro, realizando pulsaciones de plasma con ángulo de salida tangente para ver la respuesta del sistema. Se espera que las turbulencias del cilindro desaparezcan con el plasma aplicado y vuelvan a generarse tan pronto se retire el efecto del actuador.
Durante esta simulación se agregaron algunos esquemas de actuación periódicos, oscilantes y asiméticos para poder relevar información que será usada en la identificación del sistema mediante el modelo ARX.


domingo, 13 de noviembre de 2011

Usando variables 'static' en Fortran para cargar un CSV sólo una vez

Para completar la idea de utilizar un vector desde un archivo dentro del código de code saturne es importante dar ciertas garantías de performance. La subrutina a utilizar se encuentra dentro del archivo usclim.f90 que es utilizado de forma repetitiva durante la simulación.

Una forma de garantizar que la carga del CSV se realiza una sóla vez es utilizando código global y variables globales. Lamentablemente, el archivo disponible no cuenta con las secciones comunes con lo que se supone una composición de ese código dentro del programa de simulación general.

Otra alternativa es utilizar variables en el data segment persistentes entre llamadas a la subrutina. Optamos por ese camino y utilizamos la palabra reservada SAVE para realizar un breve programa de prueba:

program read_csv_once
  implicit none
 
  call readAndPrint(1000)
  call readAndPrint(100)

end program read_csv_once

subroutine readAndPrint(test_offset)
  integer:: test_offset
  integer, parameter:: values_qty = 2
  real, save, pointer:: pValues(:)
  logical, save:: pValues_initialized= .false.

  if (.not.pValues_initialized) then
    allocate(pValues(values_qty))
    open(unit=99, file="act1.csv", action="read")
    read(99, *) pValues
    close(99)
    pValues = pValues + test_offset
    pValues_initialized = .true.
  endif

  write(*,*) "Values from file:", pValues
  write(*,*) "Values Quantity:", values_qty
end subroutine


Nuevamente tenemos una cantidad fija de valores a leer del CSV dado que por el momento es un número conocido.
Notar que el puntero estático pValues NO posee un valor por defecto definido, es por ese motivo que se utiliza pValues_initialized para controlar si se hizo la lectura del CSV o no.
El resultado de correr este programa de pruebas es:

 Values from file:   1003.6295       1003.6688    
 Values Quantity:           2
 Values from file:   1003.6295       1003.6688    
 Values Quantity:           2

Demostrando que la inicialización se realizó una vez.


Algo en fortran: leyendo un array desde CSV

Ante la preguntá de ¿qué motivos tengo para hacer eso? contesto: necesito levantar un array dinámico desde un archivo para hacer uso de los valores dentro de una simulación con code saturne.

Vamos a un código muy sencillo:

program readcsv
  implicit none
 
  call readAndPrint()
end program readcsv

subroutine readAndPrint()
  integer, parameter:: values_qty = 20
  real, pointer:: pValues(:)

  allocate(pValues(values_qty))
  open(unit=99, file="act1.csv", action="read")
  read(99, *) pValues
  close(99)

  write(*,*) "Values from file:", pValues
  write(*,*) "Values Quantity:", values_qty
end subroutine

La cantidad de valores a leer está fija por el momento. Una posible solución sería hacer una prelectura de la primer fila y contar la cantidad de caracteres ',' para luego hacer una alocación con la cantidad exacta de valores. Otra posibilidad es leer toda la fila y agregarlos al array realocando en cada iteración.

martes, 1 de noviembre de 2011

Primeros Resultados de la Simulación

Luego de muchas corridas fallidas, donde el mallado de la experiencia y los distintos parámetros numéricos fueron ajustados, llegamos a una corrida limpia de la simulación.

El CodeSaturne demostró ser un tanto complejo en su puesta a punto pero muy potente y versátil. Permite cálculos simples como la experiencia del cilindro o mucho más complejos usando nociones químicas o de termodinámica.

En este caso nos limitamos a mostrar unas imágenes del cilindro luego de alcanzar un régimen estacionario respecto de las oscilaciones de la estela. En la imagen se puede notar el desprendimiento de vórtices que se pretende controlar en el sistema de control y para el cual será necesario agregar unos actuadores al mallado del cilindro.




martes, 25 de octubre de 2011

Mallado para Simulación de la Experiencia

A fin de conseguir una simulación realista de la experiencia mediante el uso de Code Saturne armamos un mallado de la experiencia que nos permite definir la simetría en juego. Existen varias herramientas de mallado, optamos por GiD por la pequeña curva de aprendizaje y la disponibilidad de información dentro del laboratorio. Como agregado este software tiene la posibilidad de hacer simulaciones con y la posibilidad de hacer simulaciones con módulos codificados por colaboradores del laboratorio.
En este caso no utilizaremos GiD para simular o realizar post-procesamiento sino simplemente para el mallado de la experiencia.
La simetría del experimento es muy sencilla y los mallados preliminares no fueron problemáticos. A continuación se incluyen distintos acercamientos de la experiencia, donde se puede observar el cilindro sobre el lado izquierdo dejando buena parte del mallado derecho para la estela que debe producir.







A lo largo de las imágenes se puede observar que el grado de precisión necesario para los nodos del mallado implica una cantidad enorme de nodos. Por otro lado se intenta minimizar la precisión en zonas donde no se requiere gran nitidez durante el procesamiento o el análisis de resultados.