Sistema de Control - Error con Suerte

Algoritmo de Control

Ya identificado el modelo se plantea un algoritmo de control generalizado (o Generalized Predictive Control - GPC) donde se plantea un funcional a minimizar y se itera desde un punto de inicio hacia una solución óptima tomando un paso del gradiente en cada iteración.  Trataremos de dar una explicación más detallada en otro post.
En este caso, la función a minimizar es la norma de los valores medidos de salida (los sensores colocados a cada lado del cilindro) más un costo por la utilización del actuador de plasma. El resultado esperado es minimizar la magnitud medida en los sensores, consiguiendo que los trazadores escurran por el eje central, utilizando la menor cantidad de energía posible en los actuadores.

Primera Implementación

La primer implementación del sistema requirió de las técnicas de comunicación ya descriptas entre Fortran y C++. En efecto, fue necesario detener el algoritmo de Code Saturne (Fortran) para informar y esperar un valor de actuación del algoritmo controlador (C++).
Como agregado, dados los cambios y correcciones frecuentes en la minimización del funcional y la identificación ARX del modelo, parte del código se realizó en Matlab por lo que se empleó MEX para su ejecución.

Un Error en el Código con Resultados Prometedores

Aún luego de emplear varias horas en revisar los cálculos y comunicación entre procesos, no se pudo evitar que se filtrara un error de programación.
En este caso el error fue tan simple como iniciar a cero la variable de acumulación para cada iteración. El efecto inmediato fue que el algoritmo no avanzara hacia una solución del funcional sino que mantuviera una relación de proporción con la magnitud de los trazadores. Por fortuna, esta relación también provee un esquema de control que, aunque rudimentario, resuelve el objetivo final. De esta forma el valor del actuador aumenta si la norma de los trazadores aumenta y disminuye si así lo hacen los trazadores.

Resultados de la Simulación

El resultado final es un actuador oscilante con frecuencia similar al experimento. A continuación se muestran los valores de actuación computados por el algoritmo de control. En este caso se plantea un control con inicio en el paso 6000 de la simulación, luego de entrar en régimen estacionario.

El algoritmo fuerza una actuación discontínua al inicio pero luego retoma un función de oscilación suave y con la frecuencia del experimento.
Los valores medidos en los sensores se pueden observar en el siguiente gráfico:

Luego de la actuación, y omitiendo el caso de discontinuidad, se puede observar una gran reducción en el valor medido, producto del soplido de plasma.

El siguiente video ilustra los resultados previos:

Comunicación C++ - Matlab

Otro punto de contacto entre las plataformas fue encontrado en el algoritmo de control a realizar por lo que fue necesario investigar las librerías MEX.

El escenario es sencillo, se posee código C/C++ que debe comunicarse con Matlab enviando matrices y efectuando procesamiento con funciones de toolkits conocidos. Una vez concluido el procesamiento, se debe recuperar los valores de las matrices resultantes.

Para ilustrar este tipo de comunicación dejamos un ejemplo simple sobre cómo se efectúa una operación básica en Matlab desde C/C++.

#include <string.h>
#include <stdio.h>
#include "engine.h"
#define BUFSIZE 1024

int main(int argc, char* argv[]) {
        Engine* ep;
        if (!(ep = engOpen("\0"))) {
                fprintf(stderr, "\nCan't start MATLAB engine\n");
                return 1;
        }
        else {
                char buffer[BUFSIZE+1];
                memset(buffer, 0, sizeof(buffer));
                engOutputBuffer(ep, buffer, sizeof(buffer) - 1);
                engEvalString(ep, "pi^2");
                printf("Evaluating: pi^2. \nResult: %s\n", buffer);
                engClose(ep);
                return 0;
        }
}

 Existen ciertas compilaciones para compilar y ejecutar el código anterior. La manera más fácil es sin duda realizarlo por la consola de Matlab:

>> mex('-v', '-f', fullfile(matlabroot,'bin','engopts.sh'),'main.cpp')
 -> g++-4.0 -c  -I/Applications/MATLAB_R2010a.app/extern/include[...] -DMACI64  -DMX_COMPAT_32 -O2 -DNDEBUG  "main.cpp"
 -> gcc-4.0 -O -Wl,-twolevel_namespace -undefined error -arch x86_64[...] -mmacosx-version-min=10.5 -o  "main"  main.o  -L/Applications/MATLAB_R2010a.app/bin/maci64 -leng -lmx -lstdc++

