MIMO MPC
4.2 (84.44%) 9 votes

En esta entrada veremos un poco el funcionamiento de un controlador predictivo basado en modelo (GPC) utilizados para controlar procesos multivariables o por sus siglas en ingles MIMO (Multi-Input Multi-Output). Una representación basica en diagrama de bloques para un sistema MIMO 2×2 puede verse en la siguiente figura:

sistema MIMO

En este ejemplo especifico vemos un sistema que tiene 2 entradas (u1 y u2) y dos salidas (y1 y y2). Cada una de esas entradas puede afectar en cierta proporción a cada una de las salidas del sistema.

 

CARACTERÍSTICAS DEL MPC MIMO

 

  • No se preocupa por realizar desacoplo de variables ya que trata el problema de forma general
  • Se deben escoger adecuadamente los pares entrada-salida (Una técnica para decidir esto es la matriz RGA)
  • Es importante escalar adecuadamente TODAS las variables del sistema, debido a que normalmente se trabaja con distintas variables (Temperatura, Presión, Nivel, Flujo, etc)

 

FUNCION OBJETIVO

 

La función objetivo viene siendo igual a como hemos venido trabajando en el caso SISO. Donde se tiene el seguimiento de referencia y las acciones de control, todo eso ponderado. La diferencia es que ahora las ponderaciones no son escalares, dado a que tengo más de una variable en mi sistema. Y  para cada salida voy a tener horizontes diferentes.

 

J = \mathop \sum \limits_{k = {N_1}}^{{N_2}} ||\hat y\left( {t + k} \right) - W(t + k)||{^2}{Q_1} + \mathop \sum \limits_{K = 1}^{{N_3}} ||\Delta u(t + k - 1)||{^2}{Q_2}

 

La propiedad aplicada en la función objetivo es una Norma 2.

 

||X||_{{Q_1}}^2 = {X^T}{Q_1}X

 

En otras palabras como para cada salida voy a tener diferentes horizontes, suponiendo a modo de ejemplo que tengo dos salidas (y1 que va desde el horizonte N11 hasta N12) y (y2 que va desde el horizonte N21 hasta N22). Mi matriz de predicciones “y” va a contener las dos salidas, es decir que la salida estará dividida por sectores y que para cada sector voy a tener que ponderar de forma diferente, quiere decir que mi vector de ponderación también estará dividido por sectores:

ponderacion MIMO

Y se hace exactamente lo mismo con las ponderaciones y con el vector del incremento de control, el cual también estará dividido por sectores:

MIMO GPC

En forma general de la función objetivo podemos ver que los tamaños de cada vector o matriz viene dado por:

 

J = \mathop \sum \limits_{k = {N_1}}^{{N_2}} ||\hat y\left( {t + k} \right) - W(t + k)||{^2}{Q_1} + \mathop \sum \limits_{K = 1}^{{N_3}} ||\Delta u(t + k - 1)||{^2}{Q_2}

 

  • Vector de predicción tamaño = Suma(N_{2i}-N_{1i})
  • Vector de control tiene tamaño = Suma(N_{ui})
  • Matriz de Ponderación Q_1 y Q_2 bloque diagonal

Para montar la matriz de coeficientes del sistema, se deben dejar todas las entradas del sistema constantes y aplicar un escalón únicamente en una entrada y ver como se comporta cada una de las salidas del sistema. Para la primera variable manipulada podría calcular los primeros sectores de la matriz G:

Matriz G GPC

Y asi sucesivamente hasta llenar toda la matriz (Este procedimiento sirve para determinar la matriz G tanto para el DMC cuanto para el GPC).

 

Caso GPC

Analicemos como se estructura nuestro modelo MIMO expresado en su forma CARIMA:

 

y\left( t \right) = P\left( {{z^{ - 1}}} \right)u(t)

 

Donde “y” son las salidas, “P” es la función de transferencia MIMO y “u” son las entradas.

 

Cada función de transferencia que conforma el elemento “P” podria expresarse de la siguiente manera:

 

