Discretización: cómo programar un PID en C++

La discretización es el proceso de convertir ecuaciones continuas en pasos discretos para que puedan ejecutarse en un microcontrolador o un ordenador.

En los artículos teóricos hemos hablado con mucha alegría de “integrar el error” y “calcular la derivada”. El problema es que los microcontroladores que usamos habitualmente, como los de un Arduino o un ESP32, no tienen ni idea de cálculo diferencial. Ellos solo saben sumar, restar, multiplicar y dividir números en un bucle secuencial que se repite una y otra vez.

Hoy vamos a bajar de las nubes de la teoría de Laplace y vamos a mancharnos las manos de silicio. Vamos a ver cómo convertir esas matemáticas en un código C++ limpio y eficiente.

Del modelo continuo al modelo digital

Para refrescar la memoria, la ecuación clásica del controlador PID en el dominio del tiempo continuo es:

En esta ecuación, el tiempo t fluye de manera suave y constante. Pero en nuestro microcontrolador, el tiempo avanza “a saltos”. Leemos el sensor, calculamos, escribimos en el actuador, y volvemos a empezar. A cada uno de estos ciclos lo llamamos tiempo de muestreo o Δt (dt).

Nuestra misión es sustituir la función continua e(t) por una lectura instantánea en un momento exacto e[k], donde k es el ciclo actual del bucle, y k−1 es el ciclo anterior.

Acción proporcional (P)

Esta es la parte fácil. La acción proporcional no tiene memoria, solo depende del instante presente. Por tanto, su traducción a código es directa. Simplemente multiplicamos el error actual por la ganancia:

Acción integral (I): sumar rectángulos

Aquí la cosa se pone interesante. La integral de una señal es, por definición geométrica, el área bajo la curva que forma esa señal.

Como nuestro microcontrolador no puede ver la curva perfecta, solo ve un punto nuevo cada Δt milisegundos. Para calcular el área, usamos una técnica numérica llamada Integración de Euler. Consiste en asumir que durante el tiempo Δt, el error se ha mantenido constante, formando un pequeño rectángulo de base Δt y altura e[k].

El área total acumulada es simplemente el área que ya teníamos guardada más el área de este nuevo rectángulo:

Acción derivativa (D): calcular pendientes

La derivada representa la velocidad de cambio o la pendiente de la curva en un instante dado.

En nuestro mundo de pasos discretos, si queremos saber a qué velocidad está cambiando el error, simplemente restamos el error actual menos el error que tuvimos en el ciclo anterior, y lo dividimos por el tiempo que ha pasado entre ambas medidas (otra vez nuestro querido Δt):

Para que esto funcione en código, al final de cada ciclo, tenemos que acordarnos de hacer error_anterior = error_actual, guardándolo para el futuro.

Implementación en C++ para Arduino/ESP32

Ahora sí, vamos a plasmar esto en código real. Para calcular el tiempo transcurrido (Δt), en el ecosistema Arduino solemos usar la función millis() (o micros() para sistemas muy rápidos).

Una solución tentadora sería poner un delay(10). En un programa que comparte el procesador con comunicaciones, interfaz u otras tareas, conviene evitar delay() dentro del lazo de control, porque bloquea su ejecución y dificulta controlar el periodo real.

Vamos a crear una estructura básica usando el paradigma de «evaluación por marcas de tiempo»:

// Variables globales de sintonización del PID
float Kp = 2.0;
float Ki = 5.0;
float Kd = 1.0;

// Variables de estado interno del PID (Memoria)
float integral_acumulada = 0.0;
float error_anterior = 0.0;
unsigned long tiempo_anterior = 0;

// Variables del sistema
float setpoint = 100.0;

void setup() {
    // Inicialización de hardware
    tiempo_anterior = millis();
}

void loop() {
    unsigned long tiempo_actual = millis();
    unsigned long dt_ms = tiempo_actual - tiempo_anterior;

    // Solo ejecutamos el cálculo si ha pasado un tiempo mínimo
    // Por ejemplo, ejecutamos el PID cada 10ms (100 Hz)
    if (dt_ms >= 10) {
        
        // 1. Convertir dt a segundos (fundamental para que las constantes Ki y Kd tengan sentido físico)
        float dt = (float)dt_ms / 1000.0; 
        
        // 2. Medir la variable de proceso
        float pv = leerSensor();
        
        // 3. Calcular el error
        float error = setpoint - pv;
        
        // --- CÁLCULO DEL PID ---
        
        // Proporcional
        float P = Kp * error;
        
        // Integral (Euler integration)
        integral_acumulada += error * dt;
        float I = Ki * integral_acumulada;
        
        // Derivativa (Diferencia finita)
        float derivada = (error - error_anterior) / dt;
        float D = Kd * derivada;
        
        // Salida total
        float salida = P + I + D;
        
        // --- FIN DEL CÁLCULO ---
        
        // 4. Actuar sobre la planta
        escribirActuador(salida);
        
        // 5. Guardar estado para el próximo ciclo
        error_anterior = error;
        tiempo_anterior = tiempo_actual;
    }
    
    // Aquí podemos poner otras tareas no bloqueantes
    // leer_puerto_serie();
    // actualizar_pantalla();
}

Notas críticas sobre la implementación

1. La conversión a segundos: Fíjate que dividimos dt_ms entre 1000.0. Esto es importante. Si no lo haces y multiplicas directamente por milisegundos, tu término integral se volverá inmenso en fracciones de segundo y tu derivada será minúscula, obligándote a usar valores de Ki y Kd extrañísimos e inestables.
2. Determinismo: El enfoque con millis() en el loop() está bien para robótica básica. Sin embargo, si tu microcontrolador está haciendo muchas cosas (como manejar un servidor web con WiFi en un ESP32), el tiempo dt_ms podría fluctuar entre 10 ms y 15 ms. Si necesitas reducir el jitter, usa un temporizador de hardware o una tarea periódica. Mantén la ISR breve y deja el cálculo y el acceso a periféricos en una tarea con prioridad y tiempo de ejecución controlados.

El elefante en la habitación

El código que acabamos de escribir reproduce la estructura matemática básica, pero todavía no está listo para una planta real. Le faltan límites de salida, anti-windup, filtrado de la derivada y un tratamiento explícito del primer ciclo para evitar un pico derivativo.

El mundo real añade además varios problemas. ¿Qué pasa si calculamos que el motor necesita “2000 Voltios” pero nuestra batería solo da 12V? La integral empezará a acumular un error gigantesco que nunca podrá borrar. ¿Y qué pasa si nuestro sensor tiene ruido eléctrico? Que la resta e[k]−e[k−1] de la derivada generará unos picos de cálculo salvajes que harán que nuestro motor chille y vibre como una lavadora estropeada.

Estos son los problemas reales de la implementación digital, y son exactamente lo que vamos a resolver en los dos últimos artículos del curso. Empezando por el temible “Windup” de la integral.