>> !./main
Evaluating: pi^2.
Result: >>
ans =

    9.8696

Para realizar la compilación por línea de comando, notar el uso de la variable matlabroot. En su directorio bin se puede encontrar el comando mex y la librería de configuración engopts.sh.

$ /Applications/MATLAB_R2010a.app/bin/mex -v -f /Applications/MATLAB_R2010a.app/bin/engopts.sh main.cpp -o main

Si bien la compilación no se realiza con el comando gcc sino con mex, la ejecución sólo requiere que Matlab y csh estén instalados sin cambiar la línea de comando.

Comunicación Fortran - C++

Para realizar una simulación realista del efecto que provoca el control sobre el modelo físico es necesario aplicar la actuación en simultaneo con la resolución del sistema. Para ello se plantea un esquema de comunicación entre el Code Saturne y nuestro programa en C++ donde podamos definir el valor de actuación para cada paso de la simulación.
Se cuenta con varias posibilidades concretas de comunicación IPC: sockets, pipes, archivos. En este caso se opta por la comunicación mediante archivos FIFO (es decir, implementando el concepto de pipes mediante archivos virtuales) donde un programa lee del archivo mientras el otro escribe en el mismo.
Para analizar la viabilidad se decide realizar un programa simple donde dos procesos se comuniquen mutuamente. En este caso contamos con un proceso que lee e imprime los valores de actuación a utilizar (hecho en Fortran 90) y otro proceso que genera y escribe los mismos (C++).

Lector (Fortran 90) - lector.f90

program readcsv
  implicit none
  call readAndPrint()
end program readcsv

subroutine readAndPrint()
  integer status
  double precision value

  open(unit=99, file="act1.csv", action="read")
  do
    read(99, *, IOSTAT=status) value
    if (status > 0) then
      !Error
      write (102,*) 'An error has been detected in the file'
      exit
    else if (status < 0) then
      !EOF
      exit
    else
      print *, 'Value:',value
    end if
  end do
  close(99)
end subroutine

Escritor (C++) - escritor.cpp

#include <iostream>
#include <fstream>
#include <ctime>

int main(int argc, char* argv[]) {
        std::ofstream arch("act1.csv", std::fstream::app);
        for (int i = 0; i < 2; ++i) {
                arch << 1.0 + i / 100.0 << std::endl;
                sleep(2);
        }
        return 0;
} 

Para ejecutar el programa es necesario abrir un archivo FIFO del sistema operativo. Esto es muy sencillo en sistemas Linux/Unix:

mkfifo act1.csv

La secuencia real de compilación y ejecución es la siguiente:

$ gfortran -o lector lector.f90
$ g++ -o escritor escritor.cpp
$ mkfifo act1.csv
$ ./lector &
[1] 892 

$ ./escritor

 Value:   1.0000000000000000     
 Value:   1.0100000000000000     
[1]+  Done                    ./lector

Reducción de Variables y agregado de Trazadores

Sensores

Luego de las primeras pruebas con el algoritmo ARX propio del laboratorio, y teniendo en cuenta los pobres resultados de performance y precisión para grandes cantidades de variables, se decide simplificar el problema de estudio a un conjunto reducido de sensores.
Como primer aproximación se plantea la utilización de 4 únicos puntos sensores y la medición de la presión en esos puntos.

Los puntos se encuentran simétricos al eje X por lo que se espera una función espejada en cada par de los mismos. Los puntos colocados en la segunda línea de sensores deberían mantener cierta similitud con los de la primer línea pero con un decaimiento en la intensidad debido a la difusión más un retardo en el tiempo por la velocidad del escurrimiento.

Trazadores

Si bien esta información se puede obtener en el simulador CodeSaturne, es necesario idear un método de procesamiento que permita calcularla de las imágenes a obtener en el experimento real. Por tal motivo se decide agregar trazadores a la simulación que serán implementados con fuentes de humo (u otras sustancias que queden en suspención sobre el fluido) al montar el experimento.
A tal fin se agregan escalares a la simulación del Code Saturne como se observa a continuación:

ARX - Resultados del Algoritmo con Información Real