{p_{ij}} = \dfrac{{{z^{ - 1}}{B_{ij}}'({z^{ - 1}})}}{{{A_{ij}}({z^{ - 1}})}}{z^{ - {d_{ij}}}}

 

Y como es un sistema MIMO, es muy probable que cada una de las funciones de transferencia presenten retardos diferentes. Por eso lo más conveniente es tomar para cada salida del sistema MIMO (Cada fila de la matriz P) el mínimo retardo que se presente en esa salida (esa fila):

 

{d_i} = \mathop {\min }\limits_j [{d_{ij}}]

 

Por ejemplo si tenemos la siguiente matriz de transferencia 2×2:

 

\left[ {\begin{array}{*{20}{c}}  {{y_1}(t)}\\  {{y_2}(t)}  \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}  {\dfrac{{0.3776{z^{ - 1}}}}{{1 - 0.9705{z^{ - 1}}}}}&{\dfrac{{ - 0.4447{z^{ - 1}}}}{{1 - 0.9765{z^{ - 1}}}}{z^{ - 2}}}\\  {\dfrac{{0.2959{z^{ - 1}}}}{{1 - 0.9552{z^{ - 1}}}}{z^{ - 2}}}&{\dfrac{{ - 0.6621{z^{ - 1}}}}{{1 - 0.9659{z^{ - 1}}}}{z^{ - 1}}}  \end{array}} \right]\left[ {\begin{array}{*{20}{c}}  {{u_1}(t)}\\  {{u_2}(t)}  \end{array}} \right]

 

Podemos extraer el retardo mínimo de cada salida y agruparlo en la diagonal principal de otra matriz, es decir llevarlo a la siguiente representación:

 

P\left( {{z^{ - 1}}} \right) = D\left( {{z^{ - 1}}} \right)G\left( {{z^{ - 1}}} \right)

 

Donde “P” es la matriz de transferencia, “D” es la matriz con los retardos mínimos y “G” es la matriz de transferencia que contiene los retardos efectivos del sistema:

 

P\left( {{z^{ - 1}}} \right) = \left[ {\begin{array}{*{20}{c}}  1&0\\  0&{{z^{ - 1}}}  \end{array}} \right]\left[ {\begin{array}{*{20}{c}}  {\dfrac{{0.3776{z^{ - 1}}}}{{1 - 0.9705{z^{ - 1}}}}}&{\dfrac{{ - 0.4447{z^{ - 1}}}}{{1 - 0.9765{z^{ - 1}}}}{z^{ - 2}}}\\  {\dfrac{{0.2959{z^{ - 1}}}}{{1 - 0.9552{z^{ - 1}}}}{z^{ - 1}}}&{\dfrac{{ - 0.6621{z^{ - 1}}}}{{1 - 0.9659{z^{ - 1}}}}}  \end{array}} \right]

 

Ahora se debe aplicar el mismo procedimiento que se aplica en el caso de las perturbaciones. Si analizamos bien el problema, veremos que el sistema MIMO es un sistema que presenta interacciones entre variables y que cada interacción, nosotros podriamos considerarla como si fuera una perturbación medible. Osea que podriamos pensar en un desacoplo por FeedForward. Entonces aplicando el mismo concepto, debemos obtener el mínimo común multiplo de cada salida como se muestra a continuación:

 

{A_1}\left( {{z^{ - 1}}} \right) = \left( {1 - 0.9705{z^{ - 1}}} \right)\left( {1 - 0.9765{z^{ - 1}}} \right) = 1 - 1.9470{z^{ - 1}} + 0.9477{z^{ - 2}}

 

{A_2}\left( {{z^{ - 1}}} \right) = \left( {1 - 0.9552{z^{ - 1}}} \right)\left( {1 - 0.9659{z^{ - 1}}} \right) = 1 - 1.9210{z^{ - 1}} + 0.9226{z^{ - 2}}

 

Con el minimo común multiplo, pasaría a calcular los nuevos polinomios “B” los cuales incluiran los retardos efectivos.

 

\left[ {\begin{array}{*{20}{c}}  {{y_1}(t)}\\  {{y_2}(t)}  \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}  1&0\\  0&{{z^{ - 1}}}  \end{array}} \right]\left[ {\begin{array}{*{20}{c}}  {\dfrac{{0.3776\left( {1 - 0.9765{z^{ - 1}}} \right)}}{{{A_1}\left( {{z^{ - 1}}} \right)}}}&{\dfrac{{ - 0.4447\left( {1 - 0.9705{z^{ - 1}}} \right){z^{ - 2}}}}{{{A_1}\left( {{z^{ - 1}}} \right)}}}\\  {\dfrac{{0.2959\left( {1 - 0.9659{z^{ - 1}}} \right){z^{ - 1}}}}{{{A_2}\left( {{z^{ - 1}}} \right)}}}&{\dfrac{{ - 0.6621\left( {1 - 0.9552{z^{ - 1}}} \right)}}{{{A_2}\left( {{z^{ - 1}}} \right)}}}  \end{array}} \right]\left[ {\begin{array}{*{20}{c}}  {{u_1}(t - 1)}\\  {{u_2}(t - 1)}  \end{array}} \right]

 

