Ajuste del PID y filtrado de señales

01 · Punto de partida

De las ganancias a ojo a un método

En el tema 43 llegamos a unas ganancias para el motor probando en un orden razonable. Funciona, pero tiene dos problemas: depende de la experiencia de quien prueba, y no dice cuándo parar. Los métodos de sintonía invierten el camino: primero se mide una característica del sistema, y de esa medición salen las ganancias con una fórmula. No dan el ajuste óptimo, pero sí uno sensato, con el que se parte y se refina.

La segunda mitad del tema es el filtrado. Los sensores reales entregan ruido (la lectura tiembla), cuantización (cuentas enteras) y, a veces, espigas (una lectura absurda de vez en cuando). La derivada del PID amplifica todo eso. Filtrar ayuda, pero cada filtro retrasa la señal, y un retardo extra en el lazo es lo que desestabiliza (tema 41). El tema muestra los números de ese compromiso.

La planta es el lazo de velocidad de un motor: una constante de tiempo de 0,5 s, una ganancia de 2 rad/s por unidad de orden, y un retardo de 0,1 s (el driver, el muestreo y la medida). La orden está limitada a ±1, el muestreo es de 10 ms, y a los 5 s aparece una carga constante.

Todos los números del tema salen de simulaciones reproducibles (el código está al final), y el sketch de Arduino se probó compilándolo en una computadora con aritmética de 32 bits.

02 · Medir antes de ajustar

Un modelo para sintonizar

Muchos procesos (un motor con carga, un horno, un tanque) se parecen, vistos desde afuera, a un sistema de primer orden con retardo: cuando se cambia la orden, nada ocurre durante un tiempo, y después la salida se acerca exponencialmente a un valor final. Tres números lo describen:

• K, la ganancia: cuánto cambia la salida por cada unidad de cambio de la orden, una vez asentada.
• L, el retardo (tiempo muerto): cuánto tarda en empezar a responder.
• T, la constante de tiempo: qué tan rápido se acerca al valor final, una vez que empezó.

Se miden con un experimento de lazo abierto: se parte del reposo, se aplica una orden constante (aquí 0,5, lejos de la saturación) y se registra la salida. Con dos puntos de la curva se obtienen L y T: el tiempo en que la salida llega al 28,3 % y al 63,2 % de su valor final.

Las fórmulas: T = 1,5 · (t63 − t28) y L = t63 − T. Con los tiempos de la figura, T = 1,5 · (0,598 − 0,266) = 0,498 s y L = 0,598 − 0,498 = 0,100 s; la ganancia es el cambio de la salida dividido por el de la orden, K = 1,0 / 0,5 = 2,0 (1,999 medido). El modelo identificado coincide con el real (K = 2, T = 0,5 s, L = 0,1 s), porque el simulador es exactamente de ese tipo. En un sistema real, el modelo es una aproximación y los números salen un poco distintos cada vez.

Dos precauciones para el experimento. Debe hacerse con el sistema en un punto de trabajo normal y sin perturbaciones que se muevan durante la prueba. Y debe repetirse con un escalón de otro tamaño: si K cambia mucho, el sistema no es lineal (saturación, zona muerta), y las ganancias calculadas valdrán sólo cerca de ese punto.

Una cantidad útil es la relación L / T (aquí 0,2). Mide qué tan «difícil» es controlar el sistema: con L/T pequeña, el retardo casi no molesta y se admiten ganancias altas; con L/T cerca de 1 o más, el retardo domina y hay que ir con cuidado.

03 · Lo que se sabe sin fórmulas

Ajuste manual en orden

Antes de las fórmulas, el resumen de lo que hace subir cada ganancia, que se verificó en los temas 41 a 43:

El orden recomendado en el tema 43 sigue siendo válido: Kp para la rapidez, Kd para amortiguar, Ki al final y poco. Para este lazo, un P solo con Kp = 2 no tiene sobreimpulso, pero la carga lo deja con un error permanente de 0,28 rad/s y una caída de 0,307: ese es el punto de partida contra el que se miden los métodos que siguen.

El ajuste manual tiene su lugar: es el que se usa para refinar lo que dan las fórmulas, mirando la respuesta real. Lo que las fórmulas aportan es un punto de partida que no requiere adivinar.

04 · La curva de reacción

Ziegler–Nichols por curva de reacción

Ziegler y Nichols publicaron en 1942 reglas empíricas basadas en K, L y T. Se diseñaron para que el lazo cerrado rechace bien las perturbaciones, con una respuesta que decae en cada oscilación a un cuarto de la anterior. Las fórmulas para el modelo identificado:

Las ganancias en la forma que usa el código se obtienen de Ki = Kp / Ti y Kd = Kp · Td. Con T = 0,498, K = 1,999 y L = 0,100, la base es T / (K · L) = 2,49: un sistema con más retardo (L mayor) o más ganancia (K mayor) recibe menos Kp, como corresponde. Un P puro con 2,49 dejaría un error permanente de 1 / (1 + K · Kp) = 17 % en un escalón.

La ventaja de este método es que no obliga a llevar el sistema cerca de la inestabilidad. La desventaja es que la identificación con dos puntos es sensible a la precisión con que se lee la curva, y que está pensado para sistemas con L/T entre 0,1 y 1.

05 · La ganancia crítica

Ziegler–Nichols por oscilación sostenida

La otra variante no necesita un modelo: se hace en lazo cerrado con un P puro. Se sube Kp de a poco hasta que la salida oscila de forma sostenida, ni creciente ni decreciente. Esa es la ganancia crítica Ku, y el período de la oscilación es Tu.

En el simulador, la oscilación sostenida aparece con Ku = 4,06, con un período Tu = 0,390 s. El Ku teórico de la planta continua es 4,25; el muestreo de 10 ms agrega, en efecto, medio período de retardo y lo baja a 4,06. La prueba reveló el efecto del muestreo sin que nadie lo modelara, que es justamente lo valioso del método.

Cuidados indispensables en un sistema real:

• Llevar un sistema a oscilar puede dañarlo. Hay que usar una referencia pequeña, límites de la salida que contengan la amplitud, y estar listo para cortar.
• Con saturación, la oscilación sostenida es en realidad un ciclo límite, y su amplitud depende de los límites, no de Ku. Lo que importa para Ku es en qué ganancia deja de amortiguarse.
• Ku depende fuertemente del retardo. Con el mismo motor, al variar el retardo del lazo:

Duplicar el retardo reduce Ku a la mitad aproximadamente y duplica Tu. Esto da una regla para después: un retardo extra en el lazo (por ejemplo, el de un filtro) baja la ganancia máxima admisible; y las ganancias ya calculadas pueden pasarse del nuevo límite.

06 · Probar las reglas

Comparar las reglas

Las reglas anteriores y una más conservadora, la de SIMC (Skogestad internal model control): para un PI, Kp = T / (2 · K · L) y Ti = mín(T, 8 · L). Está pensada para no sobrepasar. Resultados con el escalón de 1 rad/s y la carga a los 5 s:

Qué se aprende de la tabla:

• Todas las reglas con integral anulan el error permanente que deja el P solo. Eso es lo que comprobaba el tema 43, ahora sin tanteo.
• Ziegler–Nichols rechaza mejor la carga: caída de 0,080 y recuperación en 0,24 s, contra 0,119 y 1,13 s de SIMC. Para eso fue diseñado. SIMC, en cambio, es lo más suave: sin sobreimpulso y con la orden casi sin ruido.
• Las dos variantes de Z-N dan resultados parecidos (3,1 % y 5,5 % de sobreimpulso), aunque se obtuvieron de experimentos distintos: una prueba en lazo abierto y otra en lazo cerrado. Es una buena señal de consistencia.
• La derivada tiene un costo de ruido: el desvío de la orden pasa de 0,012 (PI) a 0,140 (PID) con sólo 5 mrad/s de ruido. Es un factor 12. Si el ruido no se resuelve, la ventaja del PID en rechazo de carga se paga con un actuador que tiembla.
• Z-N en este sistema da poco sobreimpulso. En otros, con más retardo relativo, el sobreimpulso típico es de 25 % o más: por eso se dice que sus ganancias son un punto de partida, y se las suele reducir un 20–30 % si se quiere más suavidad.

07 · Limpiar la medida

Filtrar la medida

Una medición real combina la señal con tres tipos de problemas:

• Ruido de banda ancha: la lectura tiembla alrededor del valor verdadero (el ruido eléctrico del sensor, la cuantización del encoder).
• Espigas: lecturas absurdas aisladas (un eco falso del ultrasónico, una interferencia). Muy distintas del ruido: una sola lectura vale mucho más que las demás.
• Deriva y offset: se resuelven calibrando, no filtrando.

Para los dos primeros hay dos filtros muy sencillos, que se programan en pocas líneas:

Filtro exponencial

Es un pasabajos de primer orden. Cada muestra nueva se mezcla con la salida anterior: y = y + α · (x − y), con α = Ts / (tf + Ts). El parámetro es la constante de tiempo tf: cuanto mayor, más suave y más lento. Con Ts = 10 ms, tf = 50 ms da α = 0,167; y α = 0,1 corresponde a tf = Ts · (1 − α) / α = 90 ms. Pesa muy poco en memoria (una variable) y en cálculo (una multiplicación). Su gran desventaja es el retardo: la salida llega, aproximadamente, tf segundos después que la entrada.

Filtro de mediana

La mediana de las últimas tres lecturas descarta una espiga aislada, porque queda siempre en un extremo. No suaviza el ruido continuo. Su retardo es mínimo (una muestra) y no distorsiona las variaciones reales: un escalón pasa casi sin cambios, apenas una muestra después. Es el que se usó en el ultrasónico del tema 46.

Los dos se pueden combinar: primero la mediana para sacar las espigas y luego el exponencial para suavizar. Una tercera opción, la media móvil de N muestras, suaviza como el exponencial, con un retardo de (N − 1) · Ts / 2, pero necesita guardar N valores. La figura compara el exponencial con tres constantes de tiempo sobre una señal ruidosa:

Cuanto más suave la curva, más tarde empieza a subir. Esa es la ley del filtrado: no se puede quitar ruido sin agregar retardo. La sección siguiente mide cuánto cuesta ese retardo dentro del lazo.

08 · El otro lado del filtro

El precio del filtro

El filtro está dentro del lazo: la medida filtrada es la que el PID compara con la referencia. Su retardo se suma al del sistema y consume el margen de estabilidad. Con el PID de Z-N (2,99; 14,95; 0,149) y un ruido de 10 mrad/s en la medida, al variar la constante del filtro:

La tabla tiene una forma muy clara:

• Hasta 20 ms (una quinta parte del retardo L = 100 ms) el filtro es casi gratis: el ruido en la orden cae a la tercera parte (de 0,231 a 0,066), y el sobreimpulso sube sólo de 3,1 a 5,5 %.
• A 50 ms el ruido en la orden es de 0,032, pero el asentamiento se alarga a 1,38 s y la carga tarda 1,03 s en recuperarse: se pagó el ruido con respuesta.
• A partir de 100 ms (igual al retardo del sistema) la cosa se rompe: el sobreimpulso es de 17 %, el sistema no se asienta, y el desvío de la orden vuelve a subir, porque ahora oscila. Más filtro dejó más ruido en la orden, no menos.

La explicación está en la ganancia crítica. El filtro de 100 ms agrega, aproximadamente, 100 ms de retardo al lazo, que pasa de 100 a unos 200 ms equivalentes. Para ese retardo la tabla de la sección 05 da Ku ≈ 2,24; las ganancias de Z-N (Kp = 2,99) superan ese límite, y el lazo queda al borde de la inestabilidad. Regla práctica: la constante del filtro debe ser mucho menor que el retardo del lazo (aquí, menos de L/2, y mejor L/5); si el ruido requiere más filtro, hay que bajar las ganancias.

Espigas: la mediana, no el exponencial

Con 1 % de las lecturas valiendo +1,5 rad/s (espigas) y un ruido pequeño de 2 mrad/s:

El filtro exponencial no sirve contra las espigas: las reparte en el tiempo en lugar de eliminarlas, y la orden sigue dando saltos de hasta 2,0 (el rango completo, de −1 a +1) con una espiga; encima empeora el sobreimpulso. La mediana de tres las elimina por completo (el salto máximo baja a 0,11). La lección es elegir el filtro según el problema: espigas, mediana; ruido continuo, exponencial; los dos, mediana primero.

La derivada y el ruido

La derivada es la que más sufre. Con el ruido de 10 mrad/s, Kp = 3 y Ki = 15 sin filtrar, el desvío de la orden según Kd:

Resultado curioso: sin derivada el desvío no es el menor (0,161). Con Kp = 3 y Ki = 15 el lazo está poco amortiguado y resuena, amplificando el ruido de su propia frecuencia. Una pizca de derivada (0,05) lo amortigua y baja el desvío a 0,098; más derivada ya sólo amplifica ruido. El mínimo de ruido está en una Kd chica, no en cero ni en el valor de Z-N. Y la derivada debe aplicarse siempre a la medida filtrada, como hacen el código y la figura.

09 · Poner todo en orden

Receta completa de ajuste