Luego de escribir el algoritmo ARX que soporte grandes vólumenes de información se decide ponerlo en práctica con los datos experimentales obtenidos de la simulación del Code Saturne.
De esta forma se pretende obtener resultados más cercanos a los esperados de la experiencia real y comprobar que los algoritmos escritos hasta el momento sean efectivos y eficaces.
Se practican varias pruebas con los valores de vorticidad obtenidos de la simulación obtenidos mediante una grilla definida en la zona de la estela del cilindro. Se colocan varias líneas de obtención de información y por cada línea se relevan aprox. 50 puntos en 500 pasos de tiempo a lo largo de 20 líneas.

Los resultados obtenidos con el total de la información no fueron satisfactorios. Las simulaciones divergen y se demoran mucho tiempo.
A fin de determinar un buen número de muestras a utilizar se realizan distintas comparaciones planteando configuraciones de: cantidad de líneas utilizadas, cantidad de muestras por línea, cantidad de pasos de simulación a predecir (antes de poder relevar los valores reales y poder realimentar al predictor). Los resultados fueron positivos siempre que se utilizaran menos de 200 variables de entrada. Se puede ver el porcentaje de ajuste para las distintas configuraciones en el siguiente gráfico:

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.

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)]);

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.

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.

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.

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.

Code Saturne: Realizando una simulación numérica del experimento

Ya con varias piezas dispuestas para la experiencia del control del cilindro, nos es necesario datos más cercanos a la realidad para probar los módulos por separado. Si bien se podría plantear el experimento y obtener imágenes del mismo, resulta más simple (y flexible) realizar una simulación numérica. A tal efecto, elegimos la librería de código abierto Code Saturne.


Dado que la instalación es un tanto tediosa y se lleva a cabo frecuentemente en el laboratorio, creamos una wiki para el laboratorio y tomamos nota de los pasos utilizando el archivo por lotes: install_saturne.py. Para más información, ver aquí.

La ejecución no deja de ser algo compleja:

1. Se debe definir la experiencia con su simetría y parámetros físicos
2. Se debe utilizar un programa de mallado para generar la cuadrícula para simular la física elegida. En este caso se utiliza el programa GID.
3. Se debe correr el sistema de simulación Code Saturne en alguna de sus variantes
1. Mediante el entorno gráfico.
2. Realizando la configuración por archivos y corriendo el sistema por línea de consola. Esta opción es muy común para poder correr los trabajos de forma remota. Se agrega un instructivo a la flamante wiki del laboratorio. Para más información, ver aquí.

En el proceso nos topamos con un blog muy interesante de un estudiante de fluidodinámica del Reino Unido que vale la pena visitar:
http://code-saturne.blogspot.com/

Probando el código C/CUDA (o C/C++) usando Matlab

Dado que comunicar Matlab con el código C/CUDA (o C/C++) no es la principal preocupación para la tesis, se aplica un poca investigación al tema y se implementa rápidamente un protocolo que permite comunicar entradas y salidas. En este caso, la comunicación consta de varias matrices de entrada y varias de salida. En todos los casos se trata de números flotantes.
Para que la ejecución del programa escrito en C sea transparente, se crea una función Matlab con el mismo nombre que el programa C, esta función será la encargada de grabar en disco las variables Matlab, llamar al ejecutable escrito en C, esperar que concluya exitosamente y leer las variables de salida desde el disco.
Las variables en disco se crean en una ubicación definida por el código Matlab (generalmente en una carpeta temporal) y se comunican al programa C mediante argumentos

Este esquema tiene un payload asociado a la lectura/escritura en disco de las variables que podría evitarse con algún método más directo...
A futuro, se recomienda utilizar la compilación MEX de Matlab:
http://www.mathworks.com/help/techdoc/matlab_external/f23224.html

Librerías sobre Operaciones matemáticas en C/C++

A la hora de realizar los programas de prueba en C/C++, previo a codificarlos en C/CUDA fue necesario definir ciertas funciones de soporte matemáticas.

Como todo sistema de operaciones básicas comenzó con una clase Matrix<T> propia, con operaciones simples como multiplicar o sumar. A medida que el sistema creció, fue necesario implementar transposiciones y vistas de columnas o filas de la matriz. En ese punto, y buscando una implementación de inversión de matrices, resultaba obvia la necesidad de buscar un mejor soporte en librerías existentes.
A continuación, una colección de links relacionados:

• C/C++
• Gnu Scientific Library, con muchas operaciones y TDAs escritos en C
• http://www.gnu.org/software/gsl/
• Projectos VS2010 para compilar GSL en Windows: una contribución de Brian Gladman donde se descomprimen los proyectos dentro de la carpeta de descarga del GSL y con una lectura del readme se obtienen las librerías
• http://gladman.plushost.co.uk/oldsite/computing/gnu_scientific_library.php
• Clase Matrix con operaciones de inversión de Mike Dinolfo
• http://users.erols.com/mdinolfo/matrix.htm
• Numerical Recipes, se trata de librerías matemáticas muy usadas y conocidas. Son pagas pero con un bajo costo
• http://www.nr.com 
• Matlab
• Utilización de código Matlab desde C/C++
• http://www.mathworks.com/help/techdoc/matlab_external/f38569.html
• Pequeño instructivo para utilizar GSL en Matlab
• http://www.kellogg.northwestern.edu/faculty/papanikolaou/htm/gsl.htm
• Python 
• Plataforma para python
• http://coan.burnham.org/downloads
• SAGE, otra librería para python
• http://www.sagenb.org/

Luego de evaluar todas las posibilidades, optamos por los dos primeros links: utilizar GSL y realizar una compilación custom para Windows.

Simulando modelos conocidos

Matlab provee varias funciones que son de utilidad para simular modelos del tipo SISO, MIMO, etc. Una de las más útiles es sim que, dado un conjunto de datos de entrada y/o salida más un modelo ya definido, aplica el modelo a los datos y retorna la salida ulterior.

Veamos un ejemplo:

%genero un modelo según y(n) = u(n) - 2u(n-1)

A = [1] ; B = [1 -2];

modelo = idpoly(A, B);

%genero los datos para la simulación:

% ejemplos:

% iddata([], [1:100]') para salida nula y entrada del 1 al 100

% iddata([], idinput(100)) para salida nula y entrada aleatoria -1;1

data = iddata([],[1:10]');

%simulo

simulacion = sim(modelo, data);

plot(simulacion.y);

Bastante sencillo, si se cuenta con los datos de entrada e información sobre el modelo podemos calcular rápidamente la salida. 

Obviamente, se pueden aplicar relaciones mucho más complejas cambiando los coeficientes de A y B. Tanto es así que son esos vectores los que definen la transferencia G(s) que relaciona Y(s) y U(s) (salida con entrada). Entonces, el numerador evaluado en el operador de regresión q se convierte en A(q) mientras que el denominador hace lo propio en B(q). Si observamos con atención, podemos ver esa información en la variable modelo:

>> modelo

Discrete-time IDPOLY model: y(t) = B(q)u(t) + e(t)

B(q) = 1 - 2 q^-1                                                                                   

This model was not estimated from data.           

Sampling interval: 1 

Y dónde está cutil64.dll??

Muchos de los ejemplos que circulan por internet utilizan esta DLL. Tiene utilitarios generales para CUDA y viene dentro del "GPU Computing SDK examples". Para facilitar las cosas, es conveniente tomar esa librería (la versión correcta y para la plataforma correcta) de su carpeta en el SDK y colocarla dentro de la carpeta de instalación de CUDA. Lo mismo ocurre con los correspondientes archivos .lib.
Entonces, para un sistema 64 bits, copiamos:

• <NVIDIA GPU Computing SDK Folder>\C\inc\cutil.h en <CUDA Folder>\include
• <NVIDIA GPU Computing SDK Folder>\C\bin\win64\Release\cutil64.dll en <CUDA Folder>\bin64
• <NVIDIA GPU Computing SDK Folder>\C\bin\win64\Debug\cutil64D.dll en <CUDA Folder>\bin64
• <NVIDIA GPU Computing SDK Folder>\C\common\lib\cutil64.lib en <CUDA Folder>\lib64
• <NVIDIA GPU Computing SDK Folder>\C\common\lib\cutil64D.lib en <CUDA Folder>\lib64

A modo anecdótico, dejo unos comentarios sobre la cutil.dll del release notes de una versión de CUDA:
      - Most samples link to a utility library called "cutil" whose source code
        is in "NVIDIA CUDA SDK\common". The release and emurelease versions of
        these samples link to cutil[32|64].lib and dynamically load
    cutil[32|64].dll. The debug and emudebug versions of these samples link
    to cutil[32D|64D].lib and dynamically load cutil[32D|64D].dll.
        To build the 32-bit and/or 64-bit, release and/or debug configurations
    of the cutil library, use the solution files located in
        "NVIDIA CUDA SDK\common". The output of the compilation goes to
        "NVIDIA CUDA SDK\common\lib":
         - cutil[32|64].lib and cutil[32D|64D].lib are the release and debug
       import libraries,
         - cutil[32|64].dll and cutil[32D|64D].dll are the release and debug
           dynamic-link libraries, which get also copied to
           "NVIDIA CUDA SDK\bin\win[32|64]\[release|emurelease]" and
           "NVIDIA CUDA SDK\bin\win[32|64]\[debug|emudebug]"
       respectively;

Volviendo el tiempo atrás... CUDA 2.3

Lamentablemente, y aunque la sintaxis de CUDA 4 está muy buena, debemos volver a la versión 2.3 sobre Windows y compilar el código del Onera [1] con una resolución del Lucas-Kanade en GPU.

Lamentablemente, como suele suceder, aparecen nuevos problemas a salvar:

• NVCC sin un compilador no funciona. Necesitamos el Visual Studio...
• Visual Studio 2008 Express no tiene soporte para procesadores x64 con lo que hay que optar por una versión Trial del Professional.
• Para compilar por línea de comandos nos quedamos con un ejemplo pequeño de archivo .cu y logramos compilarlo sólo luego de leer un post interesante [2] en el foro de nvidia: es necesario ejecutar un batch que setea el entorno del VS en modo amd64 y luego indicar algunos parámetros extras al NVCC.
• Dentro del VS2008, NVCC sigue sin compilar. Es necesario actualizar el archivo nvcc.profile para que tenga un include extra del visual:

INCLUDES += "-I$(TOP)/include" "-I$(TOP)/include/cudart" "-IC:/Program Files (x86)/Microsoft Visual Studio 9.0/VC/include" $(_SPACE_)

[1] Onera - FOLKI GPU
[2] Post compilación x64

Eligiendo la cámara - Conclusiones

Ganador momentaneo: PixelFly

Conclusiones generales:

1. La Pulnix es demasiado lenta. Además no ofrece ninguna ventaja significativa.
2. La PixelFly tiene un buen tiempo de transmisión con binning activado. Según las pruebas de concepto se obtendrían 10 ms entre frame y frame con lo cual se puede usar esta cámara.
3. La PixelFly tiene además modo de doble shutter que no fue analizado pero brinda más alternativas.
4. La SpeedCam, aún siendo la de mejor calidad y velocidad no tuvo buenos resultados. La transferencia resulta muy lenta con lo que queda descartada por el momento.
5. Es importante aclarar que el análisis se basa en la necesidad de hacer un algoritmo lo más veloz posible y poder realiar cálculos en el momento. En caso de no ser esto así se puede reconsiderar los resultados y optar por la SpeedCam en lugar de la PixelFly, por ejemplo. Esto podría ocurrir en caso que se necesite tiempo extra para realizar procesamiento, no requerir 1 frame cada 10 ms sino 2 frames espaciados en 10 ms para, más tarde, repetir la toma de fotos y el proceso.

Tiempos - Usando la SpeedCam programáticamente

En este punto, la experiencia se convirtió en algo bastante pesado. La SpeedCam no tiene una API de acceso programática. Tampoco permite sacar fotogramas de forma unitaria, sólo se puede sacar una ráfaga que es almacenada en el buffer interno de la cámara y luego transmitido a la PC.
La cámara posee código hecho en Java del cual se pueden obtener los .jar para emplearlos en un código de prueba.
Consultamos entonces a los fabricantes: nos responden que la compañía cambió de manos. Consultamos a los nuevos responsables y obtenemos como información que tenemos tres caminos:

• utilizar una API a nivel registro con la cual se accede a la cámara enviando códigos de operación
• revisar los archivos minivis.jar y minivis.dll
• comprar otra cámara o adquirir el nuevo software.

Obviamente nos quedamos con la segunda por "simplicidad" (si es que se puede decir que esto sea simple). Entonces, armamos un archivo .Java y tratamos de adivinar cómo se usa la API que contiene el .jar.
Luego de varios días de revisar los resultados del Java Decompiler. Obtenemos una sucesión de pasos que permite obtener las imágenes con resultados bastante malos:

• Snapshots tomados: 200
• fps entre Snapshots: 2500
• Tiempo Total: 7800ms (sólo en la transferencia, ya que la toma de imágenes se debe hacer previa a la este paso)
• Tiempo Promedio: 39ms (es decir, 25 fps)

Conclusión: muy mala velocidad para una cámara tan potente. Claramente, no está optimizada para utilizarla de forma programática en entornos real-time. Por el contrario, el buffer parece consumir mucho tiempo de la transferencia.
Pero entonces surge una duda ¿Y si los tiempos son lentos por estar utilizando java? Para tratar de falsear esta teoría, recurrimos a código C++ que intente levantar las DLLs directamente y no utilizar minivis.jar. Nuevamente tenemos que decompilar, en este caso decompilamos el archivo minivis.dll y hacemos un matcheo contra lo observado en minivis.jar. Igualmente, debemos notar que ambos archivos se comunican gracias a que se empleó JNI en C++. Lamentablemente esto simplemente complica la tarea de decompilación y reuso de la DLL original.
Con una sola prueba de concepto nos queda claro que ese no es el mejor camino:

Run-Time Check Failure #0 - The value of ESP was not properly saved across a function call.  This is usually a result of calling a function declared with one calling convention with a function pointer declared with a different calling convention

Aquí el código Java:

   1:  class DummyConsumer implements CallbackConsumer {

   2:      public void modeChanged(int paramInt){

   3:      }

   4:      public void statusChanged(int paramInt) {

   5:      }

   6:      public void shutterChanged(boolean paramBoolean, short paramShort) {

   7:      }

   8:  }

   9:   

  10:  public class Application {

  11:      public static void main(String[] arguments) {

  12:          int numberOfImages = Integer.parseInt(arguments[0]);

  13:          

  14:          MinivisFactory factory = MinivisFactory.getInstance();

  15:          try {

  16:              factory.discover();

  17:              TypedNetAddress[] addresses = factory.getKnownDevices();

  18:              String mac = addresses.length > 0 ? addresses[0].getEntryName() : "00-50-C2-1D-7E-AB";

  19:              DummyConsumer callbackConsumer = new DummyConsumer();

  20:              Proxy camera = factory.getCamera(mac, callbackConsumer);

  21:              try {

  22:                  camera.connect();

  23:                  try {

  24:                      takeSnapshots(camera, numberOfImages, debug);

  25:                  }            

  26:                  finally {

  27:                      camera.disconnect();

  28:                  }

  29:              }

  30:              finally {

  31:                  camera.release();

  32:              }

  33:          }

  34:          catch(Exception e) {

  35:              e.printStackTrace();

  36:          }

  37:      }

  38:      private static void takeSnapshots(Proxy camera, int numberOfImages, boolean debug) throws IOException, MDriverException, Exception {

  39:          camera.setMode(4); //0 LIVE 1 RECORD 2 TRIGGERED 3 DRAM 4 LOWLIGHT

  40:          camera.trigger();

  41:          

  42:          //missing code here: wait for trigger to finish

  43:          

  44:          MFrameInfo frameInfo = new MFrameInfo();

  45:          for(int i = 0; i < numberOfImages; ++i) {

  46:              byte[] frameBytes = camera.getFrame(-1, -1  , frameInfo);

  47:          }

  48:      }

  49:  }

Aquí el código C++:

   1:   

   2:  #include "stdafx.h"

   3:  #include <iostream>

   4:   

   5:   

   6:  int _tmain(int argc, _TCHAR* argv[])

   7:  {

   8:      HMODULE hMod = LoadLibrary("minivis.dll");

   9:      typedef long (*MinivisFactory_Create)();

  10:      typedef long (*GetCamera)(long, char*, void*);

  11:      MinivisFactory_Create n_MinivisFactory_Create = (MinivisFactory_Create)GetProcAdress("_Java_com_artho_visart_plugins_minivis_internal_jni_MinivisFactory_n_1MinivisFactory_1Create@8");

  12:      GetCamera n_GetCamera = (GetCamera)GetProcAdress("_Java_com_artho_visart_plugins_minivis_internal_jni_MinivisFactory_n_1GetCamera@24");

  13:      

  14:      

  15:      long adapterPtr = n_MinivisFactory_Create();

  16:      void* consumer = NULL;

  17:      long pointer = n_GetCamera(adapterPtr, "00-50-C2-1D-7E-AB", consumer);

  18:      

  19:      std::cout << "Camera pointer:" << pointer << std::endl;

  20:      std::cout << "Press any key to continue..." << std::endl;

  21:      char character;

  22:      std::cin >> character;

  23:      return 0;

  24:  }