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.

Temas de actualidad

Yo que soy un derechairo-liberal declarado tengo la sospecha que el atentado en París ha sido un evento de falsa bandera. De ser verdad podría ayudarnos a nosotros, ratas de ese laboratorio sociológico llamado Facebook, a decidirnos de que lado estar en esta batalla virtual. Pero ahora, ¿hasta que punto una bandera representa al Leviatán o representa al pueblo? Aún resolviendo eso queda otro problema, ¿a que pueblo se le debe mostrar empatía? A todos en el caso ideal, por supuesto. Pero esta solución genera otro problema; por más utopías que quieran sacarse de la manga, los seres humanos no podemos vivir sin fronteras (en todos los contextos posibles). Siempre he tenido como lema personal esta frase de Huygens, "El mundo es mi país, la ciencia mi religión", pero me he dado cuenta que en un contexto no político, esto no es verdad para mí. Existen muchos grupos sociales que respeto (siendo terriblemente sincero, otros no) más no puedo identificarme. Cargo en mi mente, como todos, fronteras invisibles definidas por reglas de pertenencia social. Y aquí una triste verdad. Aunque lo correcto sea mostrar empatía por la humanidad entera, la paradoja de la elección viene a arruinaros la paz. Elegir todo no es diferente a elegir nada.

Guanajuato

Hace unos días estuve una escuela-taller en el CIMAT, lo que fue el pretexto perfecto para conocer el instituto y la ciudad. Me gustó bastante y fue una lastima tener tan pocos días para recorrerla. Si bien el taller me dejó un poco frío me llevé buena idea de las áreas en las que trabaja el CIMAT y bastante información del posgrado. Durante la estancia me quedé en la casa del asesores de un amigo mío, que por cierto conocí en un congreso en Guadalajara el año pasado (aquí aquel breve episodio). Su asesor es una persona realmente agradable, un astrofísico brasileño de la UG. El primer día me recibió dando una clase de tango a la que no pude rechazar la invitación de unirme pues era algo que había querido hacer desde hace mucho. Conocí también al housemate de mi amigo, uno de sus mejores amigos que estudió junto con el en la BUAP que igual me cayo muy bien. En las tardes, después de que regresaba de las sesiones taller, bajábamos caminado al centro. Algunas partes de la ciudad me recordaron bastante a Venecia. El penúltimo día dimos un recorrido por los cerros que rodean a la ciudad, algo que suelen hacer. En algún momento tuvieron que irme a rescatar cuando me quedé aferrado como gato asustado a una pared de piedra sin poderme mover. Lo más gracioso fue que en aquel momento sentí que estaba en risco de 300 metros cuando en realidad no estaba tan alto. Ya en lo alto la vista era increíble. Debo decir que mi amigo lleva una vida de estudiante envidiable. Para mi fue un corta pero buena aventura.

Existencialismo

Si me preguntaran cuales han sido los 3 libros que marcaron mi vida, simplemente diría que sólo hay un ensayo que realmente lo ha hecho, y lo ha hecho más tarde de lo que hubiera querido. Me refiero a "El Existencialismo es un Humanismo" de Sartre. Pasé casi toda mi vida culpando al mundo de mis fallas, pasando por alto una idea tan fundamental; la existencia del hombre precede a su esencia y somos nosotros mismos quienes nos definimos. Quizá en fondo lo sabemos y nos aterra la idea de haber nacido en un mundo indiferente que existe y debe definirse desde cero. Que tarde lo he conocido señor Sartre.

Convertir tiempo UTC a GMST (Tiempo Sideral Medio de Greenwich) en Python

 -*- coding: utf-8 -*-
"""
Created on Tue Oct 27 18:37:49 2015

@author: Rodolfo Escobar et al
"""
from datetime import datetime
from astropy.time import Time
from numpy import floor

# Reductor de rango 0-24h
def to24(x):
 while x >= 24:
  if x < 0:   
   x = x + 24
  if x > 24:
   x = x - 24
 return x
#

# Conversor de formato de hora
def tohms(x):
 time_hours = x
 time_minutes = time_hours * 60
 time_seconds = time_minutes * 60

 horas    = int(floor(time_hours))
 minutos  = int(floor(time_minutes % 60))
 segundos = int(floor(time_seconds % 60))
 return (horas,minutos,segundos)
#

#Programa principal
ut = Time(datetime.utcnow(), scale='utc')
JD= ut.jd # Fecha Juliana
D = JD - 2451545.0
GMST = 18.697374558 + 24.06570982441908*D
GMST = to24(GMST)
(horas,minutos,segundos) = tohms(GMST)
#