1. Escribí las especificaciones (tema 41): sobreimpulso máximo, tiempo de asentamiento, caída tolerable ante la carga, ruido admisible en la orden.
2. Identificá el modelo con un escalón en lazo abierto, dentro de la zona lineal. Repetí con otro tamaño y comprobá que K no cambia.
3. Calculá las ganancias con una regla: Z-N si importa rechazar perturbaciones; SIMC o Z-N reducido si importa no sobrepasar. Empezá por un PI y agregá la derivada sólo si hace falta.
4. Probá con el sistema real con la referencia y con la carga más desfavorable. Medí sobreimpulso, asentamiento, caída y ruido de la orden.
5. Revisá el ruido. Mirá el desvío de la orden en reposo. Si el actuador tiembla, filtrá: mediana para espigas, exponencial para ruido continuo, con tf de una quinta parte del retardo.
6. Volvé a medir con el filtro puesto. El filtro cambia la sintonía: lo que valía antes puede ya no valer. Si el sobreimpulso crece, bajá las ganancias un 20–30 %.
7. Refiná a mano en el orden del tema 43, un cambio por vez, anotando cada resultado.
8. Probá los extremos: batería baja, carga máxima, escalón grande. Un ajuste sólo vale en el rango en que se probó.

El tema 48 aplica esta receta a los lazos del robot.

10 · Código para la placa

Sintonía y filtro en Arduino

El sketch identifica el modelo con un escalón en lazo abierto sobre el motor simulado, calcula las ganancias de Ziegler–Nichols de curva de reacción y ensaya el PID con cuatro constantes del filtro exponencial. El filtro es una estructura de cuatro líneas, y el resto es el PID del tema 43. Con un motor real, la identificación se haría con los datos del sensor.

// Sintonía de un PID por Ziegler-Nichols de lazo abierto y filtro de la medida, sobre un motor simulado.
// Motor: TAU * w' = -w + K * (u(t - L) - carga). El sketch identifica K, L y T con un escalón y calcula las ganancias.
const float DT = 0.001, TAU = 0.5, K_PLANTA = 2.0, CARGA = 0.2, T_CARGA = 5.0;
const int N_COLA = 101;                       // retardo de 0,1 s = 100 pasos de 1 ms (se guarda uno más)
const long N_RET = 100;
float colaU[N_COLA];

// Un paso de 1 ms de la planta con retardo. k es el número de paso; carga en N·m equivalentes (orden).
float pasoPlanta(float w, float u, long k, float carga) {
  colaU[k % N_COLA] = u;
  float ua = (k >= N_RET) ? colaU[(k - N_RET) % N_COLA] : 0.0;
  return w + DT * (-w + K_PLANTA * (ua - carga)) / TAU;
}

// Filtro exponencial de primer orden: alfa = Ts / (tf + Ts).
struct FiltroEMA {
  float tf, y = 0;
  bool listo = false;
  FiltroEMA(float t) : tf(t) {}
  float filtrar(float x, float ts) {
    if (!listo || tf <= 0) { y = x; listo = true; return y; }
    y += ts / (tf + ts) * (x - y);
    return y;
  }
};

struct PID {
  float kp, ki, kd;
  float integral = 0, previa = 0;
  bool primera = true;
  PID(float p, float i, float d) : kp(p), ki(i), kd(d) {}
  float calcular(float ref, float medida, float dt) {
    float e = ref - medida;
    float d = primera ? 0.0 : -(medida - previa) / dt;
    primera = false; previa = medida;
    float nueva = integral + ki * e * dt;
    float salida = kp * e + nueva + kd * d;
    bool satura = salida > 1.0 || salida < -1.0;
    if (!satura || (salida > 1.0 && e < 0) || (salida < -1.0 && e > 0)) integral = nueva;
    return constrain(salida, -1.0, 1.0);
  }
};

// Respuesta al escalón en lazo abierto con u = 0,5; dos puntos (28,3 % y 63,2 %) dan K, L y T.
void identificar(float& k, float& l, float& t) {
  float w = 0;
  for (long i = 0; i < 4000; i++) w = pasoPlanta(w, 0.5, i, 0);
  float fin = w;
  w = 0;
  float t28 = -1, t63 = -1;
  for (long i = 0; i < 4000; i++) {
    w = pasoPlanta(w, 0.5, i, 0);
    if (t28 < 0 && w >= 0.283 * fin) t28 = i * DT;
    if (t63 < 0 && w >= 0.632 * fin) t63 = i * DT;
  }
  t = 1.5 * (t63 - t28);
  l = t63 - t;
  k = fin / 0.5;
}

// Ziegler-Nichols de lazo abierto (curva de reacción), PID: Kp = 1,2 T / (K L), Ti = 2 L, Td = L / 2.
void zieglerNicholsPID(float k, float l, float t, float& kp, float& ki, float& kd) {
  kp = 1.2 * t / (k * l);
  ki = kp / (2.0 * l);
  kd = kp * 0.5 * l;
}

// Escalón de 1 rad/s y carga a los 5 s, con filtro de constante tf sobre la medida. Control cada 10 ms.
void ensayar(float kp, float ki, float kd, float tf) {
  for (int i = 0; i < N_COLA; i++) colaU[i] = 0;
  PID pid(kp, ki, kd);
  FiltroEMA filtro(tf);
  float w = 0, u = 0, maximo = 0, minimo = 10;
  long ultimo = -1;
  for (long k = 0; k < 10000; k++) {
    if (k % 10 == 0) u = pid.calcular(1.0, filtro.filtrar(w, 0.01), 0.01);
    w = pasoPlanta(w, u, k, k * DT >= T_CARGA ? CARGA : 0.0);
    if (k < 5000) {
      maximo = max(maximo, w);
      if (fabs(w - 1.0) > 0.02) ultimo = k;
    } else {
      minimo = min(minimo, w);
    }
  }
  Serial.print("tf="); Serial.print(tf * 1000, 0);
  Serial.print(" ms  Mp="); Serial.print(max(0.0f, (maximo - 1.0f) * 100.0f), 2); Serial.print(" %");
  if (ultimo >= 4999) Serial.print("  ts=no se asienta");
  else { Serial.print("  ts="); Serial.print((ultimo + 1) * DT, 2); Serial.print(" s"); }
  Serial.print("  caida="); Serial.println(1.0 - minimo, 3);
}

void setup() {
  Serial.begin(115200);
  float k, l, t, kp, ki, kd;
  identificar(k, l, t);
  Serial.print("K="); Serial.print(k, 3); Serial.print("  L="); Serial.print(l, 3);
  Serial.print(" s  T="); Serial.print(t, 3); Serial.println(" s");
  zieglerNicholsPID(k, l, t, kp, ki, kd);
  Serial.print("Kp="); Serial.print(kp, 2); Serial.print("  Ki="); Serial.print(ki, 2);
  Serial.print("  Kd="); Serial.println(kd, 3);
  const float filtros[] = {0.0, 0.02, 0.05, 0.1};
  for (byte i = 0; i < 4; i++) ensayar(kp, ki, kd, filtros[i]);
}

void loop() {}

La salida repite las tablas del texto:

K=1.999  L=0.100 s  T=0.498 s
Kp=2.99  Ki=14.95  Kd=0.149
tf=0 ms  Mp=3.12 %  ts=0.80 s  caida=0.080
tf=20 ms  Mp=5.45 %  ts=0.66 s  caida=0.088
tf=50 ms  Mp=8.04 %  ts=1.38 s  caida=0.096
tf=100 ms  Mp=17.16 %  ts=no se asienta  caida=0.176

Se compiló y ejecutó en una computadora con aritmética de 32 bits (la misma de la placa) y reproduce los resultados de Python. Con un sensor real, el filtro se llama en cada ciclo de control con la lectura nueva, y su salida es la que entra al PID: u = pid.calcular(ref, filtro.filtrar(lectura, Ts), Ts). Como en el tema 46, para las espigas se antepone la mediana de tres.

11 · Exploración interactiva

Laboratorios de sintonía

Probá este recorrido: cargá Z-N lazo abierto PID y poné un ruido de 10 mrad/s (desvío de la orden de 0,23); subí el filtro a 20 ms (0,07) y a 50 ms (0,03); seguí hasta 100 ms y mirá cómo se rompe la respuesta. Después activá las espigas y compará exponencial contra mediana. Por último, cargá SIMC y mirá la diferencia de rechazo de la carga.

Subí Kp desde 3,5 hasta que el cociente pase de 1: ocurre cerca de 4,06. Con ganancias mayores la oscilación crece hasta que la orden satura y queda un ciclo límite. Medí el período entre dos picos de la curva: unos 0,39 s.

12 · Código ejecutable

Todo el ajuste en Python

El simulador del motor con el PID filtrado, la identificación, la búsqueda de la ganancia crítica, las cuatro reglas de sintonía, y las experiencias con filtros, espigas y ruido. El código reproduce las tablas del tema.

from math import sqrt

