Multiplicación en punto fijo Pt. 2 [formato fixdt(signed,n,s,b)]

Esta entrada es una continuación del tema de multiplicación punto fijo y tocará el tema de un formato de mapeo lineal que permite ajustar números enteros a un rango lineal arbitrario como puede ser los rangos de voltaje estándar para DAC's y ADC's. El formato es el siguiente:

fixdt(signed,n,s,b) 

• signed: toma el valor de 1 si el entero almacenado es signado y 0 en caso contrario.
• n:  longitud en bits del entero almacenado. 
• s: pendiente de la función de mapeo lineal. 
• b: bias de la función de mapeo lineal.

En esta convención un número con parte entera y decimal se representa de la siguiente manera:

 El rango de este formato queda definido por sus parámetros de la siguiente manera:

 La multiplicación queda definida cómo:

 La multiplicación es mas complicada que en el formato fixed(signed,n,l) pero se simplifica cuando b = 0 siendo en este caso similar a la multiplicación en el otro formato. El siguiente paso en definir el procedimiento para computar los valores de n, s y b según el intervalo continuo en que el que deseemos trabajar. Al igual que el formato anterior ,el parámetro n define los valores máximos para los enteres almacenados. Por tanto, teniendo definidos los límites del intervalo real y del intervalo entero, los valores de s y b se calculan resolviendo el siguiente sistema de ecuaciones:
 

 Realizaremos un ejemplo práctico para demostrar la utilizad de este formato su ventaja con respecto a fixed(signed,n,l). Implementaremos en un FPGA un generador de señal cosenoidal (he escrito una entrada sobre los detalles de estos generadores). Utilizaremos una DAC LTC1661 que tiene una resolución de 10 bits y un rango de salida de 0-Vreff*(1023/1024). Las salidas Vcc de la tarjeta que voy a utilizar son de 3.3 volts por lo que este intervalo será aproximadamente de 0-3.3 V. Ya que esta DAC no puede dar salidas de voltaje negativo el valor para signed sera '0'. Con esta información podemos encontrar los valores para establecer nuestro formato. Podemos obtenerlos mediante este código de Matlab:

minimo = 0; %Volts
maximo = 3.3; % Volts
is_signed = false;
word_length = 10;

[E_min, E_max] = range(fi([], is_signed, word_length, 0));

A = double ([E_min, 1; E_max, 1]);
b = double ([minimo; maximo]);
x = A\b;
format long g
slope = x(1)
bias = x(2)

La ejecución del código nos da una pendiente s = 0.0032258064516129 y un bias b = 0. La señal para este ejemplo tendrá una frecuencia de 0.05 Hz y un periodo de muestreo de 1 segundo (si pueden usar un osciloscopio pueden probar con una frecuencia más alta y un periodo de muestro más corto). Generaremos un código VHDL con HDL Coder a partir de este modelo de Simulink:

Los subsitemas que he llamado Fixed-Point Multiplication contienene la implementación de la multiplicación que hemos descrito aquí. Dado que b = 0 la multiplicación de punto fijo simplemente se define como el producto de los enteros almacenados y la pendiente s:

 Algo que podría resultar confuso y difícil de entender al principio es el por qué el producto con s se realiza en formato fixed(signed,n,l). La respuesta es que ese formato se utiliza para principalmente para la multiplicación de números fraccionales en un entorno que solo permita trabajar con números enteros de forma interna (que es el caso de nuestro sistema de hardware). El formato de mapeo lineal fixdt(signed,n,s,b) sólo se encarga de que las operaciones entre enteros tengan el sentido físico que esperamos de forma externa. Esto tambien explica un segundo hecho que tambien puede resultar confuso: las operaciones entre enteros se realizan en modo signed (para facilitar las restas) aunque se ha definido el modo unsigned. 

Para generar el código nos vamos a Code > HDL Code > Generate HDL. Para la comunicación SPI con la LTC1661 he modificado un código que había publicado anteriormente. Los módulos son los siguiente:

• Divisor
• DAC_SPI
• Generador_Coseno_HDL (Generado por HDL Coder) 

El esquemático del sistema es el siguiente:

Nota: durante la generación del código se indica que el reloj para el modulo del generador sea el doble del establecido. En este caso el divisor tiene una salida de 2 Hz.

Según la simulación, y si implementamos correctamente las cosas, las secuencia de voltaje que deberíamos ver experimentalmente a la salida del DAC es la siguiente:

Aquí muestro un video donde se puede ver que no tengo un buen multímetro y tuve que aumentar a 4 segundos el periodo de muestreo para dar tiempo de estabilizar las lecturas:
Referencias
Compute Slope and Bias, MathWorks
Fixed-Point: Rage and Precision, MathWorks 

Multiplicación en punto fijo [formato fixdt(signed,n,l)]

Aunque al principio pueda chocar con la intuición, es posible realizar de forma digital una multiplicación de dos números con parte entera y decimal utilizando únicamente números enteros. De hecho, para aplicaciones de instrumentación científica o simulaciones numéricas dónde la precisión sea la prioridad es incluso aconsejable evitar el uso del famoso punto flotante. Recomiendo enormemente ver este capítulo del canal PBS Infinite Series para entender desde el punto de vista matemático el problema fundamental de los tipos double y single. En pocas palabras este problema radica en la en las limitaciones en la accesibilidad sobre la recta numérica computable. La resolución del punto variable no es homogénea a lo largo de todo su rango. Una segunda limitación de las operaciones de punto flotante tiene que ver a su costo de implementación en hardware; construir circuitos digitales de punto flotante consumen más recursos en un FPGA.

El formato de punto fijo con el trabajaremos en esta entrada es:

 fixdt(signed,n,l)

• signed : toma el valor de 1 si el entero almacenado es signado y 0 en caso contrario
• n : longitud en bits del entero almacenado 
• l : longitud en bits de la parte fraccional

En esta convención un número con parte entera y decimal se representa de la siguiente manera:

Los parámetros n y signed nos da el rango del formato:

A partir de la relación entre el Valor del Mundo Real (como se le llama en la jerga de punto fijo y denotaremos como Val) y el entero almacenado (que denotaremos por E) podemos definir la multiplicación de punto fijo de la siguiente forma:

Esto significa que la multiplicación se realiza como el producto de dos enteros mas un bitshift a la derecha de l locaciones. Para comprender mejor este formato realizaremos un ejemplo práctico concreto. Supongamos que queremos digitalizar una señal de voltaje sinusoidal de amplitud de 2 V con una frecuencia de 1 Hz mediante un ADC de 12 bits de resolución y un rango de entrada de 0-4 V. Supongamos que una vez digitalizada la señal queremos realizar un escalamiento de 0.321 y enviar esta nueva señal a un DAC de la misma resolución y rango. ¿Cómo resolveríamos el problema sin utilizar operaciones de punto flotante?

Bien, usaremos el formato fixdt(signed,n,l) para este ejemplo. Lo primero que debemos hacer es establecer los parámetros del formato. Ya que nuestro ADC no puede leer voltajes negativos, el parámetro signed será igual a 0. Dado que la resolución es de 12 bits entonces n = 12. Sabiendo que el voltaje máximo de lectura es 4 V encontramos el valor l de la siguiente manera:

Tenemos que el valor de l es 10 dando un rango de 0-3.99 V. Entonces el formato para nuestro ejemplo queda definido así: fixdt(0,12,10). El siguiente paso en convertir la ganancia de 0.321 en formato fixdt(0,12,10):

El equivalente de la ganancia resulta ser: 329. Podemos realizar ahora una simulación de nuestro ejemplo en Simulink y comparar el desempeño de la operación de punto fijo contra la de punto flotante:

Click para agrandar