print "Hora Sideral:","{h}:{m}:{s}".format(h=horas, m=minutos, s=segundos)

Referencias:
Approximate Sidereal Time, USNO
Local Sidereal Time, Durham University
Astronomy With Your Personal Computer, Peter Duffett-Smith
Astropy.org

winvid en Matlab 2014/2015

Por alguna razón no viene instalada por defecto ninguna interfaz de video en las ultimas versiones de Matlab. Para resolver esto solo hay que escribir en la consola supportPackageInstaller e instalar el siguiente paquete:

Con esto podrás adquirir la imagen de tu webcam sin problemas después de reiniciar Matlab.

Interrupciones externas en Raspberry Pi con Python

A partir de la versión 0.5.6 de la librería RPi.GPIO es posible utilizar rutinas de interrupción externas. Como los usuarios de microcontroladores experimentados deben saber, una rutina de interrupción permite ejecutar una o una serie de instrucciones sin perder tiempo en bucles de lectura. Es decir que si el procesador se encuentra ejecutando una tarea, éste pausará el proceso en curso y ejecutará la rutina de interrupción ante un flanco de subida o bajada ocurrido en un pin GPIO. Si no has comprado tu Raspberry recientemente y no has actualizado su software es recomendable hacerlo:

sudo apt-get update 

sudo apt-get upgrade #(Esto podría tardar hasta más de una hora)

El siguiente código es un ejemplo sencillo sobre el uso de interrupciones. Tenemos dos rutinas de interrupción CuentaA() y CuentaB() asociadas a GPIO23 y GPIO24 respectivamente. Cada una incrementan una variable de conteo y muestran el resultado en consola. El programa se detiene cuando el contador A es mayor o igual a 5. Para más detalles sobre interrupciones pueden revisar la documentación de la librería RPi.GPIO. Para evitar el el efecto de rebote se utilizan capacitores de entre 22 o 10 uF conectados en paralelo con las resistencias pulldown de 10K como se muestra en el circuito [aunque también puede utilizarse un debounce vía software agregando el argumento bouncetime=200 a la función add_event_detect()]:

import RPi.GPIO as GPIO
import os

contaA = 0
contaB = 0

#Setup
GPIO.setmode(GPIO.BCM)
GPIO.setup(23,GPIO.IN)
GPIO.setup(24,GPIO.IN)

#Callbacks
def CuentaA(channel):
    global contaA
    contaA += 1
    os.system("clear")
    print "Contador A: ", contaA
    print "Contador B: ", contaB

def CuentaB(channel):
    global contaB
    contaB += 1
    os.system("clear")
    print "Contador A: ", contaA
    print "Contador B: ", contaB

#Interrupciones
GPIO.add_event_detect(23, GPIO.RISING, callback = CuentaA)
GPIO.add_event_detect(24, GPIO.RISING, callback = CuentaB)

print "Contador A: ", contaA
print "Contador B: ", contaB

#Bucle principal
while(contaA < 5):
    pass

GPIO.cleanup()

Simulación de una señal ECG y filtro de Notch en Matlab

Descargar programa comentado

Filtro de Notch de 60 Hz

Debido a la restricción de los valores comerciales para los resistores y capacitores, este circuito tiene una frecuencia central de 58.94 Hz. Sin embargo, como se puede ver en una simulación, se tiene una atenuación de -34 dB (un factor de de 0.0004) para la componente de 60 Hz.

Para más detalles sobre este tipo de filtros pueden consultar éste artículo.

Circuito para electrocardiografía (ECG/EKG)

Este circuito está basado en el que se puede encontrar en la hoja de datos del INA114. En el circuito original se asume que el rango de voltaje de la señal de entrada es de 1 mV ~ 5 mV más un offset de DC de 300 mV. Sin embargo, en la práctica estos valores resultan ser aún más pequeños (posiblemente por las condiciones poco controladas en un laboratorio de ingeniería universitario o el uso de caimanes y protoboard) así que es necesario incrementar la ganancia del amplificador de instrumentación. Para esta modificación se tiene una ganancia de 2501 y estos fueron los resultados:

En este video utilizan un circuito con ganancia de 100 con buenos resultados (el cual no funcionó en mi caso), lo que me hace sospechar que existen varios factores que afectan el voltaje de entrada. Posiblemente este circuito no funcionará en condiciones diferentes por lo que recomiendo experimentar con diferentes valores de ganancia.

Un segundo circuito que también funciona muy bien es una modificación del que pueden encontrar en The Biosignal How-To: