Optimización de funciones de coste para control predictivo de máquinas de inducción multifásicas

www.revista-riai.org Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 Manuel R. Arahala, Agnieszka Kowal G.b,∗, Federico Barrerob, Mar´ ıa del Mar Castillaa aDpto. de Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Camino de los descubrimientos s/n, 41092, Sevilla, Espa˜na. bDpto. de Tecnolog´ıa Electr´onica, Universidad de Sevilla, Camino de los descubrimientos s/n, 41092, Sevilla, Espa˜na. Resumen El control de m´ aquinas de inducci´ on multif´ asicas avanza con la incorporaci´ on de controladores predictivos que act´ uan directamente sobre el inversor sin mediar moduladores. Se lleva a cabo en tiempo discreto, mediante una optimizaci´ on exhaustiva de los estados posibles del inversor. La optimizaci´ on utiliza una funci´ on de coste que penaliza desviaciones de la referencia y otros efectos indeseados, por lo que su dise˜ no tiene gran importancia en el comportamiento en bucle cerrado. En este art´ ıculo se presenta como contribuci´ on una nueva estrategia para el dise˜ no de funciones de coste y su optimizaci´ on. La estrategia define un ´ ındice ponderado de comportamiento basado en el funcionamiento del sistema. Las funciones de coste se eligen de manera que se optimice dicho ´ ındice. Se presentan resultados de simulaci´ on correspondientes a una m´ aquina de cinco fases y se compara con el m´ etodo tradicional, todo ello en el contexto de control predictivo de corrientes de est´ ator. Palabras Clave: Control predictivo, M´ aquinas de inducci´ on, Inversores, M´ aquinas polif´ asicas, Optimizaci´ on global Cost Function Optimization for Multi-phase Induction Machines Predictive Control Abstract Control of multi-phase induction machines is evolving with the incorporation of predictive controllers that operate directly with the power converter without modulators. It is realized in discrete-time through an exhaustive optimization of the inverter states. The optimization uses a cost function that penalizes deviations from the reference and other undesired factors. The cost function design has a big impact on closed loop performance. In this paper, a new strategy for cost function design is presented as a contribution. The strategy defines and a performance index relative to the observed behavior. The cost functions are derived from an optimization of said index. Experimental and simulation results are included for a five-phase induction machine with a comparison with previous methods in the context of predictive current control. Keywords: Predictive control, Induction machines, Inverters, Multi-phase machines, Global optimization 1. Introducci´ on Las m´ aquinas de inducci´ on multif´ asicas han recibido recientemente un fuerte impulso frente a las convencionales trif´ asicas debido por un lado a sus cualidades inherentes (robustez, disminuci´ on de arm´ onicos, mayor densidad de potencia, reconfigurabilidad (Zoric et al., 2018)) y por otro a la aparici´ on de nuevos m´ etodos y equipos de control que pueden lidiar con la complejidad a˜ nadida por las fases extra. En este ´ ultimo cap´ ıtulo juega un papel importante el desarrollo de controladores digitales que act´ uan directamente sobre el convertidor de potencia, sin usar moduladores intermedios. La eliminaci´ on del modulador permite un control potencialmente mucho m´ as r´ apido (Lim et al., 2014), a cambio se necesita alguna manera de ligar la respuesta deseada a los estados del convertidor. El control predictivo basado en modelo (Camacho y Bordons, 2013) es una herramienta muy adecuada para esta tarea y ha sido empleado en la forma descrita recientemente. Con anterioridad se hab´ ıa usado control predictivo en tiempo continuo ∗Autor para correspondencia: [email protected] To cite this article: Manuel R. Arahal, Agnieszka Kowal G., Federico Barrero, María del Mar Castilla. 2019. Cost Function Optimization for Multi-phase Induction Machines Predictive Control. Revista Iberoamericana de Automática e Informática Industrial 16, 48-55. https://doi.org/10.4995/riai.2018.9771 Attribution-NonCommercial-NoDerivatives 4,0 International (CC BY-NC-ND 4,0) Optimización de Funciones de Coste para Control Predictivo de Máquinas de Inducción Multifásicas usando un modulador de ancho de pulso para convertir la actuaci´ on calculada en estados del inversor (Holtz y Stadtfeld, 1983). Con el desarrollo de procesadores de se˜ nal m´ as r´ apidos se pudo atacar el problema de resolver en tiempo real la optimizaci´ on necesaria (Holmes y Martin, 1996), (Abu-Rub et al., 2001). En el contexto del control de m´ aquinas de inducci´ on polif´ asicas, el control predictivo sin usar modulaci´ on ha recibido la denominaci´ on de control predictivo de estado finito. Los autores Arahal y Barrero aplicaron esta estrategia en (Arahal et al., 2009) a una m´ aquina de seis fases. Para ello usaron un algoritmo de tiempo discreto que selecciona mediante una optimizaci´ on exhaustiva el estado del inversor para su uso en el siguiente periodo de muestreo. La optimizaci´ on utiliza una funci´ on de coste que penaliza desviaciones respecto de las corrientes de referencia y que puede incluir otros efectos indeseados como distorsi´ on arm´ onica y p´ erdidas por conmutaci´ on. El dise˜ no de la funci´ on de coste tiene gran importancia en el comportamiento del sistema en bucle cerrado. En el ya citado art´ ıculo de Arahal y Barrero (Arahal et al., 2009), se abord´ o por primera vez el problema del dise˜ no de la funci´ on de coste mostr´ andose la flexibilidad que proporciona para acomodar los objetivos de seguimiento en ejes α−β(relacionados con la generaci´ on de par) y en ejes x−y(relacionados con p´ erdidas). Continuando con esta l´ ınea, en (Barrero et al., 2009) se present´ o una funci´ on de coste que incluye penalizaci´ on del n´ umero de conmutaciones del inversor, lo que permite incluir las p´ erdidas por conmutaci´ on como criterio en el controlador. A partir de ah´ ı se han propuesto diversos m´ etodos para el dise˜ no de funciones de coste. El lector interesado puede consultar el reciente art´ ıculo de revisi´ on (Mamdouh et al., 2018) y las citas que incluye. En particular (Rojas et al., 2013) donde se presenta un m´ etodo que evita la sinton´ ıa de par´ ametros de la funci´ on de coste para el caso trif´ asico de control de par. En contraposici´ on, la presente propuesta utiliza expl´ ıcitamente dichos par´ ametros y propone una nueva estrategia para el dise˜ no de funciones de coste para control predictivo de corriente de m´ aquinas polif´ asicas. La estrategia define un ´ ındice de comportamiento basado en el funcionamiento observado del sistema. Las funciones de coste se eligen de forma que optimicen dicho ´ ındice. La novedad reside en considerar los coeficientes de la funci´ on de coste como variables que dependen del punto de operaci´ on. Este cambio conlleva la necesidad de considerar el comportamiento de la m´ aquina de inducci´ on en intervalos de funcionamiento. Se pasa as´ ı, como se ver´ a, de una situaci´ on en que los coeficientes de la funci´ on de coste afectan de manera poco clara al comportamiento global a otra en la que se puede intuir su efecto. De este modo se consiguen dos cosas: por un lado una relaci´ on m´ as intuitiva entre par´ ametros y comportamiento y por otra un comportamiento optimizado para cada punto de funcionamiento. La t´ ecnica propuesta se presenta a continuaci´ on en el contexto general de m´ aquinas polif´ asicas. Posteriormente, en la Secci´ on 3 se aplica a una m´ aquina de inducci´ on de cinco fases. M´ as adelante, en la Secci´ on 4, se presentan los resultados obtenidos, en los que se compara la nueva estrategia con el m´ etodo tradicional de dise˜ no de funciones de coste. 2. Nueva Estrategia para el Dise˜ no de Funciones de Coste A fin de presentar con rigor la t´ ecnica propuesta se usar´ an una gran cantidad de s´ ımbolos que se resumen en la tabla 1. Tabla 1: S´ ımbolos empleados. Variables electromec´ anicas iIntensidad IAmplitud corriente vVoltaje Vcc Alimentaci´ on inversor feFrecuencia TPeriodo α−βPlano principal x−yPlano de p´ erdidas ωmVelocidad angular ∆Distorsi´ on arm´ onica uEstado inversor UConjunto de estados posibles del inversor NfN´ umero de fases Sub´ ındices y super´ ındices sEst´ ator ∗Referencia rRotor ˆ Predicci´ on 1 Componente fundamental o´ Optimo Control predictivo eError de control KT´ ermino penalizable derivado del error LT´ ermino penalizable no derivado del error JFunci´ on de coste λPeso de ponderaci´ on de la funci´ on de coste ΛVector de pesos λ TmPeriodo de muestreo k´ Indice de tiempo discreto zEstado medible del modelo Ac,BcMatrices ecuaciones de estado Operaci´ on xPunto de funcionamiento XConjunto de los puntos de funcionamiento XmMallado o muestra del conjunto X NxN´ umero de muestras en Xm (Fi,Ij) Mallado para experimentos ´ Indices y medidas ECM Error cuadr´ atico medio NC Conmutaciones en el inversor en un periodo FMC Frecuencia media de conmutaci´ on γFactor el escalado de magnitudes IC ´ Indice de comportamiento δFactor para ponderar los IC IPC ´ Indice ponderado de comportamiento 2.1. Antecedentes y Definiciones En la figura 1 se muestra el diagrama de bloques del control predictivo de corriente de una m´ aquina de inducci´ on cuyos elementos principales son: la m´ aquina (MI), el inversor y el controlador digital (CPCE). Dicho controlador utiliza las mediciones de corrientes de est´ ator (is∈RN f siendo Nfel n´ umero de fases) y de velocidad angular ωmjunto con una referencia para las corrientes i∗ sque proviene del bucle externo de control de velocidad y/o par. El controlador debe calcular en cada instante de muestreo el vector u∈BN f que contiene el estado del inversor. Cada componente de ues un valor binario que indica si el correspondiente interruptor del inversor est´ a abierto (valor 0) o cerrado (valor 1). Para cada fase del inversor, el interruptor M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 49 conjugado toma l´ ogicamente el valor negado. Los voltajes vresultantes de la disposici´ on tomada por el inversor imprimen en el sistema el´ ectrico un cierto comportamiento. Este comportamiento es indicado por las trayectorias de corrientes de est´ ator, que a su vez son responsables de la aparici´ on de otras magnitudes de inter´ es como el flujo, el par generado y, finalmente, la velocidad mec´ anica. Por tanto, controlando la forma de onda de las corrientes de est´ ator se consigue controlar el comportamien- to del sistema. Predicción MI ~ FC u(k) CPCE Optimización is(k) Vcc V is*(k+2) Nf Figura 1: Diagrama del control predictivo de corriente (CPCE) de una m´ aquina de inducci´ on (MI) de Nffases. En cada periodo de muestreo el control predictivo de corriente calcula u(k) mediante la optimizaci´ on combinatoria de los posibles estados del inversor, es decir: uo(k)=argmin u∈U J(k,u) (1) siendo Uel conjunto de estados del inversor y Jla funci´ on de coste que depende de la acci´ on tomada. En el caso de las m´ aquinas de inducci´ on, las isse convierten mediante la transformaci´ on de Clarke en varios subespacios: un plano α−β(relacionado con la generaci´ on de par), varios planos x−y(relacionados con p´ erdidas) y un eje z. Esta descomposici´ on depende de las disposici´ on interna de la m´ aquina y su conexionado externo con el inversor (Duran et al., 2007). En el caso de la m´ aquina pentaf´ asica s´ olo existe un plano x−y. La funci´ on de coste suele incluir como t´ erminos principales aquellos que penalizan la desviaci´ on entre los valores de referencia y los predichos para cada uno de los ejes. En el caso del plano α−βla referencia viene fijada por el bucle exterior de control de velocidad/par. En el caso de los planos x−yla referencia suele ser cero excepto en casos en los que se requiera la inyecci´ on de arm´ onicos. ´ Este es el caso del art´ ıculo (Arahal y Dur´ an, 2009) donde se aborda la sinton´ ıa de los PI para inyecci´ on de tercer arm´ onico. Se define el error de control de corriente predicho para cada uno de los ejes h∈ {α, β, x1,y1,x2,y2, ...} como ˆeh(k+2|k)=i∗ s,h(k+2) −ˆ is,h(k+2|k,u(k)) (2) donde se ha hecho constar expl´ ıcitamente que la predicci´ on a dos pasos es una funci´ on de la acci´ on de control (decidida en el instante k) que se calcula mediante un modelo del sistema. La doble predicci´ on es una t´ ecnica habitual que resulta del uso de un tiempo de muestreo cercano al tiempo de c´ alculo necesario para obtener u(v´ ease (Arahal et al., 2009)). Consiste en aplicar la se˜ nal de control tras dejar transcurrir todo el tiempo de muestreo actual, por lo tanto dicha se˜ nal no producir´ a efectos medibles hasta el siguiente periodo de muestreo por lo que se necesita una predicci´ on a dos pasos. Con las anteriores definiciones se puede ya presentar la estructura t´ ıpica de una funci´ on de coste del siguiente modo: J(k,u)=Kαβ(k,u)+λK1Kx1y1(k,u)+... +λL1L1(k,u)+... (3) donde los distintos t´ erminos Klm son errores (de control predichos) en distintos planos l−m(de la transformaci´ on de Clarke) y se definen mediante Klm(k,u)=k(ˆel(k+2|k),ˆem(k+2|k))k2(4) mientras que los t´ erminos Licorresponden a otras magnitudes penalizables, como por ejemplo el n´ umero de conmutaciones. Los coeficientes λson pesos de ponderaci´ on que se emplean habitualmente en optimizaci´ on multiobjetivo pues permiten ponderar magnitudes diversas dentro de una misma funci´ on de coste. La tarea del controlador consiste en resolver (1) y aplicar y mantener la configuraci´ on del inversor uodurante un periodo de muestreo, repitiendo el proceso al siguiente periodo de muestreo mediante la t´ ecnica del horizonte deslizante (Camacho y Bordons, 2013). El empleo de esta funci´ on de coste con un modelo suficientemente preciso ((Rodas et al., 2016), (Mart´ ın et al., 2017)) permite minimizar el error de control y los otros factores considerados. Ha de tenerse en cuenta que en la mayor´ ıa de casos la soluci´ on de (1) es un compromiso entre criterios opuestos, pues no siempre existe un estado del inversor que reduce a la vez todos los t´ erminos KyL. 2.2. Dise˜no Habitual de Funciones de Coste La forma habitual de dise˜ nar la funci´ on de coste consiste en seleccionar qu´ e t´ erminos KyLse van a tener en cuenta y luego asignar un valor a los pesos de ponderaci´ on. La primera parte obedece a criterios espec´ ıficos de la aplicaci´ on de la m´ aquina (generaci´ on, tracci´ on, etc.) los cuales determinan hasta qu´ e punto es importante penalizar el error en distintos planos y otras variables (conmutaciones por ejemplo). Adem´ as puede ocurrir que se deba acotar la complejidad de la funci´ on de coste para garantizar que los c´ alculos se llevar´ an a cabo en un periodo de muestreo con el procesador disponible (Arahal et al., 2009). La segunda parte (valor de los pesos de ponderaci´ on) se ha realizado normalmente mediante exploraci´ on (experimental o por simulaci´ on) de algunos puntos de operaci´ on. El valor hallado se considera un par´ ametro del controlador y por tanto se mantiene fijo durante la operaci´ on. La idea de que los pesos de ponderaci´ on no son ´ optimos para todos los puntos de operaci´ on ha aparecido en la literatura, sin embargo no se ha presentado a´ un un control predictivo de corriente que utilice otra t´ ecnica diferente de la habitual. A continuaci´ on se detalla la propuesta de este trabajo. 2.3. Ponderaci´on Optimizada para cada Punto de Operaci´on La t´ ecnica para el dise˜ no de funci´ on de coste que se propone consiste en usar la forma funcional indicada en (3) pero tomando los pesos de ponderaci´ on como variables que dependen del 50 M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 punto de operaci´ on en la forma que se indicar´ a m´ as adelante. Por ahora puede considerarse que el punto de operaci´ on queda definido por el estado electro-mec´ anico: x(k)=(is(k),ir(k), ωm(k),s(k))(5) siendo irla corriente de rotor, ωmla velocidad angular y sel deslizamiento. De este modo se llega a la expresi´ on para la funci´ on de coste propuesta J(k,u)=Kαβ(k,u)+λK1(k)Kx1y1(k,u)+... +λL1(k)L1(k,u)+... (6) donde la dependencia del punto de operaci´ on del instante kse hace constar mediante dicho ´ ındice; es decir, se ha abreviado λ(x(k)) como λ(k). En la funci´ on de coste se deben incluir todos los t´ erminos que se consideren importantes, con la ´ unica limitaci´ on del tiempo de c´ alculo disponible. El proceso de determinaci´ on de las funciones λ(x(k)) asignar´ a valores que ponderen dichos t´ erminos incluyendo la posibilidad de anular alguno de ellos. Para determinar los pesos de ponderaci´ on se tiene en cuenta que el comportamiento deseado de la m´ aquina se revela al considerar, no muestras, sino intervalos de funcionamiento. A modo de ejemplo consid´ erese la distorsi´ on arm´ onica (∆) de la corriente de est´ ator, definida como la raz´ on de los valores eficaces de la componente principal I1y de la corriente total (componente fundamental +arm´ onicos); es decir: ∆ = I1/I. Por su propia definici´ on el valor eficaz no tiene sentido en un ´ unico periodo de muestreo, m´ as a´ un, las componentes arm´ onicas (incluyendo la fundamental) precisan de algunos periodos de la frecuencia fundamental (T=1/fe) para poder ser estimadas. Por tanto, para hablar con propiedad del comportamiento de la m´ aquina hay que considerar distintas variables cuya definici´ on y c´ alculo abarca varios T. Los anteriormente citados ´ ındices de comportamiento pueden estimarse mediante simulaciones o mediante experimentaci´ on. En todo caso los experimentos (reales o simulados) deber´ ıan ser representativos de todo el conjunto de puntos de funcionamiento del sistema. Adem´ as, los ´ ındices de comportamiento pueden ser contradictorios entre s´ ı. Por ejemplo, al mejorar el contenido arm´ onico posiblemente se empeore la frecuencia media de conmutaci´ on. Teniendo esto presente, parece l´ ogico incorporar los distintos ´ ındices de comportamien- to en una funci´ on de coste que tenga en cuenta la distinta importancia que se le concede a cada uno. Para evitar confusiones con la funci´ on de coste del control predictivo de corriente se va a nominar esta nueva funci´ on de coste como ´ ındice ponderado de comportamiento. Al igual que la funci´ on de coste, el ´ ındice ponderado de comportamiento no es m´ as que un modo de lidiar con un problema de optimizaci´ on multi-objetivo. En este caso lo que se intenta determinar son los pesos de ponderaci´ on. En la literatura de control existen propuestas con ideas similares, por ejemplo la utilizaci´ on de diferentes par´ ametros de un PID en funci´ on del punto de operaci´ on. Los ´ ındices de funcionamiento como integral del error, sobreoscilaci´ on, etc. permiten en estos casos sintonizar dichos par´ ametros (Wittenmark, 1975). Para describir el problema en forma matem´ atica se van a agrupar todas las funciones λen un vector funcional Λ = (λK1, λK2, ..., λL1, ...)|(7) y se buscar´ a el vector que optimiza el ´ ındice ponderado de comportamiento IPC definido como IPC(Λ)= Ni X n=1 δnICn(Λ) (8) donde los ICnson los ´ ındices de comportamiento que se desean tener en cuenta y los δnson factores de ponderaci´ on. As´ ı las cosas parece natural seleccionar los pesos de ponderaci´ on mediante Λo=argmin Λ IPC(Λ) (9) Para ello en primer lugar se han de definir los par´ ametros δ que ponderan los ´ ındices de comportamiento en el ´ ındice ponderado de comportamiento. Esta tarea requiere conocer cu´ anto se est´ a dispuesto a sacrificar en un ´ ındice para mejorar otro. Mediante una correcta normalizaci´ on de variables se puede obtener una aproximaci´ on para los δque puede posteriormente retocarse. De manera intuitiva un aumento de un cierto δnprovocar´ a una mejora de ICna expensas de un posible empeoramiento de los otros ´ ındices (ICr/r,n). Es importante resaltar que se ha pasado de una situaci´ on en que los pesos de ponderaci´ on afectan de manera poco clara al comportamiento (dificultando su elecci´ on), a esta nueva situaci´ on en que los δest´ an ligados a un determinado ´ ındice por lo que se puede intuir su efecto sobre el comportamiento del sistema en bucle cerrado. Para obtener Λomediante (9) se ha de resolver un problema de optimizaci´ on para cada punto de operaci´ on. Se hace necesario por tanto muestrear el espacio Xde los puntos de operaci´ on x. Para ello se puede recurrir a dise˜ nar el muestreo para recorrer una cierta cuadr´ ıcula que abarque el intervalo de variaci´ on de cada variable o se puede recurrir a un algoritmo como k-medias. En el primero de los casos los valores quedan fijados por los nodos de la cuadr´ ıcula mientras que en el segundo la distribuci´ on es el resultado de un algoritmo auto-organizado. En el caso de que las variables a tener en cuenta no sean un n´ umero muy alto (1 a 4) el primer m´ etodo es factible y m´ as intuitivo. En el caso de las m´ aquinas polif´ asicas las diferentes cantidades est´ an relacionadas entre s´ ı por medio de la din´ amica del sistema. Por tanto el punto de operaci´ on queda definido con unas pocas cantidades y puede aplicarse el m´ etodo basado en la cuadr´ ıcula. En el siguiente apartado se ilustrar´ a este aspecto con el caso de una m´ aquina pentaf´ asica. Sea cual sea el m´ etodo de muestreo elegido, se obtendr´ a un conjunto de pares (xi,IPC(Λ,xi)) donde los puntos de operaci´ on xipara i=1, ..., Nxforman el conjunto Xm={xi}con Xm⊂X. Restringiendo la optimizaci´ on a estos pares, se obtiene Λo i=argmin Λ IPC(Λ,xi) (10) que permite la creaci´ on de un nuevo conjunto de pares (xi,Λo i) para i=1, ..., Nx. A partir de estos pares es posible crear una funci´ on que abarque todos los puntos de operaci´ on. Para ello existen muchas t´ ecnicas posibles, siendo la m´ as simple la interpolaci´ on multidimensional (Kirkland, 2010). Como ejemplo de otras t´ enicas cabe citar la reconstrucci´ on de funciones por M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 51 polinomios aproximantes, mediante series de Fourier, mediante redes de neuronas artificiales, etc. (Rivlin, 2003). Sea cual sea la t´ ecnica empleada se llega a unas expresiones para las componentes del vector Λoque dependen del punto de operacion actual x(k). Por tanto se pueden utilizar en tiempo real para obtener los pesos de ponderaci´ on del control predictivo de corriente. Estos valores sustituir´ ıan a los valores gen´ ericos consignados en (3). 3. Aplicaci´ on a M´ aquina de Inducci´ on Pentaf´ asica Tomando las asunciones t´ ıpicas para el modelado de m´ aquinas de inducci´ on, a saber: uniformidad del entrehierro, distribuci´ on magn´ etica senoidal, bobinados sim´ etricamente distribuidos, p´ erdidas en el n´ ucleo cero y no existencia de saturaci´ on magn´ etica, el modelo en tiempo continuo de la m´ aquina toma la forma ˙zc=Ac(ωr)·zc+Bc·u(11) en el que se ha usado como estado zc zc=is,α,is,β,is,x,is,y,ir,α,ir,β|(12) siendo isel vector de corrientes de est´ ator, irel vector de corrientes del rotor y ωrla velocidad el´ ectrica del rotor. Las matrices AcyBcdependen de los par´ ametros el´ ectricos del sistema, como resistencias e inductancias. Para mayor informaci´ on acerca del modelado se remite al lector a (Mart´ ın et al., 2017). La anterior descripci´ on (11) contiene magnitudes no medibles pues pertenecen al rotor. A fin de derivar un modelo predictivo se han de estimar dichas magnitudes (o su efecto) como se describe en (Mart´ ın et al., 2016). De este modo se llega a un nuevo modelo en tiempo discreto z(k+1) =F(ωr)·z(k)+G(k)+H·u(k) (13) en el que el vector Gtiene en cuenta el efecto de las variables del rotor sobre el estado medible, definido mediante z(k)=is,α,is,β,is,x,is,y|(k) (14) Las predicciones a dos pasos necesarias para (3), indicadas mediante ˆ is,h(k+2|k,u(k)) for h∈ {α, β, x,y}, se obtienen de manera inmediata a partir de (13) como se muestra en (Arahal et al., 2009). 3.1. Factores de la Funci´on de Coste Para la aplicaci´ on a la m´ aquina de cinco fases se va a seguir el caso presentado en (Barrero et al., 2009) por lo que la funci´ on de coste toma la forma J(k)=kˆeαβ(k)k2+λ1kˆexy(k)k2+λ2NC(k) (15) siendo NC el n´ umero de conmutaciones producidos al aplicar u(k) (partiendo del estado previo del inversor dado por uo(k)), que puede calcularse como NC(k)=γ N f X i=1ui(k+1) −uo i(k)(16) siendo γun factor para tener en cuenta que los valores t´ ıpicos del n´ umero de conmutaciones es muy diferente de los valores t´ ıpicos de los errores y se elige como. γ= Ehkˆek2i EhPN f i=1ui(k+1) −uo i(k)i(17) siendo Eel operador esperanza matem´ atica. Con esta funci´ on de coste se tienen en cuenta varios de los factores que afectan al comportamiento de la m´ aquina pues los arm´ onicos est´ an relacionados con kˆeαβkmientras que las p´ erdidas lo est´ an con kˆexyk y con NC. 3.2. ´ Indices de Comportamiento Para el caso tomado como ejemplo se van a considerar, sin p´ erdida de generalidad, los indicadores t´ ıpicos derivados del error cuadr´ atico medio de control en el plano α−β, del plano x−yy de la frecuencia media de conmutaci´ on del inversor. El error cuadr´ atico medio de control se calcula para un eje gen´ erico hmediante ECMh=v u t1 N N X k=1i∗ s,h(k)−is,h(k)2(18) y la frecuencia media de conmutaci´ on (por fase y por ciclo) como FMC =1 N·Tm·fe N X k=1 NC(k) (19) siendo Nel n´ umero de muestras considerado para el c´ alculo. Las anteriores cantidades se han de agrupar por planos y normalizar al intervalo [0,1], dividiendo por los valores m´ aximos, obteni´ endose IC1=ECMα+ECMβ1 2 m´ ax ECMα,β (20) IC2=ECMx+ECMy1 2 m´ ax ECMx,y (21) IC3=FMC m´ ax FMC (22) Con estas definiciones se puede escribir el ´ ındice ponderado de comportamiento usado para este ejemplo como IPC =δ1IC1+δ2IC2+δ3IC3(23) donde los factores de ponderaci´ on δiser´ an especificados en cada experimento a fin de ilustrar las capacidades del m´ etodo propuesto. T´ engase presente que la elecci´ on realizada permite desarrollar un caso t´ ıpico para ilustrar la t´ ecnica propuesta. El lector interesado es remitido a (Vukosavi´ c et al., 2017) donde se especifican otros ´ ındices. Conviene indicar que si se elimina un ´ ındice de comportamiento entonces su valor queda libre, pudiendo empeorar. Esto afecta a la soluci´ on pues dicha eliminaci´ on es similar a eliminar una restricci´ on de un problema de optimizaci´ on. Por contra, si se a˜ nade un nuevo ´ ındice de comportamiento se tiene un caso similar a a˜ nadir una restricci´ on adicional a un 52 M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 problema de optimizaci´ on. Se consigue acotar el ´ ındice a˜ nadido a costa de un posible empeoramiento de los que previamente estaban. 3.3. Optimizaci´on de la Funci´on de Coste Para la optimizacion de la funci´ on de coste se seleccionan en primer lugar el conjunto de variables que definen un punto de operaci´ on. Debido a que la m´ aquina puede funcionar a distintas velocidades y con distinta carga es preciso tomar en consideraci´ on la frecuencia el´ ectrica fey la amplitud de las corrientes de referencia de est´ ator I∗. Es f´ acil ver que, debido a las ecuaciones del modelo del sistema, estas variables tienen gran influencia en la din´ amica de la m´ aquina. Es posible considerar otras variables adicionales (como el deslizamiento) con lo que se consigue una discretizaci´ on m´ as fina de los puntos de operaci´ on. A cambio se producir´ ıa un aumento de la carga de c´ alculo necesaria para evaluar las funciones Λ(x(k)) =(λ1(x(k)), λ2(x(k)))|. Por tanto se van a usar ´ unicamente feeI∗para definir el punto de operaci´ on y se toma x=(fe,I∗). El muestreo del espacio de puntos de operaci´ on Xse lleva a cabo considerando una cuadr´ ıcula que discretiza dicho espacio. Por sencillez se puede definir la cuadr´ ıcula como el produc- to cartesiano Xm=SF×SIsiendo SF={F1,F2, ..., FNF }un conjunto ordenado de frecuencias y siendo SI={I1,I2, ..., INI } un conjunto ordenado de intensidades. El n´ umero de muestras considerado es card(Xm)=NF·NIy aparecen distribuidas en las coordenadas (Fi,Ij)∀i=1, ..., NF,j=1, ..., NI. Para cada punto de operaci´ on definido por los ´ ındices (i,j) (x=(Fi,Ij)) se obtiene la soluci´ on de: λo 1(i,j) λo 2(i,j)!=argmin Λ IPC(Λ,(Fi,Ij)) (24) expresi´ on que proporciona los valores ´ optimos (seg´ un el ´ ındice ponderado de comportamiento) en los puntos de operaci´ on correspondientes a la cuadr´ ıcula Xm. Durante el funcionamiento de la m´ aquina es probable que se visiten otros puntos de operaci´ on diferentes a los considerados en Xm, por lo que es necesario proporcionar un valor de los pesos de ponderaci´ on para un caso general. Esto se puede lograr mediante diversas t´ ecnicas, por ejemplo mediante la interpolaci´ on bilineal. Esta t´ ecnica no es m´ as que la extensi´ on de la interpolaci´ on lineal para el caso de funciones de dos variables. Para ello se emplea una malla regular de dos dimensiones (Kirkland, 2010). El m´ etodo constituye una t´ ecnica est´ andard que se aplica f´ acilmente del modo siguiente: si el punto de operaci´ on en el instante de tiempo discreto kes x(k)=(fe(k),I∗(k)), para interpolar se buscar´ an dos ´ ındices i1,i2tales que Fi1≤fe(k)≤Fi2con i1,i2y tal que Fi2−Fi1sea m´ ınimo. Del mismo modo se buscar´ an dos ´ ındices j1yj2tales que Ij1≤I∗(k)≤Ij2con j1,j2y tal que Ij2−Ij1sea m´ ınimo. Los pares de ´ ındices (i1,j1), (i1,j2), (i2,j2) y (i2,j1) definen los v´ ertices de la celda de la cuadr´ ıcula que contiene al punto de operaci´ on actual. La interpolaci´ on bilineal permite generar un valor para los λ(x(k)) a partir de los valores de λoen los cuatro pares antecitados. Mediante simple sustituci´ on de estos valores en la f´ ormula (Kirkland, 2010) y considerando las siguientes abreviaturas: ∆I20 =(Ij2−I(k)), ∆I01 =(I∗(k)−Ij1), ∆F20 =(Fi2−fe(k)), ∆F01 =(fe(k)−Fi1) y λo rs =λo q(ir,js), se obtiene la siguiente aproximaci´ on para el valor λ(x(k)): λ(x(k)) = (λo 11∆I20 +λo 12∆I01)∆F20 +(λo 21∆I20 +λo 22∆I01)∆F01 (∆I20 + ∆I01)(∆F20 + ∆F01)(25) Como es l´ ogico habr´ a de aplicarse (25) a cada componente λqq∈ {1,2}. 4. Simulaciones Se presentan a continuaci´ on los resultados obtenidos por simulaci´ on. El sistema usado en las simulaciones es una m´ aquina pentaf´ asica con las siguientes caracter´ ısticas: resistencia de est´ ator 19.4 (Ω), resistencia del rotor 6.8 (Ω), inductancia de est´ ator 101 (mH), inductancia del rotor 39 (mH), inductancia mutua 657 (mH), velocidad nominal 1000 (rpm), potencia 1 (kW), pares de polos 3. Se ha simulado el controlador pertinente en cada caso suponiendo un periodo de muestreo de 80 (µs) que es un valor habitualmente encontrado y bastante asequible para procesadores digitales de se˜ nal actuales. En primer lugar se presentan resultados correspondientes a un control predictivo de corriente con pesos de ponderaci´ on fijos (correspondientes a la t´ ecnica habitual) y que servir´ an de base de comparaci´ on. En la tabla 2 se muestran los resultados para un par de puntos de funcionamiento O1yO2correspondientes a una carga media y alta respectivamente. Concretamente, para O1se usa fe=40 (Hz) e I∗=1 (A), y para O2se usa fe=60 (Hz) e I∗=2,5 (A). Se puede obervar el efecto de los distintos valores de pesos de ponderaci´ on considerados (indicados en la primera columna de cada fila) sobre los ´ ındices de comportamiento y c´ omo el efecto no es id´ entico para cada punto de operaci´ on. Tabla 2: Resultados con la t´ ecnica habitual. O1O2 (λ1, λ2)IC1IC2IC3IC1IC2IC3 ( 0.02, 0.02 ) 0.13 0.77 0.91 0.27 0.87 0.52 ( 0.50, 0.02 ) 0.25 0.31 0.54 0.33 0.35 0.68 ( 0.02, 0.10 ) 0.16 1.09 0.52 0.27 0.99 0.40 ( 0.50, 0.10 ) 0.26 0.32 0.47 0.34 0.36 0.66 A modo de introducci´ on al siguiente apartado de este art´ ıculo, es interesante observar que, de los cuatro valores de pesos de ponderaci´ on considerados en la tabla 2 y tomando δ1=1, δ2=0,7 y δ3=0,4, el menor valor de ´ ındice ponderado de comportamiento para O1corresponde a (λ1, λ2)=(0.50, 0.10) mientras que para O2corresponde a (λ1, λ2)=(0.50, 0.02). Esta observaci´ on apunta a la posibilidad de optimizar los pesos de ponderaci´ on en funci´ on del punto de operaci´ on como se muestra a continuaci´ on. 4.1. Resultados Optimizando la Funci´on de Coste Se presentan a continuaci´ on los resultados obtenidos al considerar la t´ ecnica propuesta. Se utiliza la aproximaci´ on bilineal indicada anteriormente para obtener los valores de λ1yλ2acor- des a cada punto de operaci´ on. M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 53 Figura 2: Valores ´ optimos de λ1(arriba) y λ2(abajo) obtenidos mediante la resoluci´ on de (10) para cada punto de operaci´ on. N´ otese la escala logar´ ıtmica empleada. Para las pruebas que siguen se ha seleccionado SFcon valores entre 20 (Hz) y 70 (Hz) tomados cada 10 (Hz). De manera similar se ha considerado SIcon valores entre 0,4 (A) y 3 (A) tomados cada 0,15 (A). De este modo se construye Xmque cubre el espacio de operaci´ on para esta m´ aquina particular. Para la optimizaci´ on se han usado los valores δ1=1, δ2=0,7 y δ3=0,4 que corresponden a un cierto compromiso entre los distintos ´ ındices de comportamiento. En la figura 2 se muestran los resultados de (24). Los ejes corresponden a las dos dimensiones consideradas para definir el punto de operaci´ on: feeI∗, mientras que los valores de λ1yλ2 se han codificado como un color seg´ un la leyenda. Se ha usado escala logar´ ıtmica debido a la disparidad de valores hallados. As´ ımismo, los valores que no caen encima de los nodos de la cuadr´ ıcula se han aproximado mediante (25). Puede verse que los valores ´ optimos son bastante distintos de un punto de operaci´ on a otro. Esto valida la hip´ otesis inherente al m´ etodo presentado, que es la dependencia de los pesos de ponderaci´ on con el punto de operaci´ on. Tambi´ en puede verse que la dependencia no es trivial. Esto se pone de manifiesto por la no existencia de patrones claros en la variaci´ on de λ1yλ2por el espacio de operaci´ on. Este descubrimiento valida el uso de simulaci´ on intensiva y la posterior aproximaci´ on de funciones. Como siguiente resultado se presenta en la figura 3 (arriba) el valor obtenido para el ´ ındice ponderado de comportamiento para cada punto de funcionamiento al utilizar el m´ etodo propuesto. Figura 3: Valores del índice ponderado de comportamiento con la técnica propuesta (arriba) y con valores fijos de los pesos de ponderación (abajo) para cada punto de operación. Puede verse que el índice ponderado de comportamiento mantiene unos valores bajos en todo el espacio de operación con una pequeña subida a medida que se progresa hacia el límite de la máquina para cargas muy altas. Esto contrasta con el resultado obtenido al usar valores fijos de los pesos de ponderación, como puede verse en 3 (abajo) donde se han usado valores típicos λ1 = 0,5 y λ2 = 0,02. Se puede observar: un mayor y más temprano deterioro en la zona de alta carga, mayores valores de índice ponderado de comportamiento en la totalidad del espacio de trabajo y en particular en una franja diagonal hacia el centro del espacio de trabajo. Esta situación indica una mejora sustancial (en términos de rendimiento energético) y global (para todo el espacio de trabajo) a favor del método propuesto. Como complemento a la gráfica, se incluye la Tabla 3 donde se indican algunos valores que permiten comparar los resultados de la técnica propuesta y del procedimiento habitual. Los valores indicados son el mínimo, máximo, media (µ) y desviación típica (σ) al considerar todo el espacio de trabajo. Es destacable que, con la técnica propuesta, se consiguen menores valores en casi todos los casos y además con menor variación en el espacio de trabajo (menor valor de la desviación típica). Merece la pena insistir en que, modificando los valores δ se puede orientar el diseño. Por ejemplo para tener menores pérdidas en el plano x − y se puede aumentar δ2. De manera parecida se puede orientar el diseño a menores pérdidas por conmutación, siendo necesario en este caso un aumento de δ3 o a mejor seguimiento en α − β (aumentando δ1). 54 M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 Tabla 3: Comparaci´ on de los valores del ´ ındice ponderado de comportamiento para la t´ ecnica habitual (TH) y propuesta (TP). TH m´ ınimo m´ aximo µ σ ECMαβ (A) 0.023 0.229 0.032 0.025 ECMxy (A) 0.029 0.172 0.037 0.017 FMC (kHz) 0.944 6.114 4.539 1.333 IPC 0.294 1.900 0.564 0.189 TP m´ ınimo m´ aximo µ σ ECMαβ (A) 0.018 0.047 0.033 0.005 ECMxy (A) 0.020 0.194 0.037 0.034 FMC (kHz) 0.773 4.237 3.044 0.843 IPC 0.294 1.287 0.494 0.122 5. Conclusiones Se ha puesto de manifiesto la conveniencia de usar distintos factores de ponderaci´ on para cada punto de funcionamiento. Se ha planteado un m´ etodo de optimizaci´ on como medio de obtener los mejores valores para los factores de ponderaci´ on. Los resultados muestran que es posible usar la interpolaci´ on bilineal para abarcar todo el espacio de operaci´ on. Los resultados por simulaci´ on avalan la factibilidad del m´ etodo propuesto. Agradecimientos Este trabajo ha sido financiado por el proyecto DPI2016- 76144-R del ministerio de econom´ ıa y competitividad de Espa˜ na. Referencias Abu-Rub, H., Guzinski, J., Krzeminski, Z., Toliyat, H. A., 2001. Predictive current control of voltage source inverters. In: Industrial Electronics Society, 2001. IECON ’01. The 27th Annual Conference of the IEEE. Vol. 2. pp. 1195–1200 vol.2. Arahal, M. R., Barrero, F., Toral, S., Dur´ an, M. J., Gregor, R., 2009. Multi-phase current control using finite-state model-predictive control. Control Engineering Practice 17 (5), 579 – 587. th Barrero, F., Arahal, M. R., Gregor, R., Toral, S., Durán, M. J., 2009. Onestep modulation predictive current control method for the asymmetrical dual three-phase induction machine. Industrial Electronics, IEEE Transactions on 56 (6), 1974–1983. Camacho, E. F., Bordons, C., 2013. Model predictive control. Springer. Duran, M., Toral, S., Barrero, F., Levi, E., 2007. Real-time implementation of multi-dimensional five-phase space vector pulse-width modulation. Electronics Letters 43 (17), 949–950. Holmes, D., Martin, D., 1996. Implementation of a direct digital predictive current controller for single and three phase voltage source inverters. In: Industry Applications Conference, 1996. Thirty-First IAS Annual Meeting, IAS’96., Conference Record of the 1996 IEEE. Vol. 2. IEEE, pp. 906–913. Holtz, J., Stadtfeld, S., 1983. A predictive controller for the stator current vector of AC machines fed from a switched voltage source. In: JIEE IPEC-Tokyo Conf. pp. 1665–1675. Kirkland, E. J., 2010. Bilinear interpolation. In: Advanced Computing in Electron Microscopy. Springer, pp. 261–263. Lim, C.-S., Levi, E., Jones, M., Rahim, N., Hew, W.-P., Aug 2014. A comparative study of synchronous current control schemes based on FCS-MPC and PI-PWM for a two-motor three-phase drive. Industrial Electronics, IEEE Transactions on 61 (8), 3867–3878. Mamdouh, M., Abido, M., Hamouz, Z., 2018. Weighting factor selection techniques for predictive torque control of induction motor drives: A comparison study. Arabian Journal for Science and Engineering, 1–13. Mart´ ın, C., Arahal, M. R., Barrero, F., Dur´ an, M. J., 2016. Five-phase induction motor rotor current observer for finite control set model predictive control of stator current. IEEE Transactions on Industrial Electronics 63 (7), 4527– 4538. Mart´ ın, C., Berm´ udez, M., Barrero, F., Arahal, M. R., Kestelyn, X., Dur´ an, M. J., 2017. Sensitivity of predictive controllers to parameter variation in five-phase induction motor drives. Control Engineering Practice 68, 23–31. Rivlin, T. J., 2003. An introduction to the approximation of functions. Courier Corporation. Rodas, J., Barrero, F., Arahal, M. R., Mart´ ın, C., Gregor, R., 2016. Online estimation of rotor variables in predictive current controllers: A case study using five-phase induction machines. IEEE Transactions on Industrial Electronics 63 (9), 5348–5356. Rojas, C. A., Rodriguez, J., Villarroel, F., Espinoza, J. R., Silva, C. A., Trincado, M., 2013. Predictive torque and flux control without weighting factors. IEEE Transactions on Industrial Electronics 60 (2), 681–690. Vukosavi´ c, S. N., Peri´ c, L. S., Levi, E., 2017. A three-phase digital current controller with improved performance indices. IEEE Transactions on Energy Conversion 32 (1), 184–193. Wittenmark, B., 1975. Stochastic adaptive control methods: a survey. International Journal of Control 21 (5), 705–730. Zoric, I., Jones, M., Levi, E., 2018. Arbitrary power sharing among three-phase winding sets of multiphase machines. IEEE Transactions on Industrial Electronics 65 (2), 1128–1139. M. Arahal et al. / Revista Iberoamericana de Automática e Informática Industrial 16 (2019) 48-55 55 harmonic injection. Control Engineering Practice 17 (7), 787–797 Arahal, M. R., Dur´ an, M. J., 2009. PI tuning of five-phase d rives w ith third