Haciendo todas estas operaciones puedo recudir mi problema a la solución de varias ecuaciones MISO (Multi-Input Single-Output) osea multiples entradas y una sola salida, lo cual es mucho más sencillo de calcular.

 

{A_1}\left( {{z^{ - 1}}} \right)y\left( t \right) = [{B_{11}}\left( {{z^{ - 1}}} \right){u_1}\left( {t - 1} \right) + {B_{12}}\left( {{z^{ - 1}}} \right){u_2}\left( {t - 1} \right)]

 

{A_2}\left( {{z^{ - 1}}} \right){\rm{y}}\left( {\rm{t}} \right) = {z^{ - 1}}[{B_{21}}\left( {{z^{ - 1}}} \right){u_1}\left( {t - 1} \right) + {B_{22}}\left( {{z^{ - 1}}} \right){u_2}\left( {t - 1} \right)]

 

Por ultimo, podemos ver que las predicciones del sistema MIMO van a partir desde la formulación del modelo multivariable:

 

{A_i}\left( {{z^{ - 1}}} \right){y_i}\left( t \right) = {z^{ - {d_i}}}{{\bf{B}}_{\rm{i}}}\left( {{z^{ - 1}}} \right)u\left( {t - 1} \right) + \frac{1}{\Delta }e(t)

 

{{\bf{B}}_{\rm{i}}}\left( {{z^{ - 1}}} \right) = [{B_{i1}},{B_{i2}}, \ldots ,{B_{im}}]

 

Para poder implementar un GPC MIMO tenemos dos posibilidades equivalentes:

  1. Realizar la implementación de Manera Recursiva (Como lo vimos en el primer video de GPC SISO)
  2. Utilizando las Ecuaciones Diofantinas – Se repite el método aprendido, pero con cada uno de los modelos del sistema MIMO.

 

 

CALCULO DE LAS PREDICCIONES POR ECUACIONES DIOFANTINAS (Método 1)

Básicamente es similar al procedimiento hecho en el caso SISO. Para un modelo CARIMA de n salidas por m entradas puedo usar el siguiente polinomio identidad que me permite obtener la predicciones optimas de cada salida y(k+1):

 

1 = {E_{ik}}\left( {{z^{ - 1}}} \right)\Delta {A_i} + {z^{ - k}}{F_{ik}}\left( {{z^{ - 1}}} \right)

 

{\hat y_i}\left( {t + j} \right) = {F_{ik}}\left( {{z^{ - 1}}} \right){y_i}\left( t \right) + {z^{ - {d_i}}}{E_{ik}}\left( {{z^{ - 1}}} \right){{\bf{B}}_{\rm{i}}}\left( {{z^{ - 1}}} \right)\Delta u\left( {t - 1} \right)

 

La matriz A_i sera una matriz diagonal n\times n, el polinomio E_{ik} tendrá un orden de k-1 (numero de predicciones en el horizonte de predicción menos uno) y el polinomio F_{ik} tendrá un odern de n_{ai} es decir su tamaño depende del polinomio A relacionado a cada salida del sistema MIMO.

 

Retomando nuestro ejemplo con un horizonte de predicción de 4 y un horizonte de control igual a 2 Miremos las predicciones para la salida y1 comenzando desde el retardo mínimo d1

Ejemplo MIMO GPC

 

 

