Predicción de voltajes en la red eléctrica por interpolación Kriging

Revista Iberoamericana de Autom´atica e Inform´atica Industrial 22 (2025) 112-119 www.revista-riai.org Resumen En este trabajo, abordamos el problema de la predicci´ on en l´ ınea de las trayectorias de voltaje e intensidad nodales en la red de distribuci´ on. Para esto, proponemos una formulaci´ on basada en datos utilizando la interpolaci´ on Kriging, una t´ ecnica de aprendizaje autom´ atico que ha mostrado aplicaciones prometedoras en el campo del control basado en datos. Producimos un or´ aculo de predicci´ on no param´ etrico que permite inferir trayectorias futuras directamente a partir de medidas de voltaje e intensidad en tiempo real. Adem´ as, proporcionamos una implementaci´ on algor´ ıtmica simple pero efectiva basada en el conocido esquema ISTA. Demostramos la efectividad de nuestra metodolog´ ıa para la predicci´ on r´ apida (subsegundos) de la din´ amica del voltaje mediante simulaciones. Palabras clave: Aprendizaje para el control, M´ etodos no param´ etricos, Redes el´ ectricas inteligentes, Monitoreo y control de restricciones y seguridad, Control de recursos de energ´ ıa renovable, Control basado en datos. Prediction of grid voltages by means of Kriging interpolation Abstract We address the problem of the online prediction of nodal voltage and current trajectories in the distribution grid. For this, we propose a data-driven formulation based on the Kriging interpolation, a machine learning technique that recently showed promising applications in the field of data-based control. We produce a nonparametric prediction oracle which allows to infer future trajectories directly from real-time voltage and current measures. We further provide a simple but effective algorithmic implementation based on the well-known ISTA scheme. We showcase the effectiveness of our methodology for the fast (subsecond) prediction of voltage dynamics through simulations. Keywords: Learning for control, Nonparametric methods, Smart grids, Constraint and security monitoring and control, Control of renewable energy resources, Data-based control. 1. Introducci´ on En un mundo donde las fuentes de energ´ ıa renovable y las tecnolog´ ıas emergentes como los veh´ ıculos el´ ectricos est´ an transformando el panorama energ´ etico, la capacidad de predecir el voltaje y la frecuencia de la red el´ ectrica se ha vuelto crucial. El uso creciente de fuentes de generaci´ on distribuida basadas en energ´ ıas renovables ha introducido desaf´ ıos significativos para mantener la estabilidad de la red, especialmente debido a la naturaleza intermitente de estas fuentes y la carga adicional introducida por los veh´ıculos el´ectricos (Milano et al., 2018). Estos cambios din´amicos pueden causar fluctuaciones en el voltaje y la frecuencia de la red, afectando la calidad del suministro el ´ectrico (Kaur and Vaziri, 2006). Para abordar estos desaf´ıos, la predicci´on en tiempo real del voltaje y la frecuencia de la red el´ectrica es esencial. Una predicci´on precisa no s´olo mejora el monitoreo y el control de los sistemas de energ´ıa, sino que tambi´en juega un papel fundamental en la protecci´on de estos sistemas (Dobbe et al., 2020). Esta capacidad permite optimizar el estado de carga de los sis- ∗Autor para correspondencia: [email protected] Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) Predicci´on de voltajes en la red el´ectrica por interpolaci´on Kriging Carlos Moreno-Blazquez∗, Filiberto Fele, Daniel Limon, Teodoro Alamo Dpto. de Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Av. de los Descubrimientos s/n, 41092, Sevilla, Espa˜na. To cite this article: Moreno-Blazquez, C., Fele, F., Limon, D., Alamo, T. 2025. Prediction of grid voltages by means of Kriging interpolation. Revista Iberoamericana de Automática e Informática Industrial 22, 112-119. https://doi.org/10.4995/riai.2024.21923 temas de almacenamiento de energ´ıa y controlar microrredes y plantas de energ´ıa virtuales de manera ´optima, mejorando la calidad y la estabilidad del suministro el´ectrico (Zufferey et al., 2020). En el ´ambito de la ingenier´ıa de control, los problemas de estimaci´on y predicci´on se han abordado mediante m´etodos param´etricos, que requieren la definici´on d e u n m odelo basado en principios f´ısicos o datos observados. Sin embargo, estos pueden ser limitados por la necesidad de conocer y representar con precisi´on las leyes f´ısicas subyacentes y los par´ametros espec´ıficos del sistema (Car et al., 2021). Recientemente, han surgido m´etodos no param´etricos que ofrecen una alternativa poderosa y flexible, siendo ideales para entornos con alta variabilidad. Estos no requieren una definici´on expl´ıcita de un modelo, pues se basan directamente en los datos observados para realizar predicciones y tomar decisiones de control (Nadales et al., 2023; Ordonez et al., 2021; Merch´anRiveros et al., 2024). Entre estos m´etodos, el Kriging se ha destacado en recientes aplicaciones en el campo de la clasificaci´on y el control (Carnerero et al., 2023). El presente estudio se centra en la aplicaci´on de Kriging para la predicci´on en tiempo real de voltaje y corriente en la red el´ectrica. Este enfoque no solo facilita la gesti´on y control de redes el´ectricas en tiempo real, sino que tambi´en ofrece una robustez superior frente a las fluctuaciones y el ruido en los datos, comparado con los m´etodos param´etricos tradicionales. La r´apida intervenci´on en la red el´ectrica ante eventos inesperados a menudo requiere acciones en escalas de tiempo inferiores a un segundo (Deakin et al., 2021). Por lo tanto, es fundamental contar con la capacidad de prever con precisi´on la evoluci´on del voltaje de la red en intervalos de tiempo muy peque˜nos, del orden de d´ecimas de segundo, como se ha destacado en estudios previos (Gupta and Milanovic, 2006). La predicci´on de trayectorias del voltaje de la red se vuelve fundamental para mantener su equilibrio din´amico, gestionando eficazmente la demanda y la inserci´on de nuevas cargas. Adem ´as, las directivas vigentes, como las limitaciones de tensi´on a corto plazo (ENA Task Group on Statutory Voltage Limits, 2017), enfatizan la importancia de contar con pron´osticos precisos en escalas de tiempo muy peque˜nas para garantizar la estabilidad y fiabilidad del sistema el´ectrico. Para predecir la trayectoria de la red en peque˜nas escalas de tiempo es crucial comprender su comportamiento de forma precisa. Esto requiere datos detallados de la din´amica del sistema, utilizando excitaciones que no alteren su funcionamiento normal. El m´etodo de barrido de frecuencia (Francis et al., 2011; Huang et al., 2009), que inyecta se˜nales sinusoidales de diferentes frecuencias para estimar la impedancia de la red, es com´un en aplicaciones reales. Sin embargo, este m´etodo necesita un dispositivo adicional y es lento por requerir la inyecci´on de m´ultiples se˜nales. Otros m´etodos, que utilizan se˜nales de banda ancha, reducen el tiempo de medici´on. Por ejemplo, los m´etodos de inyecci´on de impulsos (Liu et al., 2020; C´espedes and Sun, 2012), aunque r´apidos, pueden excitar respuestas no lineales y da˜nar la red debido a su alta perturbaci´on. M´etodos avanzados, como el uso de ruido blanco limitado en banda (Xiao et al., 2007), PRBS (Martin et al., 2013; Roinila et al., 2017), se˜nales chirp (Shen et al., 2013) y se˜nales multitonales (Xiao et al., 2022), minimizan la perturbaci´on pero a´un enfrentan problemas, como comentan (Haberle et al., 2023), para estimar de forma param´etrica la impedancia de la red. Estos m´etodos, aunque requieren perturbaciones secuenciales y varios ciclos de medici´on, permiten generar datos no lineales ricos al excitar el sistema de manera controlada en frecuencia y amplitud, sin alterar demasiado la red y capturando su comportamiento no lineal inherente. Este trabajo se enfoca en la propuesta de un m´etodo de predicci´on no param´etrico, Kriging, para las variables de inter´es en la red de distribuci´on. De acuerdo con este enfoque, la operaci´on de la red el´ectrica se perturba con peque˜nas inyecciones de corriente en el punto de acoplamiento com´un (PCC), al fin d e e xcitar s u d in´amica y o btener u n c onjunto d e datos adecuado para la descripci´on del sistema. Los datos de operaci´on recolectados en tiempo real se transforman de forma directa en predicciones de trayectorias futuras, permitiendo su uso en aplicaciones de monitoreo de restricciones y control con requisitos estrictos de tiempo de c´alculo. Como segunda contribuci´on, se describe una posible implementaci´on pr´actica, produciendo una versi´on ad hoc del algoritmo ISTA (del ingl´es Iterative Shrinkage-Thresholding Algorithm), especializada para la estructura del problema Kriging (1). El documento esta´ organizado de la siguiente manera. En la secci´on 2 presentamos el enfoque del problema, y su formulaci´on dual se detalla en las secci´on 3.1. El algoritmo propuesto se presenta en las secci´on 3.2. Posteriormente, se muestra la estructura de la base de datos utilizada, presentando los modelos NARX en la secci´on 4 y destacando la importancia de utilizar el marco dq en la secci´on 4.1. El an´alisis num´erico se presenta en la secci´on 5, antes de las observaciones finales. 2. Formulaci´ on del problema Kriging El Kriging (Cressie, 1990) es una t´ecnica avanzada dentro de la estad´ıstica espacial, que permite realizar predicciones ´optimas en espacios geogr´aficos a partir de datos dispersos. Esta metodolog´ıa, introducida por Matheron en la d´ecada de 1960 (Matheron, 1967), se fundamenta en el concepto de interpolación ponderada, tal como fue refinado y popularizado por Krige en el contexto minero sudafricano (Krige, 1981). Dada su ca-pacidad para modelar la variabilidad espacial con precisi´on ha sido ´util no s´olo en el campo de la miner´ıa, donde las estimacio-nes iniciales basadas en promedios simples eran insuficientes debido a la variabilidad local del dep´osito, sino tambi´en en me-teorolog´ıa, climatolog´ıa y control basados en datos (Carnerero et al., 2023). Seg´un Hemyari and Nofziger (1987), Kriging es una forma de promediado ponderado en la cual los pesos se eligen de manera que el error asociado con el predictor sea menor que para cualquier otra suma lineal. Sea D = {z1, z2, . . . , zN } un conjunto de datos compuesto por las muestras zi ∈ Rn, i = 1, . . . , N. Dado un punto z¯ ∈ Rn, estamos interesados en determinar z¯ en funci´on de una combinaci´on lineal de las zi del conjunto de datos D. Con este objetivo, podemos definir u na funci´on Jγ : Rn × D → [0, ∞) que eval´ua la disimilitud entre z¯ y todos los zi ∈ D. Esta funci´on devolvera´ un valor alto cuando z¯ sea significativamente diferente del conjunto de datos D; por otro lado, valores peque˜nos indicar´an alta similitud. Esto se traslada Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 113 formalmente al siguiente problema (Carnerero et al., 2022): Jγ(¯z,D)Bm´ ın λ1,...,λN N X i=1 wiλ2 i+γi|λi|(1a) s.t. ¯z= N X i=1 ziλi(1b) 1= N X i=1 λi,(1c) donde γi ≥ 0 y wi > 0, i = 1, ..., N, se utilizan para ponderar la informaci´on proporcionada por cada muestra de datos. Esto permite, por ejemplo, filtrar puntos de datos ruidosos o irrelevantes. Los t´erminos γi|λi|i N =1 promueven la esparsidad. Por otro lado, la restricci´on (1b) define z¯ como una combinaci´on lineal de las muestras {zi}N i=1ponderadas por {λi}N i=1; la restricci´on (1c) normaliza la soluci´ on. Ahora, supongamos que est´an disponibles N pares de entrada-salida observados (ui, yi), i = 1, . . . , N, y sea u¯ una nueva entrada, cuya salida y¯ asociada sea desconocida. Definiendo zi = ui (por lo tanto, D = {u1, u2, . . . , uN }), y z¯ = u¯, resolvemos (1) para obtener λ∗ B {λ∗ 1, λ∗ 2, . . . , λ∗ N }. Con esto, podemos obtener una estimaci´on yˆ de y¯ mediante la interpolaci´on Kriging (Carnerero et al., 2022; Matsui and Yamakawa, 2023) ˆy= N X i=1 λ∗ i(¯z,D)·yi.(2) La funci´on (1) es un componente clave en muchos desarrollos novedosos en el contexto de enfoques basados en datos. Formulaciones similares se han utilizado para identificaci´on de sistemas, como en el conocido data-enabled predictive control (v´ease por ejemplo la rese˜na bibliogr´afica de Markovsky et al. (2023)), donde estos m´etodos directos de identificaci´on a partir de trayectorias observadas se fundamentan en el enfoque conductual, o behavioural approach, de Willems (1986). Del mismo modo, enfoques similares como los conocidos m´etodos de optimizaci´on directa de pesos (Roll et al., 2005), permiten encontrar la combinaci´on lineal de puntos en la base de datos que se asemeje lo mejor posible a un punto espec´ıfico. 3. ISTA para Kriging La formulaci´on (1) constituye un programa de optimizaci´on estrictamente convexo sujeto a restricciones convexas; como tal, tiene una soluci´on ´unica λ∗ B (λ∗ 1, λ∗ 2, . . . , λ∗ N ) entre todas las posibles combinaciones lineales de las muestras de datos en D que producen el punto de consulta z¯. El problema (1) se enmarca en la estructura est´andar de un problema cuadr´atico, tras la introducci´on de variables auxiliares para obtener una formulaci´on lineal equivalente del valor absoluto en (1a), resultando en un problema cuya variable de decisi´on tiene dimensi´on 2N. De esta manera la soluci´on de (1) se puede calcular utilizando una amplia variedad de solvers disponibles comercialmente o en acceso abierto. i 3.1. Enfoque dual La adopci´on del enfoque primal en problemas de predicci´on y control esta´ limitada por el hecho de que el n´umero de variables de decisi´on en (1) crece linealmente con el tama˜no N de la base de datos. Observamos que la soluci´on de (1) se puede calcular eficientemente confiando en su reformulaci´on dual (Merhy et al., 2018; Matsui and Yamakawa, 2023). En este caso, el n´umero de variables de decisi´on duales es igual al n´umero de restricciones de igualdad, es decir, n + 1, una cantidad que puede ser significativamente menor que el n´umero de variables primales N. Definiendo r¯ = [z¯⊤1]⊤, ri = [z⊤1]⊤, i = 1, . . . , N, podemos reescribir (1) como J∗ γ=m´ ın λ1,...,λN N X i=1 wiλ2 i+γi|λi|(3a) s.t. ¯r= N X i=1 riλi,(3b) y consideramos la siguiente suposici´on a lo largo del art´ıculo Suposici´on 1. El problema de optimizaci´on (3) es factible, es decir, la matriz R B [r1, r2, . . . , rN ] tiene rango de fila completo. Ahora, sea µ ∈ Rn+1 la variable dual asociada a las restricciones de igualdad (3b). Esto conduce a la funci´on dual ϕ(µ)=−µ⊤¯r+m´ ın λ1,...,λN N X i=1 wiλ2 i+γi|λi|+(µ⊤ri)λi.(4) Fijado µ, el problema (4) es separable, y los minimizadores pueden calcularse para i = 1, . . . , N como λ∗ i(µ)=argm´ ın λi∈Rwiλ2 i+γi|λi|+(µ⊤ri)λi.(5) Este problema de optimizaci´on escalar posee una soluci´on en forma cerrada (Beck, 2017, Example 6.8), como se afirma en la siguiente proposici´on. Proposici´on 2. Consideremos el problema de optimizaci´on z∗=argm´ ın z∈Rwz2+γ|z|+cz,(6) con w >0,γ≥0, y c ∈R. Definamos ψ:R→Rcomo ψ(w, γ, c)B         sign(c) γ− |c| 2w!si |c|> γ, 0en otro caso. (7) Entonces, z∗ = ψ(w, γ, c). 3.2. Algoritmo K-ISTA Bajo la Suposici´on 1, el problema (1) cumple las condiciones para la dualidad fuerte (Bertsekas, 2009, Prop. 5.3.3), lo que significa q ue Jγ ∗ = m´axµ ϕ(µ). A dem´as, definiendo W B diag([w1, w2, . . . , wN ]), es posible demostrar, v´ease por ejemplo (Beck, 2017), que ϕ(µ) satisface ϕ(µ+ ∆µ)≥ϕ(µ)+ ∆µ⊤g(µ)−1 4∆µ⊤RW−1R⊤∆µ, (8) para cada µ, ∆µ∈Rn+1, donde g(µ)B−¯r+ N X i=1 riλ∗ i(µ).(9) 114 Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 Esta propiedad proporciona una forma intuitiva de maximizar la funci´on dual ϕ(µ). Partiendo de cualquier µ fijo, podemos maximizar el lado derecho de (8) con respecto a ∆µ a lo largo de iteraciones sucesivas; esto nos permitira´ recuperar el maximizador µ∗ de (4) a medida que ∆µ∗ → 0. Este ´ultimo puede obtenerse expl´ıcitamente al diferenciar el lado derecho en (8) y resolver para ∆µ, como ∆µ∗= Ω−1g(µ),(10) con ΩB1 2RW−1R⊤. Teniendo en cuenta que W∈RN×N, con conjuntos de datos grandes, una forma eficiente en t´ erminos de memoria para calcularlo es Ω = PN i=1 1 2wrir⊤ i. La implementaci´on algor´ıtmica propuesta para la soluci´on de (1) se basa en el conocido algoritmo de contracci´on y umbral iterativo (ISTA) (v´ease (Beck and Teboulle, 2009; Alamo et al., 2019) y sus referencias) que especializamos a la reformulaci´on dual particular de (1) adoptada en este art´ıculo, utilizando el resultado de la Proposici´on 2. Algorithm 1 K-ISTA i: Entradas ϵ > 0. ii: Inicializaci´ on k=0, µ0=0n+1. iii: repeat iv: k=k+1. v: for i=1, . . . , Ndo vi: λk i=ψ(wi, γi,r⊤ iµk−1). vii: end for viii: gk=−¯r+PN i=1riλk i. ix: µk=µk−1+ Ω−1gk. x: until ∥gk∥< ϵ xi: Salidas λ∗← {λk i}N i=1 1 k k2 En el Algoritmo 1, primero i) se define la tolerancia ϵ > 0 para establecer el criterio de parada con el cuasi cumplimiento de las restricciones de igualdad. ii) Se inicializan el contador de iteraciones k = 0 y la variable dual µ0 = 0n+1. iii) El algoritmo itera hasta que la norma del gradiente gk sea menor que la tolerancia impuesta ϵ. iv) En cada iteraci´on, se incrementa el contador k en uno. v) Para cada i : i = 1, . . . , N, en vi) se actualiza λi k usando la ecuaci´on (7). viii) Se recalcula el gradiente gk como en (9), actualizando seguidamente en ix) µk con (10). Finalmente, en xi) el conjunto de coeficientes que satisface la condici´on de parada forma el vector soluci´on λ∗. Este algoritmo hereda la simplicidad del m´etodo ISTA, demostrando una tasa de convergencia de O( ). Aunque el algoritmo propuesto llega a ser compatible con computaciones en tiempo real, es posible obtener variantes optimizadas que disminuyen el n´umero de iteraciones y, por ende, el tiempo de c´omputo. Esto se demuestra en la secci´on 5, donde los resultados han sido obtenidos utilizando K-ISTA y una versi´on del mismo basada en el esquema FISTA con restart, que acelera la tasa de convergencia a O( 1 ) (v´ease (Alamo et al., 2019) para m´ as detalle). 4. Composici´ on del conjunto de datos NARX El sistema a controlar se caracteriza por sus entradas, u∈ U ⊂ Rm, y sus salidas y∈ Y ⊂ Rp, donde UeYse suponen convexos. La ´ unica informaci´ on disponible de este sistema es un conjunto hist´orico adecuado de N pares observados de entrada-salida {(u1, y˜1), . . . , (uN , y˜N )}, donde y˜ denota la medida posiblemente ruidosa de la salida, con el ruido asumido acotado por ε. Asumimos que la din´amica del sistema puede describirse mediante un modelo autoregresivo no lineal con entradas ex´ogenas (NARX) (Chen and Billings, 1989; Leontaritis and Billings, 1985) como y(k+1), . . . , y(k+Np)(11) =fy(k), . . . , y(k−na),u(k), . . . , u(k−nb), donde na, nb ∈ N parametrizan el horizonte para los valores pasados de salida y entrada, respectivamente. Bajo suposiciones leves sobre la observabilidad del sistema, el modelo NARX puede describir la din´amica de un sistema no lineal para horizontes suficientemente grandes. Sin embargo, pueden existir errores de modelado debido a una descripci´on incompleta de la din´amica real del sistema representada por f . Sea n¯ el orden del sistema. Como demostrado por Levin and Narendra (1997), se puede lograr una representaci´on perfecta del sistema seleccionando los horizontes de memoria na y nb como na = nb = 2¯n. En caso de que el orden del sistema sea desconocido, los valores de na y nb deben estimarse mediante alg´un m´etodo de validaci´on cruzada. Nos referimos a los par´ametros de entrada de f en (11) como el regresor zk B (x(k), u(k)), donde x(k) B (y(k), y(k − 1), ..., y(k −na), u(k − 1), ..., u(k − nb)) es el estado en el instante de muestra k. Adem´as, sea z˜(k) el regresor ruidoso que inclu-ye las observaciones corruptas, es decir, z˜k B (x˜(k), u(k)) = (y˜(k), y˜(k − 1), ..., y˜(k − na), u(k), u(k − 1), ..., u(k − nb)), y sea D el conjunto de regresores conformado por datos entrada-salida reorganizados de la siguiente manera: D=˜zi,i=1, ..., N.(12) Asimismo, las trayectorias futuras a Nppasos de la salida, correspondientes a cada uno de los regresores ˜zi, se almacenan en otro conjunto Ξ = ξi,i=1, ..., N, donde ξi= ˜y(i+1), . . . , ˜y(i+Np). 4.1. Marco dq Una vez clara la estructura del modelo NARX, se va a emplear el marco dq para representar las medidas trif´asicas. Este es un marco de referencia que considera dos ejes, llamados eje directo y eje en cuadratura, que giran a una frecuencia dada (por ejemplo, la frecuencia de red). Esto permite convertir las componentes de corriente y voltaje en constantes en estado estacionario, lo que simplifica significativamente el dise˜no y an´alisis de controladores. En el marco dq, las componentes directas (d) y en cuadratura (q) pueden ser controladas de manera independiente. Esta capacidad de desacoplamiento es particularmente ventajosa, por ejemplo, para la implementaci´on de controladores proporcionales-integrales (PI) (Francis et al., 2011). Adem´as, las componentes de alta frecuencia, o arm´onicos, son m´as f´aciles de filtrar en este marco, lo que mejora la calidad de la se˜nal y la eficiencia del sistema (Saccomando and Svensson, 2001). La facilidad de las transformaciones directas e inversas entre los marcos abc y dq tambi´en contribuye a su preferencia, ya Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 115 Figura 1: Modelo Matlab Simulink de un prosumer convencional conect´andose a la red en el PCC. Se toman datos de corriente y voltaje en el marco dq con el objetivo de conformar el dataset D. As´ı mismo, se recogen datos del ´angulo de fase ωgt para poder utilizar las predicciones dq en el marco abc. que estas transformaciones est´an bien definidas y pueden implementarse f´acilmente en sistemas de control digital. Por tanto, haciendo uso de las ventajas del marco dq y conformando el conjunto de datos como en (12), nuestro objetivo es inferir una aproximaci´on F basada en aprendizaje (impl´ıcito) del modelo NARX en (11) que, dado un regresor del sistema z¯ = (x¯, u¯), devuelva un sucesor ξˆ ∈ RpNp . En nuestro caso, ˆ ξ=F(¯z,D)= N X i=1 λ∗ i(¯z,D)·ξi,(13) donde F toma la forma de (2). Esta expresi´on hace uso de los pesos {λ∗ 1, . . . , λ∗ N } obtenidos al resolver (1), en funci´on de z¯ y D, mediante el Algoritmo 1. 5. An´ alisis num´ erico Para probar el enfoque propuesto, consideramos un caso de simulaci´on implementado en Matlab Simulink, como puede verse en la Figura 1. Lo que se pretende es estimar la trayectoria de voltaje de la red en el marco dq a partir de medidas de corriente y voltaje obtenidas de una serie de sensores situados en el PCC. Cabe comentar que los resultados se muestran en el marco abc, lo cual es posible realizando predicciones al mismo tiempo de la frecuencia de la red ωg haciendo uso del modelo utilizado para predecir el voltaje. Tabla 1: Par´ametros El´ectricos del Experimento Num´erico Par´ ametro S´ ımbolo Valor Valores base Vb,Sb,fb380 V, 1.5 kVA, 50 Hz Carga R12 p.u. L´ ınea 1 R2,L2,C20.015 p.u., 0.15 p.u., 0.05 p.u. L´ ınea 2 R3,L3,C30.015 p.u., 0.15 p.u., 10 p.u. Frecuencia de muestreo fs1 kHz Horizonte de predicci´ on Np200 Los par´ametros el´ectricos y de medici´on del sistema representado en la Figura 1 pueden verse contenidos en la Tabla 1. Dado que el fin ´ ultimo e s e stimar t rayectorias e n e scalas de tiempo inferiores al segundo, se ha elegido un tiempo de c´ alculo de la trayectoria de 0,2 s. Dada la frecuencia de muestreo fs, se deduce el valor del horizonte de predicci´ on Np=0,2fs. Dado que el voltaje en dq tiene dos componentes, en este caso la ˆ ξ∈R2pNp, siendo en este caso 2Np=400 valores de dimensi´ on pque deben ser predichos. incluido; en este caso, el rango fue de 2,5,1,8 Figura 2: Voltajes y corrientes trif´asicos perturbados para la obtenci´on de la base de datos NARX utilizada en la regresi´on. Durante el experimento de identificaci´on de la red, consideramos condiciones de red constantes y estacionarias. Para excitar la red el´ectrica bajo prueba y poder generar D, el sistema se excito´ en corriente en el marco dq mediante una combinaci´on de se˜nales chirp (oscilaciones con frecuencia y amplitud variable a lo largo del tiempo) y PRBS (Nelles, 2020) no correlacionadas, cada una con amplitud de 0,1 p.u. como m´aximo. El rango de frecuencias de excitaci´on para la se˜nal chirp y PRBS se ha escogido de modo que queden por debajo de la frecuencia de muestreo, y as´ı, el pico de resonancia  fs del fs  sistema estuviera rad/s. De este experimento se obtuvieron 9201 pares de entrada-salida. De estos, se extrajo la base de datos de entrenamiento D. Para validar el modelo, utilizamos un conjunto de datos de 4501 elementos, obtenidos mediante se˜ nales chirp y PRBS del mismo modo. 116 Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 Figura 3: Comparaci´on entre la trayectoria generada en simulaci´on y la estimada por (1)-(2) haciendo uso del Algoritmo 1, K-ISTA, en dos esce-narios diferentes utilizando la base de datos NARX con na = nb = 4. La gr´afica de arriba representa un escenario donde las tres componentes son excitadas al mismo tiempo, mientras que abajo ´unicamente las fases a y c son excitadas, dejando a b en todo momento sin excitar. Las respuestas de corriente y voltaje resultantes en el PCC se muestran en la Figura 2. Se observa como el nivel de perturbaci´on es bastante peque˜no en comparaci´on con el caso no excitado. Dado que la excitaci´on tambi´en esta ´ limitada a unos pocos segundos, no deteriora la operaci´on continua de la red. Analizamos este aspecto m´as en detalle, para asegurar que el funcionamiento normal de la red no se vea afectado y que se cumplen los l´ımites establecidos por el est´andar IEEE (IEEE Power and Energy Society, 2021) y el est´andar Europeo (CENELEC (European Committee for Electrotechnical Standardisation), 2011). ´Estos fijan un l´ımite del 8 % de distorsi´on arm´onica total (THD) en voltaje, sobre 10 ciclos de la se˜nal. Nuestro an´alisis sobre diversos segmentos de 10 ciclos ha verificado que el THD en la se˜nal excitada no excede el 6,74 %, asegurando as´ı el cumplimiento de las normativas vigentes y la integridad operativa de la red. Esta evaluaci´on se ilustra gr´aficamente en la Figura 4, donde se presenta el espectro de frecuencias y el THD calculado utilizando el segmento de la se˜nal de voltaje asociado con la m´axima distorsi´on. D esta´ en la forma de (12), donde se ha impuesto que na = nb y se han tenido en cuenta diferentes valores, como puede verse en la Tabla 2. Cabe destacar que cada regresor z˜ tiene dimensi´on n = (na + 1)p + (nb + 1)m. Para mejorar el tiempo de c´alculo y la calidad de la predicci´on, se calcula una partici´on de D llamada D◦, la cual sera´ la base de datos que se utilizara´ para realizar la predicci´on. Esta partici´on se conforma exclusivamente para cada punto de consulta z¯, tomando los N◦ = 1000 puntos de la base de datos original D que est´an m´as cercanos, en t´erminos de distancia eucl´ ıdea, a ¯z; la cardinalidad N◦se eligi´ o basado en distintas pruebas de simulaci´ on, donde se evalu´ o la precisi´ on de la predicci´ on. El peso wiasociado a cada uno de los puntos en D◦se ha ajustado para asignar mayor importancia a los puntos m´ as cercanos a ¯r(es decir, penalizando menos los coeficientes λi resultantes para estos puntos). Figura 4: An´alisis de la Distorsi´on Arm´onica Total (THD) en el voltaje de la fase c del sistema trif´asico abc. En la gr´afica superior se muestra la forma de onda temporal del voltaje, donde se marcan en rojo los 10 ciclos seleccionados para obtener el espectro de frecuencia con el cual se calcula el THD representado en la gr´afica inferior, correspondiente al peor de los casos. Este an´alisis asegura que la se˜nal cumple con el limite máximo permitido según los datos experimentales utilizados. Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 117 En este caso, se ha optado por definir wi B exp(α∥z¯ − z˜i∥s), donde α > 0 y s ≥ 1 son par´ametros ajustables. La modificaci´on de estos par´ametros permite controlar el peso de los puntos utilizados en la regresi´on, garantizando que el problema (1) este´ siempre bien condicionado num´ericamente. Los datos experimentales se han obtenido con s = 1 y haciendo uso de normas eucl´ıdeas, tomando α = log (1, 05) /r◦, donde r◦ se define como la distancia m´axima entre z¯ y cualquier punto en D◦, es decir, r◦ = m´axi∈{1,...,N◦}∥z˜i − z¯∥2. Para validar el predictor, calculamos el error relativo normalizado como ζ = ∥ξv − ξˆ∥∞/∥ξv∥∞, donde ξv es una muestra de salida del conjunto de datos de validaci´on, y ξˆ la salida estimada por Kriging. En la Tabla 2 es posible ver una comparaci´on de los tiempos de c´alculo de la estimaci´on con horizonte Np = 200, teniendo en cuenta que se ha utilizado un ordenador equipado con un procesador Intel Core i7 de 8ª generaci´on y 16 GB de RAM. La fila que se destaca en dicha tabla hace referencia a la mejor elecci´on de par´ametros na = nb = 4, escogida por produ-cir el menor error medio de predicci´on obtenido con los datos de validaci´on. Destacar que los tiempos mostrados en la Tabla 2 se han obtenido tanto con el algoritmo K-ISTA presentado en la secci´on 3.2, como con una versi´on acelerada de dicho algoritmo basada en el esquema FISTA con restart (Alamo et al., 2019), imponiendo una tolerancia ϵ = 10 −3 para establecer el criterio de parada. Adicionalmente, en Figura 3 se muestra una comparativa de la calidad de la predicci´on realizada tras resolver (1)-(2) haciendo uso del Algoritmo 1, K-ISTA, en dos escenarios diferentes. eficientes {λ∗ 1, . . . , λ∗ N◦}que satisfacen la relaci´ on ¯z=i=1 Figura 5: Comparaci´on de dos trayectorias predichas para la fase c: modelo ARX (param´etrico) frente a predicci´on basada en Kriging (no param´etrico). En la gr´afica superior, se evidencia que el modelo ARX, entrenado con la totalidad de la base de datos, no estima adecuadamen-te la se˜nal no excitada. Por otro lado, en la gr´afica inferior, el modelo basado en Kriging demuestra una mayor precisi´on en la captura de las excitaciones de la se˜nal, superando al modelo ARX en la detecci´on de dichas perturbaciones. Como se discutio´ en la secci´on 2, esto se realiza en dos pasos: primero, se resuelve (1) para obtener el conjunto PN◦de λ∗ coi˜zi; ˆ ξ= luego, P la trayectoria de voltaje ξˆ se estima a trav´es de (2), como N◦ i=1λ∗ iξi. Del mismo modo, se predice la trayectoria de la frecuencia de la red ˆωgusando los mismos coeficientes λipara poder pasar la predicci´ on de voltaje del marco dq al marco abc. De esta forma, ˆωg=PN◦ i=1λ∗ i˜ωg,i. Asimismo, para evaluar la efectividad del predictor, se ha comparado con un modelo ARX ajustado sobre D por m´ınimos cuadrados ordinarios, utilizando como m´etrica de evaluaci´on la media de ζ sobre todo el conjunto de validaci´on. Con na, nb = 4, el predictor basado en K-ISTA alcanza una precisi´on del ζ = 3,78 %, comparado con el ζ = 5,04 % del modelo ARX, pudiendo estimar de forma m´as exacta las excitaciones a las que esta´ siendo sometido el modelo simulado. Esto puede verse en la Figura 5, donde se ha comparado la predicci´on de trayectorias para la fase c, mostrando que el modelo basado en Kriging tiene una mayor precisi´on en la captura de las excitaciones de la se˜nal en comparaci´on con el modelo ARX. Tabla 2: Tiempo de c´alculo (CT) y error relativo (ζ) para valores distintos de na y nb, tomando ϵ = 10−3. na,nbCT K-ISTA [ms] CT K-ISTA acelerado [ms] ζ[%] Max Median Min Max Median Min Mean 21135 25.80 0.47 28.40 5.97 0.26 3.92 4715 18.80 0.58 41.24 5.04 0.28 3.78 8919 19.39 0.74 64.41 5.52 0.50 3.95 16 3372 14.17 0.71 49.88 5.08 0.40 5.61 25 3686 12.96 0.35 74.79 4.81 0.37 11.3 50 4896 19.80 1.39 91.26 7.89 0.70 6.29 6. Conclusiones En este trabajo se ha propuesto un m´ etodo de predicci´ on no param´ etrico basado en Kriging para las variables de inter´ es en la red de distribuci´ on. Los resultados demuestran que la precisi´ on de la predicci´ on puede llegar a ser inferior al 5 %, con un tiempo de c´ alculo de menos de 200 ms en el peor de los casos, para grandes horizontes de predicci´ on. Esto destaca su aplicabilidad en aplicaciones de la red que requieran maniobras en escalas de tiempo inferiores al segundo. Agradecimientos Expresar nuestro agradecimiento a Mar´ ıa Camila Merch´ an Riveros por su valiosa ayuda en el modelado de redes el´ ectricas para nuestro art´ ıculo. Este trabajo ha sido realizado en el marco del Proyecto de investigaci´ on PID2022-142946NA-I00 financiado por MICIU/AEI /10.13039/501100011033 y por FEDER, UE, y del Proyecto 2023/00000487 financiado por el VIIPP-2022 de la Universidad de Sevilla. F. Fele y C. Moreno Bl´ azquez tambi´ en agradecen el soporte de la ayuda RYC2021-033960-I financiada por MICIU/AEI /10.13039/501100011033 y por la Uni´ on Europea NextGenerationEU/PRTR. Referencias Alamo, T., Krupa, P., Limon, D., 2019. Gradient based restart FISTA. In: 2019 IEEE 58th Conference on Decision and Control (CDC). pp. 3936–3941. Beck, A., 2017. First-order methods in optimization. SIAM. Beck, A., Teboulle, M., 2009. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences 2 (1), 183– 202. Bertsekas, D. P., 2009. Convex optimization theory. Athena Scientific. 118 Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 Car, M., Leˇ si´ c, V., Vaˇ sak, M., 2021. Cascaded control of back-to-back converter dc link voltage robust to grid parameters variation. IEEE Transactions on Industrial Electronics 68 (3), 1994–2004. Carnerero, A. D., Ramirez, D. R., Alamo, T., 2022. Probabilistic interval predictor based on dissimilarity functions. IEEE Transactions on Automatic Control 67 (12), 6842–6849. Carnerero, A. D., Ramirez, D. R., Limon, D., Alamo, T., 2023. Kernel-based state-space Kriging for predictive control. IEEE/CAA Journal of Automatica Sinica 10 (5), 1263–1275. CENELEC (European Committee for Electrotechnical Standardisation), 2011. Voltage characteristics of electricity supplied by public distribution systems. European Norm EN 50160. C´ espedes, M., Sun, J., 2012. Online grid impedance identification for adaptive control of grid-connected inverters. In: 2012 IEEE Energy Conversion Congress and Exposition (ECCE). IEEE, pp. 914–921. Chen, S., Billings, S. A., 1989. Representations of non-linear systems: the NARMAX model. International journal of control 49 (3), 1013–1032. Cressie, N., 1990. The origins of Kriging. Mathematical geology 22, 239–252. Deakin, M., Greenwood, D. M., Taylor, P. C., Armstrong, P., Walker, S., 2021. Analysis of network impacts of frequency containment provided by domestic-scale devices using matrix factorization. IEEE Transactions on Power Systems 36 (6), 5697–5707. Dobbe, R., van Westering, W., Liu, S., Arnold, D., Callaway, D., Tomlin, C., 2020. Linear singleand three-phase voltage forecasting and Bayesian state estimation with limited sensing. IEEE Transactions on Power Systems 35 (3), 1674–1683. ENA Task Group on Statutory Voltage Limits, 2017. Statutory voltage limits at customers’ terminals in the UK and options for future application of wider limits at low voltage. Tech. Rep. ETR 140, Energy Networks Association. Francis, G., Burgos, R., Boroyevich, D., Wang, F., Karimi, K., 2011. An algorithm and implementation system for measuring impedance in the dq domain. In: 2011 IEEE Energy Conversion Congress and Exposition. IEEE, pp. 3221–3228. Gupta, C. P., Milanovic, J. V., 2006. Probabilistic assessment of equipment trips due to voltage sags. IEEE Transactions on power delivery 21 (2), 711–718. Haberle, V., Huang, L., He, X., Prieto-Araujo, E., Smith, R. S., Dorfler, F., 2023. MIMO grid impedance identification of three-phase power systems: Parametric vs. nonparametric approaches. In: 2023 62nd IEEE Conference on Decision and Control (CDC). IEEE, pp. 542–548. Hemyari, P., Nofziger, D., 1987. Analytical solution for punctual Kriging in one dimension. Soil Science Society of America journal 51 (1), 268–269. Huang, J., Corzine, K. A., Belkhayat, M., 2009. Small-signal impedance measurement of power-electronics-based ac power systems using line-to-line current injection. IEEE Transactions on Power Electronics 24 (2), 445–455. IEEE Power and Energy Society, 2021. IEEE draft standard for harmonic control in electric power systems. IEEE P519/D5.1, January 2021, 1–30. Kaur, G., Vaziri, M., 2006. Effects of distributed generation (DG) interconnections on protection of distribution feeders. In: 2006 IEEE Power Engineering Society General Meeting. Krige, D. G., 1981. Lognormal-de Wijsian geostatistics for ore evaluation. South African Institute of mining and metallurgy Johannesburg. Leontaritis, I. J., Billings, S. A., 1985. Input-output parametric models for nonlinear systems part i: deterministic non-linear systems. International journal of control 41 (2), 303–328. Levin, A., Narendra, K., 1997. Identification of nonlinear dynamical systems using neural networks. In: Neural Systems for Control. Elsevier, pp. 129– 160. Liu, Z., Liu, J., Liu, Z., 2020. Analysis, design, and implementation of impulseinjection-based online grid impedance identification with grid-tied converters. IEEE Transactions on Power Electronics 35 (12), 12959–12976. Markovsky, I., Huang, L., D¨ orfler, F., 2023. Data-driven control based on the behavioral approach: From theory to applications in power systems. IEEE Control Systems Magazine 43 (5), 28–68. Martin, D., Nam, I., Siegers, J., Santi, E., 2013. Wide bandwidth three-phase impedance identification using existing power electronics inverter. In: 2013 Twenty-Eighth Annual IEEE Applied Power Electronics Conference and Exposition (APEC). IEEE, pp. 334–341. Matheron, G., 1967. Kriging or polynomial interpolation procedures. CIMM Transactions 70 (1), 240–244. Matsui, H., Yamakawa, Y., 2023. Sparse estimation in ordinary Kriging for functional data. arXiv preprint arXiv:2306.15537. Merch´ an-Riveros, M. C., Albea, C., Seuret, A., 2024. Data-driven control design for power converters approximated as switched affine systems and experimental validation. IEEE Transactions on Circuits and Systems II: Express Briefs, 1–1. Merhy, D., Alamo, T., Stoica Maniu, C., Camacho, E. F., 2018. Zonotopic constrained kalman filter based on a dual formulation. In: 2018 IEEE Conference on Decision and Control (CDC). pp. 6396–6401. Milano, F., D¨ orfler, F., Hug, G., Hill, D. J., Verbiˇ c, G., 2018. Foundations and challenges of low-inertia systems (invited paper). In: 2018 Power Systems Computation Conference (PSCC). pp. 1–25. Nadales, J., Carnerero, A., Moreno-Blazquez, C., Haes-Ellis, R., Limon, D., 2023. Learning-based NMPC on SoC platforms for real-time applications using parallel lipschitz interpolation. IFAC-PapersOnLine 56 (2), 6298– 6303, 22nd IFAC World Congress. Nelles, O., 2020. Nonlinear dynamic system identification. Springer. Ordonez, J. G., Nadales, J. M., Limon, D., Gordillo, F., 2021. Data-driven multirate predictive control of power inverters based on kinky inference. In: 2021 60th IEEE Conference on Decision and Control (CDC). pp. 4358– 4363. Roinila, T., Messo, T., Santi, E., 2017. MIMO-identification techniques for rapid impedance-based stability assessment of three-phase systems in dq domain. IEEE Transactions on Power Electronics 33 (5), 4015–4022. Roll, J., Nazin, A., Ljung, L., 2005. Nonlinear system identification via direct weight optimization. Automatica 41 (3), 475–490. Saccomando, G., Svensson, J., 2001. Transient operation of grid-connected voltage source converter under unbalanced voltage conditions. In: Conference Record of the 2001 IEEE Industry Applications Conference. 36th IAS Annual Meeting (Cat. No. 01CH37248). Vol. 4. IEEE, pp. 2419–2424. Shen, Z., Jaksic, M., Mattavelli, P., Boroyevich, D., Verhulst, J., Belkhayat, M., 2013. Three-phase ac system impedance measurement unit (IMU) using chirp signal injection. In: 2013 Twenty-Eighth Annual IEEE Applied Power Electronics Conference and Exposition (APEC). IEEE, pp. 2666–2673. Willems, J. C., 1986. From time series to linear system—part i. finite dimensional linear time invariant systems. Automatica 22 (5), 561–580. Xiao, D., Hu, H., Chen, S., Song, Y., Pan, P., Molinas, M., 2022. Rapid dqframe impedance measurement for three-phase grid based on interphase current injection and fictitious disturbance excitation. IEEE Transactions on Instrumentation and Measurement 71, 1–13. Xiao, P., Venayagamoorthy, G., Corzine, K., 2007. A novel impedance measurement technique for power electronic systems. In: 2007 IEEE Power Electronics Specialists Conference. IEEE, pp. 955–960. Zufferey, T., Renggli, S., Hug, G., 2020. Probabilistic state forecasting and optimal voltage control in distribution grids under uncertainty. Electric Power Systems Research 188, 106562. Moreno-Blazquez, C. et al. / Revista Iberoamericana de Automática e Informática Industrial 22 (2025) 112-119 119