# Lazo de velocidad de un motor: TAU*w' = -w + K*u(t - L). Orden limitada a [-1, 1]; carga de 0,2 a los 5 s.
DT, TAU, K, L, UMAX, CARGA, T_CARGA = 0.001, 0.5, 2.0, 0.1, 1.0, 0.2, 5.0

def azar(semilla):
    """Generador pseudoaleatorio de 32 bits (mulberry32): el ruido es repetible."""
    a = semilla & 0xFFFFFFFF
    def siguiente():
        nonlocal a
        a = (a + 0x6D2B79F5) & 0xFFFFFFFF
        t = ((a ^ (a >> 15)) * (1 | a)) & 0xFFFFFFFF
        t = (((t + (((t ^ (t >> 7)) * (61 | t)) & 0xFFFFFFFF)) & 0xFFFFFFFF) ^ t)
        return ((t ^ (t >> 14)) & 0xFFFFFFFF) / 4294967296.0
    return siguiente

def simular(kp=1.0, ki=0.0, kd=0.0, tf=0.0, ts=0.01, ref=1.0, ruido=0.0, retardo=L, t_fin=10.0,
            espigas=0.0, mediana=False, u_fijo=None, carga=CARGA, semilla=7):
    """PID con filtro exponencial (constante de tiempo tf) sobre la medida; derivada sobre la medida."""
    rnd = azar(semilla)
    w = 0.0
    u = 0.0 if u_fijo is None else u_fijo
    integral = fil = fil_prev = m1 = m2 = 0.0
    primera = True
    nc, nd = max(1, round(ts / DT)), round(retardo / DT)
    ub, W, U = [], [], []
    for k in range(round(t_fin / DT)):
        t = k * DT
        if u_fijo is None and k % nc == 0:
            m = w + ruido * (sum(rnd() for _ in range(12)) - 6.0)
            if espigas > 0 and rnd() < espigas:
                m += 1.5                                         # una espiga: lectura errónea ocasional
            if mediana:
                a, b = (m, m) if primera else (m1, m2)
                med = max(min(m, a), min(max(m, a), b))
                m2, m1, m = m1, m, med
            if primera:
                fil = fil_prev = m
            elif tf > 0:
                fil += ts / (tf + ts) * (m - fil)
            else:
                fil = m
            e = ref - fil
            d = 0.0 if primera else -(fil - fil_prev) / ts
            primera, fil_prev = False, fil
            nueva = integral + ki * e * ts
            raw = kp * e + nueva + kd * d
            sat = raw > UMAX or raw < -UMAX
            u = min(UMAX, max(-UMAX, raw))
            if not sat or (raw > UMAX and e < 0) or (raw < -UMAX and e > 0):
                integral = nueva
        ub.append(u)
        ua = ub[k - nd] if k - nd >= 0 else 0.0
        w += DT * (-w + K * (ua - (carga if t >= T_CARGA else 0.0))) / TAU
        W.append(w); U.append(u)
    return W, U

def escalon(w, ref=1.0):
    """Sobreimpulso (%) y asentamiento al 2 % durante los 5 s previos a la carga."""
    n = round(T_CARGA / DT)
    ult = max((i for i in range(n) if abs(w[i] - ref) > 0.02 * ref), default=-1)
    return max(0.0, (max(w[:n]) - ref) / ref * 100.0), (None if ult >= n - 1 else (ult + 1) * DT)

def respuesta_carga(w, ref=1.0):
    """Caída máxima y tiempo de recuperación tras la carga."""
    i0 = round(T_CARGA / DT)
    ult = max((i for i in range(i0, len(w)) if abs(w[i] - ref) > 0.02 * ref), default=-1)
    return max(0.0, ref - min(w[i0:])), (0.0 if ult < i0 else None if ult >= len(w) - 1 else (ult + 1 - i0) * DT)

def desvio(x, a, b):
    seg = x[round(a / DT):round(b / DT)]
    m = sum(seg) / len(seg)
    return sqrt(sum((v - m) ** 2 for v in seg) / len(seg))

def identificar(w, u_esc):
    """Dos puntos (28,3 % y 63,2 %) de una respuesta al escalón en lazo abierto: devuelve K, L y T."""
    fin = w[-1]
    t28 = next(i for i, v in enumerate(w) if v >= 0.283 * fin) * DT
    t63 = next(i for i, v in enumerate(w) if v >= 0.632 * fin) * DT
    T = 1.5 * (t63 - t28)
    return fin / u_esc, t63 - T, T