{\hat y_1}\left( {t + j + {d_1}} \right) = {F_{1k}}\left( {{z^{ - 1}}} \right){y_1}\left( t \right) + {z^{ - {d_1}}}{E_{1k}}\left( {{z^{ - 1}}} \right)\left[ {{{\bf{B}}_{11}}\left( {{z^{ - 1}}} \right)\Delta {u_1}\left( {t - 1} \right) + {{\bf{B}}_{12}}\left( {{z^{ - 1}}} \right)\Delta {u_2}\left( {t - 1} \right)} \right]

 

Donde las ecuaciones diofantinas E y F son:

Ecuacion Diofantina

 

{E_{1k}} = 1 + 2.9470{z^{ - 1}} + 5.7900{z^{ - 2}} + 9.4803{z^{ - 3}}

 

{F_{1k}} = \left[ {\begin{array}{*{20}{c}}  {\begin{array}{*{20}{c}}  {2.9470{z^{ - 1}} - 2.8946{z^{ - 2}} + 0.9477{z^{ - 3}}}\\  {5.7900{z^{ - 2}} - 7.5828{z^{ - 3}} + 2.7928{z^{ - 4}}}  \end{array}}\\  {\begin{array}{*{20}{c}}  {9.4803{z^{ - 3}} - 13.9673{z^{ - 4}} + 5.4870{z^{ - 5}}}\\  {13.9708{z^{ - 4}} - 21.9550{z^{ - 5}} + 8.9842{z^{ - 6}}}  \end{array}}  \end{array}} \right]

 

El producto entre el polinomio E y los polinomios B para la primera salida son:

E_1B_{11}=0.3776-0.3687z^{-1}
E_2B_{11}=0.3776+0.744z^{-1}-1.0865z^{-2}
E_3B_{11}=0.3776+0.744z^{-1}+1.0996z^{-2}-2.1346z^{-3}E_4B_{11}=0.3776+0.744z^{-1}+1.0996z^{-2}+1.4447z^{-3}-3.4951z^{-4}

E_1B_{12}=(-0.4447+0.4316z^{-1})z^{-2}
E_2B_{12}=(-0.4447-0.8789z^{-1}+1.2718z^{-2})z^{-2}
E_3B_{12}=(-0.4447-0.8789z^{-1}-1.3029z^{-2}+2.4988z^{-3})z^{-2}
E_4B_{12}=(-0.4447-0.8789z^{-1}-1.3029z^{-2}-1.7169z^{-3}+4.0914z^{-4})z^{-2}

 

Para entender bien las predicciones voy a separar cada entrada en matrices diferentes, para ordenar los términos en las predicciones de la parte forzada y la parte libre.

 

Controle preditivo MIMO GPC

Controle preditivo MIMO GPC

 

 

CALCULO DE LAS PREDICCIONES DE FORMA RECURSIVA (Método 2)

 

{A_i}\left( {{z^{ - 1}}} \right){y_i}\left( t \right) = {z^{ - {d_i}}}{{\bf{B}}_{\rm{i}}}\left( {{z^{ - 1}}} \right)u\left( {t - 1} \right) + \frac{1}{\Delta }e(t)

 

{A_1}\left( {{z^{ - 1}}} \right){y_1}\left( t \right) = {z^{ - {d_1}}}{{\bf{B}}_1}\left( {{z^{ - 1}}} \right)u\left( {t - 1} \right) + \frac{1}{\Delta }e(t)

 

\Delta {A_1}\left( {{z^{ - 1}}} \right){y_1}\left( t \right) = {z^{ - {d_1}}}{{\bf{B}}_1}\left( {{z^{ - 1}}} \right)\Delta u\left( {t - 1} \right) + e(t)

 

\left( {1 - {z^{ - 1}}} \right)\left( {1 - 1.9470{z^{ - 1}} + 0.9477{z^{ - 2}}} \right){y_1}\left( t \right) = \left( {0.3776\left( {1 - 0.9765{z^{ - 1}}} \right)} \right)\Delta {u_1}\left( {t - 1} \right) + ( - 0.4447{z^{ - 2}}\left( {1 - 0.9705{z^{ - 1}}} \right))\Delta {u_2}\left( {t - 1} \right)

 

\left( {1 - 2.9470{z^{ - 1}} + 2.8946{z^{ - 2}} - 0.9477{z^{ - 3}}} \right){y_1}\left( t \right) = \left( {0.3776 - 0.3687{z^{ - 1}}} \right)\Delta {u_1}\left( {t - 1} \right) + ( - 0.4447{z^{ - 2}} + 0.4316{z^{ - 3}})\Delta {u_2}\left( {t - 1} \right)

 

De esa forma llegamos a que nuestra ecuación de salida para la variable 1 (y1) viene dado por la siguiente ecuación:

 

{y_1}\left( t \right) = 2.9470{y_1}\left( {t - 1} \right) - 2.8946{y_1}\left( {t - 2} \right) + 0.9477{y_1}\left( {t - 3} \right) + 0.3776\Delta {u_1}\left( {t - 1} \right) - 0.3687\Delta {u_1}\left( {t - 2} \right) - 0.4447\Delta {u_2}\left( {t - 3} \right) + 0.4316\Delta {u_2}\left( {t - 4} \right)

 

Y que podriamos hacer todas las predicciones hasta el horizonte de predicción N=4, note que de forma recursiva siempre en una predicción debo usar el valor que calculé en la predicción anterior, por ejemplo en la predicción de (y1(t+2)) tengo que usar el valor numérico que obtuve en la predicción anterior osea en (y1(t+1)).

 

{{\hat y}_1}\left( {t + 1} \right) = 2.9470{y_1}\left( t \right) - 2.8946{y_1}\left( {t - 1} \right) + 0.9477{y_1}\left( {t - 2} \right) + 0.3776\Delta {u_1}\left( t \right) - 0.3687\Delta {u_1}\left( {t - 1} \right) - 0.4447\Delta {u_2}\left( {t - 2} \right) + 0.4316\Delta {u_2}\left( {t - 3} \right)

 

{{\hat y}_1}\left( {t + 2} \right) = 2.9470{y_1}\left( {t + 1} \right) - 2.8946{y_1}\left( t \right) + 0.9477{y_1}\left( {t - 1} \right) + 0.3776\Delta {u_1}\left( {t + 1} \right) - 0.3687\Delta {u_1}\left( t \right) - 0.4447\Delta {u_2}\left( {t - 1} \right) + 0.4316\Delta {u_2}\left( {t - 2} \right)

 

{{\hat y}_1}\left( {t + 3} \right) = 2.9470{y_1}\left( {t + 2} \right) - 2.8946{y_1}\left( {t + 1} \right) + 0.9477{y_1}\left( t \right) + 0.3776\Delta {u_1}\left( {t + 2} \right) - 0.3687\Delta {u_1}\left( {t + 1} \right) - 0.4447\Delta {u_2}\left( t \right) + 0.4316\Delta {u_2}\left( {t - 1} \right)

 

{{\hat y}_1}\left( {t + 4} \right) = 2.9470{y_1}\left( {t + 3} \right) - 2.8946{y_1}\left( {t + 2} \right) + 0.9477{y_1}\left( t+1 \right) + 0.3776\Delta {u_1}\left( {t + 3} \right) - 0.3687\Delta {u_1}\left( {t + 2} \right) - 0.4447\Delta {u_2}\left( t+1 \right) + 0.4316\Delta {u_2}\left( {t } \right)

 

Podemos hacer el mismo procedimiento para la salida y2 y calcular todas las predicciones que queramos

 

{y_2}\left( t \right) = 2.921{y_2}\left( {t - 1} \right) - 2.8436{y_2}\left( {t - 2} \right) + 0.9226{y_2}\left( {t - 3} \right) + 0.2959\Delta {u_1}\left( {t - 3} \right) - 0.2858\Delta {u_1}\left( {t - 4} \right) - 0.6621\Delta {u_2}\left( {t - 2} \right) + 0.6324\Delta {u_2}\left( {t - 3} \right)

 

Vemos que como en esta segunda salida tenemos un atraso minimo de 1, nuestra predicción comienza en d+1=2 hasta nuestro horizonte de predicción.

 

{y_2}\left( t+2 \right) = 2.921{y_2}\left( {t + 1} \right) - 2.8436{y_2}\left( {t } \right) + 0.9226{y_2}\left( {t - 1} \right) + 0.2959\Delta {u_1}\left( {t - 1} \right) - 0.2858\Delta {u_1}\left( {t - 2} \right) - 0.6621\Delta {u_2}\left( {t } \right) + 0.6324\Delta {u_2}\left( {t - 1} \right)

 

{y_2}\left( t+3 \right) = 2.921{y_2}\left( {t + 2} \right) - 2.8436{y_2}\left( {t+1 } \right) + 0.9226{y_2}\left( {t } \right) + 0.2959\Delta {u_1}\left( {t } \right) - 0.2858\Delta {u_1}\left( {t - 1} \right) - 0.6621\Delta {u_2}\left( {t+1 } \right) + 0.6324\Delta {u_2}\left( {t} \right)

 

{y_2}\left( t+4 \right) = 2.921{y_2}\left( {t + 3} \right) - 2.8436{y_2}\left( {t+2 } \right) + 0.9226{y_2}\left( {t +1} \right) + 0.2959\Delta {u_1}\left( {t+1 } \right) - 0.2858\Delta {u_1}\left( {t } \right) - 0.6621\Delta {u_2}\left( {t+2 } \right) + 0.6324\Delta {u_2}\left( {t+1} \right)

 

{y_2}\left( t+5 \right) = 2.921{y_2}\left( {t + 4} \right) - 2.8436{y_2}\left( {t+3 } \right) + 0.9226{y_2}\left( {t +2} \right) + 0.2959\Delta {u_1}\left( {t+2 } \right) - 0.2858\Delta {u_1}\left( {t+1 } \right) - 0.6621\Delta {u_2}\left( {t+3 } \right) + 0.6324\Delta {u_2}\left( {t+2} \right)

 

Ley de Control

Sabemos que el Calculo de la ley de control para el caso irrestricto viene dado por:

\Delta U=(G^TQ_{\delta}G+Q_{\lambda})^{-1}G^TQ_{\delta}(W-f)

donde W es la referencia y f es la respuesta libre.

 

Sabíamos que el primer termino de esta ley de control basta con calcularlo una única vez, o sea que podríamos decir que:

K=(G^TQ_{\delta}G+Q_{\lambda})^{-1}G^TQ_{\delta}

llegando a la ley de control:

\Delta U=K(W-f)

 

Y que de la matriz K solo tomabamos la primera fila, por causa del horizonte deslizante. Para un proceso Multivariable, vamos a tomar también la primera fila de cada bloque asociado a cada horizontes de control, asi:

GPC MIMO

Asi observamos que debemos tomar toda la primera fila completa de K de cada horizonte de control, para poder formar nuestra ley de control \Delta U. En la imagen, habían 3 incrementos de control, por lo tanto se tomaran 3 K diferentes para poder formar las 3 leyes de control.

\begin{bmatrix}  \Delta u_1\\  \Delta u_2\\  \vdots \\  \Delta u_{N_u}  \end{bmatrix}=\begin{bmatrix}  k_1\\  k_2\\  \vdots \\  k_{N_u}  \end{bmatrix}\begin{bmatrix}  W-fr_1\\  W-fr_2\\  \vdots \\  W-fr_{N_u}  \end{bmatrix}

 

Recordando que cada k es un vector fila.

(Agradecimientos a Deinis Muñoz Muñoz por contribuir con las imagenes y formulas de la ley de control)

 

Retomando nuestro ejemplo con N=4 y Nu=2, la matriz K seria dado por:

K=\begin{bmatrix}  0.0270& 0.0462& 0.0637& 0.0799& -0.0018& 0.0165& 0.0285& 0.0400\\  -0.0070& 0.0190& 0.0435 & 0.0669 & -0.0011 & -0.0081& 0.0107 & 0.0287\\  0.0010& 0.0026 & 0.0002 & -0.0020 & -0.0060 & -0.0106 & -0.0145 & -0.0184\\  0.0006 & 0.0015 & 0.0027 & -0.0004 & 0.0004 & -0.0050 & -0.0100 & -0.0149   \end{bmatrix}

 

Donde la matriz anterior podemos obtener las dos ganancias para el calculo de los incrementos de control. La ganancia 1 (K_1) seria la primera fila de la matriz. La ganancia 2 seria tomar la primera fila donde comienza el bloque para la segunda salida, y para poder descubrir cual fila corresponde debemos observar el horizonte de control anterior (de la ganancia 1 o primera salida). Sabemos que el horizonte de control de ambas salidas es Nu_1=2  y Nu_2=2, por lo tanto la ganancia esta en la fila: Nu_1+1=3 . Entonces debo tomar toda la fila 3 para la ganancia 2 (K_2 )

 

K_1=\begin{bmatrix}0.0270& 0.0462& 0.0637& 0.0799& -0.0018& 0.0165& 0.0285& 0.0400\end{bmatrix}

 

K_2=\begin{bmatrix}0.0010& 0.0026 & 0.0002 & -0.0020 & -0.0060 & -0.0106 & -0.0145 & -0.0184\end{bmatrix}

Ejemplo

Considere la siguiente columna de destilación la cual se muestra en la siguiente figura. (Este ejemplo es derivado de el modelamiento de Wood and Berry)

Columna de destilación

A continuación se muestra la función de transferencia que representa el modelo, se pide realizar con control MIMO GPC. (Discretice el modelo con Ts=0.5). Este ejemplo es el que hemos venido realizando en toda esta entrada, con todos los cálculos que ya vimos en la parte superior.

 

P\left( s \right) = \left[ {\begin{array}{*{20}{c}}  {\dfrac{{12.8}}{{16.7s + 1}}}&{\dfrac{{ - 18.9{e^{ - 1s}}}}{{21s + 1}}}\\  {\dfrac{{6.6{e^{ - 1s}}}}{{10.9s + 1}}}&{\dfrac{{ - 19.4{e^{ - 0.5s}}}}{{14.4s + 1}}}  \end{array}} \right]

Si te está sirviendo de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Para Brasil

Se você está em Brasil pode utilizar o botão que está aqui em baixo, para realizar o convite em Reais.




A continuación se muestra el código de implementación en MATLAB, recuerda que para ver el código debes compartir el contenido de este post por redes sociales, para que más personas puedan beneficiarse de esta información.

Ejemplo Forma Diofantina

Ejemplo Forma Recursiva

 

Para que los códigos funciones, debes bajar las siguientes librerias y colocarlas en la carpeta de MATLAB de “mis documentos” o simplemente colocarlas en la carpeta donde estas trabajando tu código MIMO GPC.

 

>> DESCARGAR CODIGOS (CLICK ACÁ) <<

 

GPC MIMO POR DIOFANTINAS:

MIMO RECURSIVO


Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Para Brasil

Se você está em Brasil pode utilizar o botão que está aqui em baixo, para realizar o convite em Reais.