El opam en el circuito añade el offset necesario para levantar la señal de la fuente y pueda ser adquirida por el ADC. Los bloques de ADC y DAC son simplemente subsistemas que realizan las operaciones de conversión de decimal a fixdt(1,12,10) y de muestro ZOH. Si bien Matlab/Simulink permite realizar las conversiones al formato fixdt(signed,n,l), decidí manejar uint16 para poder generalizar el ejemplo para ser implementado en microcontroladores y no sólo en FPGA's. Nuestro formato fixdt(0,12,10) cabe en un uint16 y eso también permite hacer una intuición más clara de que simplemente trabajamos con formatos enteros convencionales. La gráfica resultante con la comparativa es la siguiente:

Vemos que la correspondencia entre las operaciones de punto fijo y punto flotante son muy buenas. Podemos ahora a aventurarnos a escribir un codigo VHDL para nuestra operación de escalamiento de este ejemplo:

La simulación del módulo mostrando algunos valores de ejemplo para los enteros almacenados de entrada y salida es la siguiente:

 Click para agrandar

Para finalizar quisiera hacer un último comentario. Quizá notaron que tuvimos mucha suerte al calcular el valor del exponente l y nos dio un número que se ajustaba perfectamente al rango del ADC. No fue suerte, yo elegí el rango para que este ejemplo fuera más sencillo. En la practica los DAC's suelen tener rangos de 0-5.2 V o 0-3.3 V. ¿Que se hace en esos casos? Aquí es donde entra un nuevo formato: fixdt(signed,n,s,b). Pero hablaremos de él en otra entrada.

Referencias:

Fixed-Point: Rage and Precision, MathWorks 

Convertir un modelo de Simulink a código VHDL con HDL Coder [Ejemplo 2]

En una entrada anterior exponía el como generar un código VHDL de una función booleana sencilla. Este será también un ejemplo sencillo pero mucho más interesante. Generaremos un controlador P para la posición de un motor DC con caja de engranes y encoder [este modelo en particular] alimentado por un driver Pololu Qik 2s12v10. El esquema del sistema es el siguiente:

 La posición del motor será leída mediante el encoder que viene integrado en el motor. El modelo de Simulink que convertiremos a código VHDL es el siguiente:

Bloques requeridos:
HDL Coder > Commonly Used Blocks > In1
HDL Coder > Commonly Used Blocks > Out1
HDL Coder > Sources > Constant
HDL Coder >  Math Operations > Divide
HDL Coder > Discontinuities > Saturation
HDL Coder > Signal Attributes > Data Type Conversion

Diagrama:

 La operación de división por números que no sean potencias de 2 no está disponible en todos los dispositivos por lo que en este caso usaremos 1/32 como aproximación de 0.03 [es posible implementar en hardware divisiones entre números arbitrarios pero serán tratados en otro ejemplo]. Es necesario establecer en todos los bloques el tipo de dato int16 dando doble click a cada bloque y cambiando el tipo de dato la siguiente forma:

 También sera necesario dar doble click en el bloque Divide, ir a Signal Attributes y activar el modo "Saturate on integer overflow" y verificar que el modo de redondeo sea "zero" o "simplest". Habiendo hecho todo lo anterior generemos el código VHDL yéndonos a la pestaña Code > HDL Code > Generate HDL. Les generará el siguiente código:

Para complementar el controlador usaremos los siguientes módulos VHDL que ya he publicado por acá:

• Módulo UART
• Módulo Qik2s12v10
• Módulo Encoder X4 

Agregando todos los archivos .vhd requeridos al proyecto se puede proceder a conectar los bloques de forma esquemática en ISE de la siguiente manera:

Si no están familiarizados con la creación de proyectos esquemáticos en ISE aquí hay un video-tutorial en español muy bien explicado [utilizan Verilog pero el procedimiento es el mismo]. El bloque selector es simplemente una declaración when que envía una referencia de 0 si el switch esta apagado y de 2400 (90° en cuentas de encoder) cuando esta encendido. En este video muestro el funcionamiento del sistema implementado en una tarjeta Basys2:

Convertir un modelo de Simulink a código VHDL con HDL Coder [Ejemplo 1]

Para este primer ejemplo, generarémos el código VHDL para una función booleana sencilla de la forma X = AB xor C desde un diagrama en Simulink.

Bloques requeridos
HDL Coder > Commonly Used Blocks > In1
HDL Coder > Commonly Used Blocks > Out1
HDL Coder > Logic and Bit Operations > Logical Operators

Diagrama

Una vez terminado el bloque debemos irnos a la pestaña Code > HDL Code > Generate HDL. Si no hay ningún error deberá aparecer lo siguiente en la ventana de comandos:

### Generating HDL for 'Ejemplo_HDL'.
### Starting HDL check.
### Begin VHDL Code Generation for 'Ejemplo_HDL'.
### Working on Ejemplo_HDL as hdlsrc\Ejemplo_HDL\Ejemplo_HDL.vhd.
### Creating HDL Code Generation Check Report Ejemplo_HDL_report.html
### HDL check for 'Ejemplo_HDL' complete with 0 errors, 3 warnings, and 0 messages.
### HDL code generation complete. 

Por defecto el código generado aparecerá en el directorio  hdlsrc\. El código generado para este ejemplo es el siguiente:

Identificación de la función de transferencia de un motor de Lego Mindstorms NXT

Utilizando la App de System Identification Toolbox de Matlab es muy sencillo estimar experimentalmente la función de transferencia de un sistema dinámico. Básicamente, existen 3 formas de identificar la planta de un sistema dinámico: el modelo de caja negra, el modelo de caja blanca y el modelo de caja gris. En el modelo de caja blanca se obtiene la función de transferencia a partir de los principios físicos conocidos sobre el sistema (que es lo que suele hacerse en un primer curso de sistemas de control al obtener la planta de un sistema masa-resorte-amortiguador). Sin embargo, cuando se trata de sistemas reales esto o suele ser muy complicado. El método de modelado por caja negra se encarga de estimar la función de transferencia con herramientas estadísticas utilizando únicamente una función de entrada de prueba y la función de respuesta de la planta ante dicha función:

Para la que la identificación sea certera es un requisito que u(t) tengo un ancho de banda muy amplio. Esto le permite al algoritmo de estimación tener más información disponible para la estimación del modelo. Idealmente la función que cumple con ancho de banda infinito es la delta de Dirac, sin embargo, al ser imposible de generar físicamente se opta en la practica por usar secuencias pseudoaleatorias binarias (PRBS's). Estas secuencias pueden generarse con un registro de corrimiento de retroalimentación lineal (LFSR) y tienen características espectrales de ruido blanco lo que las hace ideales para ser usadas como señales de prueba. La descripción de los generadores de PRBS's es todo un tema aparte del que espero escribir a detalle más adelante, pero pueden darse muy buena idea de como funcinan con este documento.

Para esta práctica se requiere tener instalar el paquete de soporte de  Lego Mindstorms para Simulink. Pueden instalarlo ejecutando el comando supportPackageInstaller en la ventana de comandos de Matlab y buscando Lego Mindstorms NXT. El modelo implementado en el NXT es el siguiente:

 








Este diagrama nos permite adquirir los datos tanto de la función de entrada como la de salida y los envía al workspace. Para poder hacer la adquisición a la PC debe ejecutarse en modo external. [Aquí pueden ver un video tutorial para hacer esto]. Los datos obtenidos son los siguentes:

Estas dos señales son las que utilizaremos en al app para realizar la identificación. Esta aplicación tiene una interfaz gráfica muy sencilla de utilizar.[Aquí hay un videotutorial en español para aprender a usarla]. Después de experimentar con varias combinaciones de polos y ceros, el modelo preliminar que logró un acierto del 96.12 % fue el siguiente:

Este aún es el modelo crudo encontrado por el identificador y es necesario ajustarlo. Lo primero que salta a la vista el el termino 5,434e-11 que es despreciable y eliminaremos en el modelo final. Lo que resta es hacer un ajuste de ganancia. Para ello requerimos la respuesta experimental en lazo cerrado del motor ante un escalón unitario. Ejecutamos entonces el siguiente diagrama a bloques para poder obtenerlo:

Teniendo los datos de la respuesta experimental, debemos ahora compararlo con la respuesta teórica de nuestro modelo ajustando la ganancia hasta observar resultados favorables (debe haber algún procedimiento analítico para esto, pero lo he realizado de forma empírica). El diagrama del modelo teórico, junto con el modelo de caja blanca realizado por la Universidad  Politécnica de Madrid[1] como segunda referencia, es el siguiente:

 Las respuestas fueron las siguientes:

   Con esto queda verificado que la función de transferencia final de posición angular (en grados) contra potencia (%) nos queda:

 La ventaja de este proceso es que puede repetirse para caracterizar la función de transferencia de los motores con carga como puede ser un brazo robótico o un robot móvil, algo que representaría un proceso tedioso con el enfoque de caja blanca. Adelantantome a la petición de los archivos:

Descargar Modelos

Implementación de un controlador PID discreto en un PIC18F2550

En esta entrada se utilizará el método desarrollado en la entrada anterior, por lo que saltarán algunos pasos durante el proceso. En esa entrada se implementó un controlador PID para una planta descrita por un doble integrador. En esta ocasión haremos una simulación de un controlador PID implementado en un PIC para controlar una planta construida con amplificadores operacionales. Si bien es posible construir un doble integrador físicamente, resulta impráctico debido a la saturación provocada por la acumulación continua. Para resolver esto, esta planta estará caracterizada por dos integradores con reset resistivo o integradores de ganancia controlada (que una forma elegante de llamar a un filtro pasa-bajas activo de un solo polo) [1].  Eligiendo los valores de resistencia y capacitancia de tal manera que cada integrador tenga ganancia unitaria, la función de transferencia de los 2 integradores nos queda:

 Para controlar esta planta utilizaremos el siguiente controlador discreto a partir de la transformación backward Euler:

Dónde:

Obteniendo transformada inversa, la salida de control descrita como ecuación de diferencias nos queda:

Hacemos la siguiente sustitución de coeficientes:

Estos coeficientes pueden encontrarse manualmente siguiendo el procedimiento descrito aquí, o ejecutando éste programa de Matlab: [PID_Proteus.m]. Teniendo estos valores podemos realizar una verificación de la ecuación en Simulink:

Algo muy importante que se debe notar es que en esta simulación se está considerando la limitación física de la saturación de la salida del controlador real. Sin esta limitación la respuesta del controlador sería más rápida y con un sobreimpulso mucho menor. Sabiendo ahora que nuestra ecuación de diferencias es correcta podemos proceder a escribir el programa en el PIC en lenguaje C con el compilador MikroC. El programa completo se adjuntará más adelante. La función de control principal queda descrita como:

void pid(void) // Función para el controlador
{
   e2=e1; e1=e0; u2=u1; u1=u0; // update de variables
   y = ADC_Read(0);   // Lectura analógica por AN0
   yd = (double) y;  // Conversión a double
   vin = (yd/1023*5);
   y = ADC_Read(1);   // Lectura analógica por AN1
   yd = (double) y;  // Conversión a double
   reff = (yd/1023*5);
   e0 = reff - vin;  // Error
   u0 = ke0*e0 + ke1*e1 + ke2*e2 - ku1*u1 - ku2*u2; // Ecuación de diferencias
   salida = (int) floor(((u0+5)/10)*1023);
   if (salida > 1023) salida = 1023;  // Saturación de DAC
   if (salida <= 0) salida = 0;
   DAC_Output(salida); // Enviar dato
}

Es muy importante que ésta función sea ejecutada dentro de una rutina de interrupción de manera que se tenga control del tiempo de muestreo Ts en exactamente 0.01 segundos para que el controlador funcione correctamente. Para la salida analógica del controlador se esta utilizando un DAC LTC1660 [2][3]. Compilando el programa y simulando en Proteus:

 Nota: El relay en la simulación lo único que hace es evitar que la salida bipolar del DAC envíe un voltaje de saturación negativo al inicio [ver este programa].

Descargar archivos

Implementación de un controlador PID discreto

Supongamos que queremos controlar el mismo modelo de una masa móvil sin fricción que habíamos tratado en esta entrada cuya función de transferencia es la siguiente:

Dados las ganancias del controlador PID que se habían considerado (recomiendo utilizar el app PID Tuning de Matlab si no tienes estos valores para la planta de tu interés) :

La función de transferencia en tiempo continuo de nuestro controlador nos queda:

En la práctica es recomendable utilizar un LPF en la etapa derivativa a modo de reducir el ruido. En este caso la FT es:

Dónde N representa la frecuencia de corte del filtro que suele ser de 100 radianes/s. Teniendo la función de transferencia de nuestro controlador, el proceso de implementación de este controlador en un sistema digital (PIC, Arduino, FPGA, etc) es el siguiente:

La ecuación de diferencias es la que describirá el algoritmo que deberemos implementar en algún lenguaje de programación o un lenguaje de descripción de hardware. Para realizar la transformación entre la FT en dominio de Laplace al dominio de Z debemos considerar las siguientes aproximaciones numéricas para la derivada:

De aquí es fácil notar que podemos representar a la variable s en terminos de z (dónde Ts es el periodo de muestreo) de las siguientes maneras:

Estas 2 transformaciones de la variable s tienen propiedades distintas. Mientras la transformación backward Euler es estable siempre que la FT en tiempo continuo sea estable, la transformación forward Euler no lo es en todos los casos. (Para más detalles consultar la sección 19.2.3.1 de "Pasive, Active and Digital Filters", Wai-Kai Chen [1]). Para este ejemplo utilizaremos la transformación forward Euler, por lo que nuestro controlador queda representado en tiempo discreto por la TF:

Para este ejemplo se sabe de antemano que este controlador es estable pero si se arriesga a utilizar foward Euler se recomienda verificar estabilidad (puede usarse sisotool de Matlab).  Verificamos en Simulink que nuestro controlador funcione correctamente:

Ahora podemos proceder a la conversión de formatos de la FT. Para eso utilizamos el este programa de Matlab: PIDdiscreto.mat. Usando los datos que nos muestra el programa la FT nos queda:

Reordenando (click para agrandar):

Calculando transformada inversa:

Finalmente, verificamos que la ecuación de diferencias sea correcta en Simulink:

En una próxima entrada se implementará este controlador en un código en C para PIC y Arduino.

Simulación de circuitos eléctricos en Simulink

Bloques usados:

Solver Configuration: (Simscape>Utilities)
Voltage Sensor: (Simscape >Foundations Library>Electrical>Electrical Sensors)
Current Sensor: (Simscape >Foundations Library>Electrical>Electrical Sensors)
PS-Simulink Converter: (Simscape >Utilities)
Resistor  (Simscape>Foundations Library>Electrical>Electrical Elements)
Electrical Reference  (Simscape>Foundations Library>Electrical>Electrical Elements)
AC Voltage Source     (Simscape>Foundations Library>Electrical>Electrical Sources)
Mux  (Simulink >Common Used Blocks)
Scope  (Simulink>Sinks)

Para este modelo se dejó la configuración del Solver en default con la fuente de voltaje a 10 V pico-pico a 60 Hz y la resistencia es de 50 ohms. El tiempo de simulación es de 0.06 segundos.

Nota: Cómo activar las leyendas en el Scope.