def crecimiento(kp, retardo=L):
    """Cociente de amplitudes de la oscilación en dos ventanas: < 1 se amortigua, > 1 crece."""
    w, _ = simular(kp=kp, ref=1e-4, carga=0.0, t_fin=5.0, retardo=retardo)
    pp = lambda a, b: max(w[round(a / DT):round(b / DT)]) - min(w[round(a / DT):round(b / DT)])
    return pp(3, 4) / pp(2, 3)

def ganancia_critica(retardo=L, lo=3.5, hi=5.0, pasos=18):
    """Ku por bisección y Tu (período de la oscilación) midiendo los picos."""
    for _ in range(pasos):
        medio = (lo + hi) / 2
        lo, hi = (medio, hi) if crecimiento(medio, retardo) < 1 else (lo, medio)
    w, _ = simular(kp=lo * 1.0005, ref=1e-4, carga=0.0, t_fin=8.0, retardo=retardo)
    picos = [i * DT for i in range(3000, 7990) if w[i] > w[i - 1] and w[i] >= w[i + 1]]
    return lo, (picos[-1] - picos[0]) / (len(picos) - 1)

# Reglas de sintonía (devuelven kp, ki, kd)
def zn_abierto_pi(k, l, t):  kp = 0.9 * t / (k * l);  return kp, kp / (3.33 * l), 0.0
def zn_abierto_pid(k, l, t): kp = 1.2 * t / (k * l);  return kp, kp / (2 * l), kp * 0.5 * l
def zn_cerrado_pid(ku, tu):  return 0.6 * ku, 0.6 * ku / (tu / 2), 0.6 * ku * tu / 8
def simc_pi(k, l, t):        kp = t / (2 * k * l);    return kp, kp / min(t, 8 * l), 0.0

def f(v, d=2):
    return "  —  " if v is None else f"{v:.{d}f}"

if __name__ == "__main__":
    w, _ = simular(u_fijo=0.5, carga=0.0, t_fin=4.0)
    k, l, t = identificar(w, 0.5)
    print(f"Identificación en lazo abierto (escalón de 0,5): K = {k:.3f}, L = {l:.3f} s, T = {t:.3f} s")
    ku, tu = ganancia_critica()
    print(f"Prueba de oscilación: Ku = {ku:.2f}, Tu = {tu:.3f} s")

    print("\nReglas de sintonía (escalón de 1 rad/s y carga de 0,2 a los 5 s; ruido de 5 mrad/s para el desvío de u)")
    print("regla       Kp     Ki     Kd   Mp(%)  ts(s)  caída  recup.  desvío u")
    reglas = [("solo P", (2.0, 0.0, 0.0)), ("ZN lazo ab. PI", zn_abierto_pi(k, l, t)),
              ("ZN lazo ab. PID", zn_abierto_pid(k, l, t)), ("ZN lazo cer. PID", zn_cerrado_pid(ku, tu)),
              ("SIMC PI", simc_pi(k, l, t))]
    for nombre, (kp, ki, kd) in reglas:
        w, _ = simular(kp=kp, ki=ki, kd=kd)
        mp, ts = escalon(w)
        caida, rec = respuesta_carga(w)
        _, u = simular(kp=kp, ki=ki, kd=kd, ruido=0.005)
        print(f"{nombre:<17}{kp:5.2f} {ki:6.2f} {kd:6.3f} {mp:6.2f} {f(ts):>6} {caida:6.3f} {f(rec):>6}   {desvio(u, 7, 10):6.3f}")

    print("\nRuido de 10 mrad/s con el PID Kp = 3, Ki = 15 y distintas Kd (sin filtro):")
    for kd_ in (0.0, 0.05, 0.1, 0.15, 0.3):
        _, u = simular(kp=3.0, ki=15.0, kd=kd_, ruido=0.01)
        print(f"  Kd = {kd_:4.2f}: desvío de u {desvio(u, 7, 10):.3f}")

    print("\nGanancia crítica según el retardo del lazo:")
    for retardo in (0.05, 0.10, 0.15, 0.20):
        ku_, tu_ = ganancia_critica(retardo, 0.5, 10.0, 22)
        print(f"  retardo {retardo * 1000:3.0f} ms: Ku = {ku_:.2f}, Tu = {tu_:.3f} s")

    print("\nFiltro exponencial sobre la medida con el PID de Ziegler–Nichols (ruido de 10 mrad/s)")
    print("  tf(ms)  alfa   Mp(%)  ts(s)  caída  recup.  desvío u")
    kp, ki, kd = zn_abierto_pid(k, l, t)
    for tf in (0.0, 0.01, 0.02, 0.05, 0.1, 0.2):
        w, _ = simular(kp=kp, ki=ki, kd=kd, tf=tf)
        mp, ts = escalon(w)
        caida, rec = respuesta_carga(w)
        _, u = simular(kp=kp, ki=ki, kd=kd, tf=tf, ruido=0.01)
        print(f"  {tf * 1000:5.0f}  {0.01 / (tf + 0.01):5.3f} {mp:6.2f} {f(ts):>6} {caida:6.3f} {f(rec):>6}   {desvio(u, 7, 10):6.3f}")

    print("\nEspigas (1 % de las lecturas valen +1,5 rad/s) con ruido de 2 mrad/s:")
    for nombre, kw in (("sin filtro", {}), ("exponencial 50 ms", dict(tf=0.05)),
                       ("mediana de 3", dict(mediana=True)), ("mediana + exp. 20 ms", dict(mediana=True, tf=0.02))):
        w, u = simular(kp=kp, ki=ki, kd=kd, espigas=0.01, ruido=0.002, **kw)
        salto = max(abs(u[i] - u[i - 1]) for i in range(7000, 10000))
        print(f"  {nombre:<22} desvío de u {desvio(u, 7, 10):.3f}, salto máximo de u {salto:.2f}, Mp {escalon(w)[0]:.2f} %")

