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