Armónicos a la carta: eliminación selectiva por control óptimo

Publicado el

Mueve la amplitud, cambia los niveles o elige qué armónicos eliminar: cada cambio resuelve el problema de control óptimo en tu navegador. Abrir a pantalla completa ↗

Un convertidor de potencia (el que alimenta un motor desde una batería, o vuelca a la red la energía de unos paneles solares) no fabrica una sinusoide: solo sabe conmutar entre unos pocos niveles de tensión. Lo que sale es una onda escalonada, y toda onda escalonada lleva, además de la sinusoide que se quería, una colección de armónicos: componentes a 3, 5, 7… veces la frecuencia fundamental. En un motor trifásico, el 5.º y el 7.º calientan los devanados y hacen vibrar el par; en la red, los limitan las normas de calidad.

La eliminación selectiva de armónicos (SHE, selective harmonic elimination) es la técnica clásica para esto: elegir los instantes de conmutación de forma que la fundamental tenga la amplitud pedida y unos armónicos concretos salgan exactamente cero. Arriba está resuelta en directo para los armónicos 5, 7, 11 y 13 (los triples, 3, 9, 15…, se cancelan solos entre las fases de un sistema trifásico). Mueve la amplitud: la forma de la onda, el número de conmutaciones y dónde caen salen solos.

Esto es la versión divulgativa de un trabajo de 2020–2021 con Umberto Biccari y Carlos Esteve-Yagüe, dentro del proyecto CONVADP de electrónica de potencia de la Cátedra de Matemática Computacional (Deusto): el artículo Multilevel Selective Harmonic Modulation via Optimal Control (Applied Mathematics and Optimization, 2022; arXiv) y las diapositivas que lo acompañaban. El artículo resolvía el problema con IPOPT sobre una malla; el visor usa su formulación dual, que publicaron después Biccari y Zuazua, resuelta con Newton para que quepa en milisegundos.

El problema, y por qué el método clásico se atasca

La onda , en un periodo , toma valores en un conjunto finito de niveles (por ejemplo ) y tiene simetría de media onda, , así que solo hay armónicos impares. Sus coeficientes de Fourier son

y el problema es conseguir unos valores dados: (la amplitud de la fundamental), y para los armónicos que se quieren eliminar.

El planteamiento clásico fija primero la secuencia de niveles (de 0 sube a ½, luego a 1, baja a ½…) y busca después los ángulos de conmutación con Newton: con la secuencia fijada, los coeficientes son sumas de senos y cosenos de los ángulos, y es un sistema de ecuaciones no lineales. El problema es elegir la secuencia. Con cinco niveles y cinco conmutaciones por cuarto de periodo, saltando siempre a un nivel contiguo, hay 18 secuencias posibles; para una amplitud de 0,3, solo 6 tienen solución, y entre todas suman 8 soluciones distintas; para 0,8, solo 2 secuencias, con 4 soluciones. El método clásico tiene que probarlas, y sin buenos arranques Newton ni siquiera encuentra todas.

La onda como trayectoria de un sistema dinámico

El cambio de punto de vista del artículo es leer las integrales de Fourier como un sistema dinámico. Si es lo que falta por acumular de cada coeficiente,

con los valores objetivo, entonces la onda resuelve el problema exactamente cuando la trayectoria llega a . Diseñar la onda es controlar ese sistema: un problema de control óptimo,

con una penalización que hay que elegir con cuidado. El control se permite continuo en , y es la penalización la que tiene que obligarlo a quedarse en los niveles.

La penalización que fuerza la escalera

El principio del máximo de Pontryagin dice cómo es el control óptimo: en cada instante, minimiza , donde es una función suave que sale del estado adjunto. Con la penalización cuadrática de siempre, , ese mínimo varía de forma continua con : la onda tomaría cualquier valor, no una escalera. Con solo puede tomar tres valores, , y .

La idea del artículo es tomar la parábola a trozos: la interpolación lineal de entre niveles consecutivos. Es convexa, y en sus vértices (los niveles) la pendiente salta; el mínimo de una función lineal más una convexa a trozos cae siempre en un vértice. Así que el control óptimo solo puede tomar los valores de los niveles, y la secuencia y los ángulos de conmutación los decide la optimización, no quien la plantea.

Detalle técnicoPor qué con un número par de niveles no funciona del todo

Si dos niveles consecutivos son simétricos respecto al cero (con seis niveles, y ), la parábola a trozos es plana entre ellos: vale lo mismo en cualquier punto intermedio. Ahí la penalización no distingue entre los niveles y los valores de en medio, y para amplitudes pequeñas el óptimo deja de ser una escalera. Con un número impar de niveles, que incluye el 0, no hay tramos planos. Por eso el visor solo ofrece 3, 5, 7 y 9 niveles. Con dos niveles, (la SHE clásica de dos niveles), la parábola a trozos es constante y el método no dice nada.

Diez números en vez de una onda: por qué cabe en el navegador

El artículo resolvía el problema como se suele hacer: discretizando la onda en una malla de 300 puntos y pasándole a IPOPT un problema de optimización de 300 incógnitas, con la penalización suavizada para que fuera derivable. Unos segundos por amplitud.

Pero el principio del máximo da más de lo que parece, y es lo que aprovecha la formulación dual de Biccari y Zuazua (IFAC, 2022; arXiv). La función que decide el nivel en cada instante es : una combinación de los senos y cosenos de los armónicos controlados, con un coeficiente por ecuación. Con la fundamental y cuatro armónicos (cada uno con seno y coseno) son diez números. Dado , la onda entera está determinada: el nivel cambia justo donde cruza un umbral, y esos cruces se calculan con la precisión de la máquina.

Y sale de maximizar una función cóncava de diez variables (el dual del problema), cuyo gradiente es precisamente el error en los coeficientes de Fourier. Lo que añade el visor es cómo resolverlo: Newton con conmutaciones exactas: el hessiano se calcula de forma exacta a partir de los puntos de conmutación, y partiendo de la solución de la amplitud anterior bastan unas pocas iteraciones. De 1 a 20 ms por solución, en JavaScript puro, sin WebAssembly.

Hay otra ventaja que no es de velocidad: como el problema es convexo, la brecha de dualidad (coste de la onda menos valor del dual) es cero en el óptimo y solo en el óptimo. Cada solución trae su propio certificado.

Tres mapas de color, uno por número de niveles (3, 5 y 7). El eje horizontal es el ángulo en medio periodo, de 0 a 180 grados, y el vertical la amplitud de la fundamental, de 0 a 1,25. Cada fila es la onda óptima para esa amplitud, coloreada por nivel: bandas que se ramifican y se estrechan, con pulsos que nacen y desaparecen. Por encima de 1,17 la franja es oscura: sin solución.
Todas las soluciones de un vistazo: cada fila es la onda (medio periodo) para una amplitud, eliminando 5, 7, 11 y 13. Los pulsos nacen y mueren al cambiar la amplitud, y el número de conmutaciones cambia con ellos; nadie lo ha fijado. Por encima de 1,17 no hay solución, ni con 3, ni con 5, ni con 7 niveles (ni con 9, que no está dibujado).

Contra qué se comprueba

Contra el método del artículo. Reimplementado tal cual en Python con CasADi (malla de 300 puntos, penalización suavizada, ε = 10⁻², IPOPT), y comparado con el dual en 27 casos: 3, 5 y 7 niveles y nueve amplitudes entre −0,8 y 0,8. El coste coincide hasta la sexta cifra (el del dual siempre un pelo menor: es exacto en tiempo continuo, la malla no), la onda difiere en menos de un 1 % de su duración y la brecha de dualidad es cero a precisión de máquina. En Python, el dual tarda 0,2–0,4 s por caso e IPOPT 3–13 s.

Contra el método clásico. Para 5 niveles y amplitudes 0,3 y 0,8 se enumeraron las 18 secuencias de niveles posibles y, para cada una, se resolvieron las ecuaciones de los ángulos desde 60 arranques al azar. Entre todas las soluciones exactas encontradas, la que da el control óptimo es exactamente la de menor penalización (los mismos ángulos, hasta la centésima de grado), y la da sin probar ninguna secuencia.

JavaScript contra Python. Los mismos 27 casos con tres conjuntos de armónicos: misma estructura de conmutaciones en todos y ángulos iguales hasta 1,6 · 10⁻⁸ rad.

Que la eliminación es exacta. Con la penalización, el óptimo no anula los armónicos del todo: el error es proporcional a ε. Pero al bajar ε el vector converge, la onda también, y el error baja con él hasta ; el visor usa ε = 10⁻⁶. Cuando el objetivo es imposible, el dual lo dice: crece como y el error se queda donde está.

Dos gráficas en escala logarítmica frente a ε, de 0,1 a 10⁻⁸. A la izquierda, el error máximo en los armónicos: para amplitudes 0,3, 0,6 y 0,8 baja en línea recta con ε; para 1,2 se queda plano en torno a 0,016. A la derecha, el tamaño del multiplicador q: constante para las tres primeras, creciendo como 1/ε para 1,2.
Bajar ε separa los objetivos alcanzables de los imposibles: en los primeros el error cae como ε con q fijo; en el imposible, q explota y el error no baja.

Frente a los métodos habituales, a escala

En la literatura de SHE los ejemplos suelen ser pequeños: de 3 a 7 ángulos de conmutación por cuarto de periodo, es decir, eliminar de 2 a 6 armónicos. Los métodos algebraicos (resultantes, bases de Gröbner) se quedan en torno a 9 ángulos; las metaheurísticas (enjambres de partículas, algoritmos genéticos) se suelen enseñar con 3 o 5. Para ver dónde está cada uno, los tres métodos se midieron en la misma máquina, con el mismo código base y los mismos casos: la formulación clásica de tres niveles (cuarto de onda, solo senos), ángulos y los primeros armónicos impares no triples eliminados, en las amplitudes 0,2, 0,4, 0,6 y 0,8.

  • Newton desde ángulos al azar, repitiendo hasta acertar (con un máximo de 20 s).
  • Enjambre de partículas + Newton, con los parámetros de un trabajo reciente (Chen, Zhang y Yin, 2026): 300 partículas, 200 iteraciones.
  • El dual, sin ninguna información previa (en frío).
Ángulos Armónicos eliminadosNewton al azarEnjambre + NewtonDual (frío)
5hasta el 13.º4/4, < 1 ms3/4, 33–49 ms4/4, 4–32 ms
9hasta el 25.º4/4, < 11 ms2/4, 88–123 ms4/4, 14–26 ms
13hasta el 37.º4/4, 0,03–3,4 s1/4, ~180 ms4/4, 32–67 ms
17hasta el 49.º0/4 en 20 s0/4, ~290 ms4/4, 0,1–0,2 s
25hasta el 73.º0/4 en 20 s0/4, ~650 ms4/4, 0,3–0,6 s

Con pocos ángulos, Newton al azar es lo más rápido: el dual no aporta velocidad ahí. A partir de 13 ángulos la probabilidad de que Newton acierte desde un arranque al azar se desploma (miles de intentos), y con 17 o más no encuentra ninguna solución; el enjambre ya se queda corto con 9. El dual las encuentra todas, siempre con el número mínimo de ángulos ( = armónicos + 1, que nadie le ha dicho) y sin que nadie elija la secuencia ni el arranque.

Y una vez se tiene una solución, seguirla al cambiar la amplitud es otra historia: Newton, arrancando de la solución anterior en pasos de 0,01, tarda de 0,03 a 0,2 ms por paso, de 50 a 200 veces menos que el dual, y se mantiene en la misma rama (la de menor penalización). La combinación natural es esa: el dual para encontrar la primera solución, Newton para seguirla.

Detalle técnicoLo que hubo que arreglar por el camino
  • El arranque. En el dual tiene un punto anguloso y Newton no sabe salir; se arranca en la dirección del objetivo.
  • Los pulsos que nacen. Las conmutaciones se localizan muestreando en 2000 puntos y refinando con bisección. Justo cuando nace un pulso es más estrecho que una muestra y pasaba sin verse: el gradiente salía mal y Newton se clavaba (por ejemplo, en la amplitud 0,775 con cinco niveles). Ahora se localizan también los extremos de entre muestras, y ningún cruce queda escondido. Como comprobación independiente, los cruces se calcularon también sin muestrear: con , es un polinomio en y sus raíces sobre la circunferencia unidad son todos los cruces a la vez. En 2000 casos al azar coinciden hasta ; con muchos armónicos (polinomios de grado 100) las raíces son más lentas y delicadas que el muestreo, así que el visor muestrea.
  • Los pliegues del dual. Donde dos conmutaciones se juntan, el dual no es derivable; si Newton se atasca, se rescata bajando ε poco a poco o desde un objetivo ligeramente desplazado. En barridos completos de 0 a 1,25 con 3, 5, 7 y 9 niveles, todas las amplitudes alcanzables se resuelven.

Qué no demuestra esto

  • Óptimo para esta penalización, no para el convertidor. «Mejor» aquí quiere decir menor : la onda que se separa menos del cero, nivel a nivel. Un convertidor real puede preferir menos conmutaciones (pérdidas), o repartir de otra forma los niveles entre las celdas. Eso es otra penalización, y no siempre convexa.
  • La distorsión que queda. Se eliminan exactamente los armónicos elegidos; los demás quedan donde caigan. El visor muestra la distorsión total hasta el 49.º, que incluye los triples (que en una carga trifásica se cancelan entre fases).
  • No hay hardware. Ni tiempos muertos, ni retardos de conmutación, ni la tensión real de cada nivel: es la onda ideal.
  • Las ideas de base no son de aquí. La penalización lineal a trozos que fuerza valores discretos es un caso particular de penalizaciones anteriores del control con valores discretos (la SOAV de Ikeda y Nagahara, el control multi-bang de Clason y Kunisch), y el dual es el de Biccari y Zuazua. Lo propio de esta línea es aplicarlas a la SHE multinivel; lo de este artículo, el resolvedor y hasta dónde escala.
  • El banco de pruebas es uno. Una máquina, la formulación de tres niveles y cuatro amplitudes por tamaño; con otras implementaciones de Newton (mejores arranques, homotopía) o del enjambre, los números del método clásico mejorarían.
  • La existencia, a medias. El dual detecta cuándo un objetivo es imposible, pero no explica dónde está la frontera (1,17 para este conjunto de armónicos). Esa fue la otra mitad del trabajo, con programación dinámica, y queda para otro artículo.

Cómo reproducirlo

Todo está en control/she-control-optimo del repositorio de las simulaciones:

python python/she.py              # el dual en Python: residuo, brecha, conmutaciones
python python/compara_ipopt.py    # frente al método del artículo (CasADi + IPOPT)
python python/clasico.py 0.8      # frente al clásico: todas las secuencias de niveles
node web/prueba/compara.mjs       # JavaScript frente a Python
node web/prueba/banco.mjs         # Newton al azar, enjambre + Newton y el dual, N = 5…25
node web/prueba/continuacion.mjs  # seguir la solución: Newton frente al dual en caliente
cd web && python3 -m http.server  # el visor de arriba, en local

El código original del artículo (MATLAB y CasADi) está en SHE-Optimal-Control-paper.

Etiquetas: control óptimo · electrónica de potencia · Pontryagin · optimización convexa · control

← Todos los artículos