La salida repite las tablas del texto: la identificación (K = 1,999; L = 0,100 s; T = 0,498 s), Ku = 4,06 y Tu = 0,390 s, las cinco sintonías, el ruido según Kd, la ganancia crítica según el retardo, la tabla de filtros y la de espigas.

Probá cambiar el retardo L a 0,2 y recalcular todo: Z-N pide ganancias menores y las cuentas siguen valiendo. O fijá tf = L / 2 en el PID de Z-N y observá cuánto se degrada la respuesta.

13 · Diagnóstico

Errores frecuentes

1. Identificar con un escalón demasiado grande y llegar a la saturación: K y T medidos no son los del sistema lineal.
2. Aplicar Z-N y dejar las ganancias sin revisar. Son un punto de partida; en sistemas con más retardo relativo, el sobreimpulso puede ser del 25 % o más.
3. Llevar un sistema real a la oscilación sostenida sin límites ni cuidado.
4. Filtrar demasiado para «tener una señal bonita» y pasar el límite de estabilidad.
5. Usar un filtro exponencial contra espigas, o una mediana contra ruido continuo.
6. Ajustar las ganancias sin filtro y después agregar el filtro sin volver a probar.
7. Derivar la señal sin filtrar, o derivar el error en lugar de la medida (tema 43).

14 · Práctica guiada

Predecí antes de comprobar

1. Identificar

Un escalón de orden 0,5 lleva la velocidad a 1,0 rad/s. Pasa por el 28,3 % a los 0,266 s y por el 63,2 % a los 0,598 s. ¿Cuáles son K, L y T?

2. Otra planta

Otro sistema tiene K = 1,5, L = 0,2 s y T = 1,0 s. ¿Qué Kp, Ki y Kd da Z-N de lazo abierto para un PID?

3. Lazo cerrado

Con Ku = 4,06 y Tu = 0,390 s, ¿qué da Z-N de lazo cerrado para un PI?

4. Constante del filtro

Con Ts = 10 ms, ¿qué α corresponde a tf = 50 ms? ¿Y qué tf a α = 0,1?

5. Por qué se rompe con 100 ms

El filtro de 100 ms deja al PID de Z-N (Kp = 2,99) sin asentarse. Con la tabla de Ku según el retardo, ¿cómo lo explicás?

6. Elegir el filtro

Un sensor entrega ruido pequeño y, una de cada cien veces, una lectura disparatada. ¿Qué filtro usás y cuánto cambia la orden?

7. La derivada mínima

¿Por qué el desvío de la orden es menor con Kd = 0,05 que con Kd = 0 en la tabla del ruido?

15 · Síntesis

Ideas para recordar

¿Listo para el tema 45?

Antes de avanzar, identificá K, L y T de uno de los motores de tu robot con un escalón de orden 0,5 y calculá las ganancias Z-N de un PI. Probalas con y sin un filtro exponencial de 20 ms y anotá el sobreimpulso y el ruido de la orden en cada caso.

Con eso termina el bloque de control. En el próximo tema empieza la etapa final del curso: el proyecto de un robot completo, con sus requisitos, su diseño y su plan de pruebas.