Esta entrada tiene 13 comentarios

  1. Hola Sergio, quisiera preguntarte si hay algun metodo para implementar un control MPC con un plc clasico s7-300 o s7-400 de siemens. Gracias

    1. Hola Carlos, yo no trabajo mucho con PLCs, pero talvez puedas programarlo tu mismo. Yo se que para el s7-400 existe un complemento APC( Control Avanzado de Procesos) para PCS7, allí tienes una función de MPC. Saludos.

  2. Hola Sergio, el sistema SIMO del que te consultaba lo convertí en un sistema MIMO, ya que las salidas dependen también de otras variables que no fueron consideradas.
    El control MPC MIMO puede controlar cualquier sistema?, entiendo mas o menos la teoría de control MIMO, pero a pesar de que sintonicé las variables de control no puedo controlar el sistema:
    Ps=[tf(0.0083,[1 0.0258]) tf(0.000317,[1 0]);tf(0.00367,[1 0.2611 0.00818]) tf(0.000317,[1 0])];
    es un sistema de control de temparatura y humedad, sin retardos, que después de sintonizarlo:
    Ts=500
    N=[140;140];
    N1=[dmin(1)+1 dmin(2)+1];
    N2=[dmin(1)+N(1) dmin(2)+N(2)];
    Nu=[120;120];
    lambda=[10 10];
    delta=[80 80];
    se puede controlar la salida pero la señal de control “u1” siempre sale negativa. Agradecería tu comentario.

    1. La estrategia de MPC es una estrategia avanzada de control que es utilizada en procesos bastante complejos como es el caso de plataformas petroleras o industrias de refinería. En tu simulación el MPC consigue controlar tranquilamente ese sistema que expones, claro, las señales de control te dan negativo es por causa del modelo lineal que usas en la salida 2. Tienes 2 integradores Puros. Puedes interpretarlo como 2 capacitores que almacenan energia pero no tienen quien la consuma, o puedes verlo como si fueran tanques que se estan llenando de agua, pero no tienen abierta la valvula de salida, por lo tanto solo almacena energia. De esa forma el MPC calcula que debe inyectar señales negativas (disipar energia) para poder alcanzar el setpoint. Muy seguramente la implemetación real, no va a ser integrador puro pues existen muchos factores que te van a disipar la energia. Si no quieres que te salgan las señales de control negativas, puedes modelar el sistema por una FT de primer orden, o activar las restricciones de señal de control, usando el quadprog. Saludos.

    2. También tienes que tener en cuenta que aqui estamos simulando. En tu proceso real, ni tu temperatura ni tu humedad van a comenzar desde el cero absoluto, por ejemplo la temperatura puede comenzar en 25 y la humedad en 70. Puedes simular esto también en el código inicializando los vectores de salida y y los vectores de entrada u. Si ves todos esos vectores estan inicializados en cero, en tu caso puedes inicializarlos desde la condición estacionaria de tu proceso.

  3. Hola Sergio, gracias por las orientaciones, en la simulación del ejemplo mediante ecuaciones diofantinas, laley de control en una parte se hace negativa, como se puede explicar eso? o a que se debe que sea negativa?

    1. Hola Juan Carlos, el ejemplo se encuentra normalizado entre 1 y -1. Porque el objetivo es ver como funciona el MIMO MPC. la normalización de variables es un factor fundamental en procesos multivariables donde se trabajan con diferentes unidades. Si se desea el ejemplo podría trabajarse sobre los puntos de operación del proceso que sería de y1=96.25mol%, y2=0.5mol%, u1 =1.95lb/min, u2=1.71 lb/min, asi no tendrías los valores negativos. Saludos.

  4. Hola Sergio, primero agradecerte por toda la información, y consultarte como cambio todo para un sistema de una entrada y dos salidas

    1. Hola Juan Carlos, los códigos que yo coloco en la pagina NO son genéricos y yo lo hago justamente asi para evitar que las personas simplemente copien y peguen el código y no intenten entender como funciona el MPC. Para tu caso especifico creo que podrías cerar las dos funciones de transferencia de la derecha, te quedaría un sistema de dos salidas y una entrada:
      Ps=[tf(12.8,[16.7 1]) 0;tf(6.6,[10.9 1]) 0];

      ya en la sintonia, puedes ponderar cual salida deseas darle más prioridad. Por ejemplo si quiero que la salida 1 tenga 80% de prioridad y la salida 2 tenga un 20%, colocaría el Delta asi:

      delta=[0.8 0.2];

      Éxitos.

  5. Hola Sergio, Primero agradecerte por la información que publicas, es muy interesante y me ayudo a hacer mi tesis. quisiera que me ayudes en lo siguiente: tengo un proceso SIMO, un entrada y dos salidas, como puedo acondicionarlo al control MIMO que hiciste.

    1. En ese caso tu sistema de FT se reduce a un vector columna. Deberás calcular todo desde esa entrada para tus dos salidas. Como tienes menos grados de libertad tendras que priorizar cual salida es más importante para ti, o incluso trabajar por zonas, que es lo que se usa en este caso, es decir que no tendras setpoint fijos si no que tendras un setpoint por intervalos con un limite máximo y mínimo y hacer que el algoritmo de optimización encuentre el mejor camino posible. No se que tan complejo sea tu sistema, tal vez si es sencillo puede incluso lograr controlar ambas variables en setponts fijos. Estos temas todavía no los he tratado aqui en la página, por eso es recomendable que leas articulos o libros para que lo logres. Igual aqui en el blog tienes la base de como funciona el controlador y un código NO genérico base también.

  6. Usted ha hecho un trabajo muy bueno.

    Parabens.

    ¿Dónde está el código de GPC2Metodos.m?

    1. Obrigado. No final da pagina estão todos os códigos e o link de descarga também. Abs.

Deja un comentario

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.

MIMO MPC
Menú de cierre