Saltar al contenido

En esta entrada veremos un poco el funcionamiento de un controlador predictivo basado en modelo (MPC) en su resolución GPC (Generalized Predictive Control) aplicado en procesos multivariables o por sus siglas en ingles MIMO (Multi-Input Multi-Output).

Antes de comenzar, si te interesa aprender más sobre la teoría del control y en especial sobre control avanzado, te invito a que mires nuestro curso gratuito de Control Predictivo y que te suscribas al canal de YouTube.

Sistema Multivariable – MIMO

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.

O sea 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) o sea 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 o sea 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]

[magicactionbox id=»1184″]

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

[sociallocker id=948]

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 GPC for Wood an Berry column
% Sergio Andres Castaño Giraldo
% https://controlautomaticoeducacion.com/mimo-mpc/
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - 2017
% Universidade Federal de Santa Catarina
% Florianópolis - 2016
% Por Ecuaciones Diofantinas
%_____________________________________________________________________

clc
close all
clear all

%% Define o processo do reator:

Ps=[tf(12.8,[16.7 1]) tf(-18.9,[21 1]);tf(6.6,[10.9 1]) tf(-19.4,[14.4 1])];
%Ps.iodelay=[1 3;7 3];
Ps.iodelay=[0 1;1 0.5];

Ts=0.5;                     %Tiempo de Muestreo
Pz=c2d(Ps,Ts);              %Processo discreto

[Bp,Ap,dp]=descompMPC(Pz);     %Descompone num, den e atraso de Pz em celdas

%% Parametros de sintonia del GPC
%Encuentro cual es el retardo minimo en mi función de transferencia
dmin=[100 100]';
for i=1:2
    dmin(i)=min(dp(i,:));
end

N=[4;4];                              %Ventana de Predicción
N1=[dmin(1)+1 dmin(2)+1];               %Horizonte Inicial
N2=[dmin(1)+N(1) dmin(2)+N(2)];         %Horizonte final
Nu=[2;2];                               %Horizonte de control
lambda=[10 100];                        %Ponderación de la acción de control
delta=[1 1];                            %Ponderación del seguimiento de referencia

%Matrices en bloque de los parametros de poderación
Ql=blkdiag(lambda(1)*eye(Nu(1)),lambda(2)*eye(Nu(2)));  %(lambda -> Acción de Control)
Qd=blkdiag(delta(1)*eye(N(1)),delta(2)*eye(N(2)));      %(delta  -> Seguimiento de Referencia)

%% Calculo de la ecuacion Diofantina
%Se obtiene la matriz B y A (Numerador y denominador sistema MIMO)
% la matriz A es el minimo común denominador del sistema
[B,A] = BA_MIMO(Bp,Ap);
%Función que Calcula Los Polinomios de la ecuación Diofantina
[E,En,F] = diophantineMIMO(A,N,dmin);

%% Matriz de la Respuesta Forzada G
[G] = MatrizG(E,B,N,Nu,dp);

%Sistema sin Restricciones
M=inv(G'*Qd*G+Ql)*G'*Qd;    %Ganancia M
K1=M(1,:);                  %Tomo la primera columna unicamente (para control u1)
K2=M(Nu(1)+1,:);            %Tomo valores de M desde [N(1)+1] una muestra despues del primer 
                            %horizonte de control (para control u2)

%Determino el Vector con los controles pasados utilizando la funcion:
uG = deltaUFree(B,En,N,dp);

%% Loop de Controle
% inicializa parametros de Simulacion
nit=400;                        %Numero de interacciones
deltaU1(1:nit) = 0;             % Inicializa controle incremental
deltaU2(1:nit) = 0;

% Señales del processo MIMO 2x2 (Inicialización)
%Salidas de cada bloque MIMO
yp11(1:nit) = 0; yp12(1:nit) = 0; yp21(1:nit) = 0; yp22(1:nit) = 0; 
%Salidas MIMO
y1(1:nit) = 0; y2(1:nit) = 0; 
%Referencias
r1(40:nit) = 0.8; r2(150:nit) = 0.4; 
%Acciones de Control
u1(1:nit) = 0; u2(1:nit) = 0; 

% Determino el tamaño de los polinomios F
lduf1=size(F{1}); %Averiguo la longitud del polinomio F
lduf2=size(F{2}); %Averiguo la longitud del polinomio F
lduf1=lduf1(1,2);
lduf2=lduf2(1,2);
F1=F{1};
F2=F{2};

%_________ Vectores de Control Pasado _____________________________
%Determino el numero de columnas de los vectores pasados
luG11=size(uG{1,1}); 
luG12=size(uG{1,2}); 
luG21=size(uG{2,1}); 
luG22=size(uG{2,2}); 

luG11=luG11(1,2);
luG12=luG12(1,2);
luG21=luG21(1,2);
luG22=luG22(1,2);

duf11=zeros(1,luG11); %Vector incrementos de control pasados libre
duf12=zeros(1,luG12); %Vector incrementos de control pasados libre
duf21=zeros(1,luG21); %Vector incrementos de control pasados libre
duf22=zeros(1,luG22); %Vector incrementos de control pasados libre

%Comienzo el lazo de Control
for k=40:nit
     % Salida del proceso Real MIMO
     yp11(k) = transferencia2(Bp{1,1},Ap{1,1},dp(1,1),u1,yp11,k);
     yp12(k) = transferencia2(Bp{1,2},Ap{1,2},dp(1,2),u2,yp12,k);
     yp21(k) = transferencia2(Bp{2,1},Ap{2,1},dp(2,1),u1,yp21,k);
     yp22(k) = transferencia2(Bp{2,2},Ap{2,2},dp(2,2),u2,yp22,k);

     y1(k)=yp11(k)+yp12(k);
     y2(k)=yp22(k)+yp21(k);

      %Puedes simular tambien el proceso de la siguiente forma:
      %t = 0:Ts:(k-1)*Ts;
      %yest=lsim(Ps,[u1(1:k);u2(1:k)],t,'zoh');
      %y1(k)=yest(end,1);
      %y2(k)=yest(end,2);
        
     %Calculo de la respuesta libre
     free1=0;
     free2=0;
     for i=1:lduf1
         free1=free1+y1(k-i+1)*F1(:,i);
     end
     for i=1:lduf2
         free2=free2+y2(k-i+1)*F2(:,i);
     end
     
     %Adiciono vectores de controles pasados
     free1=free1+uG{1,1}*duf11'+uG{1,2}*duf12';
     free2=free2+uG{2,1}*duf21'+uG{2,2}*duf22';
        
     %Calculo del incremento de Control
     %Proyeto sin referencias futuras
     deltaU1(k)=K1*[(r1(k)*ones(1,length(free1(1,1:end)))'-free1);(r2(k)*ones(1,length(free2(1,1:end)))'-free2)];
     deltaU2(k)=K2*[(r1(k)*ones(1,length(free1(1,1:end)))'-free1);(r2(k)*ones(1,length(free2(1,1:end)))'-free2)];
     
     %Calculo de la ley de control
     if k==1
        u1(k)=deltaU1(k);
        u2(k)=deltaU2(k);
     else
        u1(k)=u1(k-1)+ deltaU1(k);
        u2(k)=u2(k-1)+ deltaU2(k);
     end
     
    %Actualizo vectores de Controles Pasados 
    aux_2=duf11(1:end-1);
    duf11=[deltaU1(k) aux_2];
    
    aux_2=duf12(1:end-1);
    duf12=[deltaU2(k) aux_2];
    
    aux_2=duf21(1:end-1);
    duf21=[deltaU1(k) aux_2];
    
    aux_2=duf22(1:end-1);
    duf22=[deltaU2(k) aux_2];

end

%Grafico Resultados
nm=nit;
t = 0:Ts:(nm-1)*Ts;
figure
subplot(2,1,1),plot(t,r1,t,r2,t,y1,t,y2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
grid on;
hold
subplot(2,1,2),plot(t,u1,t,u2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('u')
grid on;

MIMO RECURSIVO

%% MIMO GPC for Wood an Berry column
% Sergio Andres Castaño Giraldo
% https://controlautomaticoeducacion.com/mimo-mpc/
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - 2017
% Universidade Federal de Santa Catarina
% Florianópolis - 2016
% Por forma Recursiva
%_____________________________________________________________________
clc
close all
clear all
%% Define o processo do reator:
Ps=[tf(12.8,[16.7 1]) tf(-18.9,[21 1]);tf(6.6,[10.9 1]) tf(-19.4,[14.4 1])];
%Ps.iodelay=[1 3;7 3];
Ps.iodelay=[0 1;1 0.5];
Ts=0.5;                     %Tiempo de Muestreo
Pz=c2d(Ps,Ts);              %Processo discreto
[Bp,Ap,dp]=descompMPC(Pz);     %Descompone num, den e atraso de Pz em celdas
%% Parametros de sintonia del GPC
%Encuentro cual es el retardo minimo en mi función de transferencia
dmin=[100 100]';
for i=1:2
dmin(i)=min(dp(i,:));
end
N=[10;10];                              %Ventana de Predicción
N1=[dmin(1)+1 dmin(2)+1];               %Horizonte Inicial
N2=[dmin(1)+N(1) dmin(2)+N(2)];         %Horizonte final
Nu=[2;2];                               %Horizonte de control
lambda=[10 100];                        %Ponderación de la acción de control
delta=[1 1];                            %Ponderación del seguimiento de referencia
%Matrices en bloque de los parametros de poderación
Ql=blkdiag(lambda(1)*eye(Nu(1)),lambda(2)*eye(Nu(2)));  %(lambda -> Acción de Control)
Qd=blkdiag(delta(1)*eye(N(1)),delta(2)*eye(N(2)));      %(delta  -> Seguimiento de Referencia)
%% Calculo de la ecuacion Diofantina
%Se obtiene la matriz B y A (Numerador y denominador sistema MIMO)
% la matriz A es el minimo común denominador del sistema
[Bt,A] = BA_MIMO(Bp,Ap);
%Calculo de las Diofantinas
[E,En,F] = diophantineMIMO(A,N,dmin);
%% Matriz de la Respuesta Forzada G
[G] = MatrizG(E,Bt,N,Nu,dp);
%Sistema sin Restricciones
M=inv(G'*Qd*G+Ql)*G'*Qd;    %Ganancia M
K1=M(1,:);                  %Tomo la primera columna unicamente (para control u1)
K2=M(Nu(1)+1,:);            %Tomo valores de M desde [N(1)+1] una muestra despues del primer 
%horizonte de control (para control u2)
%% Loop de Controle
%% inicializa parametros de Simulacion
nit=400; %Numero de interacciones
%Inicializa vector de las respuestas libres
y_free1=[];
y_free2=[];
%Inicializa vectores de los incrementos de control
deltaU1(1:nit) = 0;            
deltaU2(1:nit) = 0;
% Señales del processo MIMO 2x2 (Inicialización)
%Salidas de cada bloque MIMO
yp11(1:nit) = 0; yp12(1:nit) = 0; yp21(1:nit) = 0; yp22(1:nit) = 0; 
%Salidas MIMO
y1(1:nit) = 0; y2(1:nit) = 0; 
%Referencias
r1(40:nit) = 0.8; r2(150:nit) = 0.4; 
%Acciones de Control
u1(1:nit) = 0; u2(1:nit) = 0; 
%Atil (Producto entre polinomio A y el integrador)
At{1,1}=conv(A{1,1},[1 -1]);
At{2,2}=conv(A{2,2},[1 -1]);
%Determino el tamaño de cada polinomio Atil
na1=length(At{1,1})-1;                       
na2=length(At{2,2})-1;   
%Determino el tamaño de cada polinomio B proveniente del calculo del minimo
%comun denominador de la funcion "BA_MIMO"
nb1_11=length(Bt{1,1});                       
nb1_21=length(Bt{2,1});                       
nb1_12=length(Bt{1,2});                       
nb1_22=length(Bt{2,2});                       
%Comienzo el Lazo de Control
for k=40:nit
% Salida del proceso Real
yp11(k) = transferencia2(Bp{1,1},Ap{1,1},dp(1,1),u1,yp11,k);
yp12(k) = transferencia2(Bp{1,2},Ap{1,2},dp(1,2),u2,yp12,k);
yp21(k) = transferencia2(Bp{2,1},Ap{2,1},dp(2,1),u1,yp21,k);
yp22(k) = transferencia2(Bp{2,2},Ap{2,2},dp(2,2),u2,yp22,k);
y1(k)=yp11(k)+yp12(k);
y2(k)=yp22(k)+yp21(k);
%Puedes simular tambien el proceso de la siguiente forma: 
%t = 0:Ts:(k-1)*Ts;
%yest=lsim(Ps,[u1(1:k);u2(1:k)],t,'zoh');
%y1(k)=yest(end,1);
%y2(k)=yest(end,2);
%% Creo los vectores que conforman la RESPUESTA LIBRE
%Creo el vector de controles pasados basandome en la ecuacion de
%salida del modelo de la planta actual y(k). O sea coloco en un vector
%todos los incrementos de control de cada una de las entradas, algo como:
%[deltaU(k-1) deltaU(k-2) ... deltaU(k-nb)] donde "nb" es
%la longitud del vector B (Numerador del modelo)
%Controles pasados de y1 en el instante actual
deltaU11_pas=deltaU1(k-dp(1,1)-1:-1:k-dp(1,1)-nb1_11);           
deltaU12_pas=deltaU2(k-dp(1,2)-1:-1:k-dp(1,2)-nb1_12); 
%Controles pasados de y2 en el instante actual
deltaU21_pas=deltaU1(k-dp(2,1)-1:-1:k-dp(2,1)-nb1_21);           
deltaU22_pas=deltaU2(k-dp(2,2)-1:-1:k-dp(2,2)-nb1_22); 
%Creo un vector de salidas pasadas, similar al vector de controles
%pasados solo que en este caso con las salidas hasta la longitud del
%polinomio A (denominador del modelo)
y_pas1=y1(k:-1:k-na1+1);                         % Vetor de saída
y_pas2=y2(k:-1:k-na2+1);                         % Vetor de saída
%Ahora realizo mis predicciones hasta el horizonte de predicción N2
n=1; %Indice para la salida libre
for j=1:N2(1)
if j<=dp(1,1) 
%Comienzo a desplazar los deltaU pasados a medida que aumento
%la ventana de predicción, hasta que la predicción en este
%vector llegue al retardo
deltaU11_pas=[deltaU1(k-dp(1,1)-1+j) deltaU11_pas(1:end-1)];
elseif j>dp(1,1) 
%Una vez la ventana de predicción es mayor que el retardo, en
%la ecuación de salida dejan de existir los deltaU pasados y
%comenzaran a aparecer los deltaU futuros, por eso comienzo a
%cerar los deltaU pasados.
deltaU11_pas=[0 deltaU11_pas(1:end-1)];
end
if j<=dp(1,2) 
deltaU12_pas=[deltaU2(k-dp(1,2)-1+j) deltaU12_pas(1:end-1)];
elseif j>dp(1,2) 
deltaU12_pas=[0 deltaU12_pas(1:end-1)];
end
%Hago el calculo de la respuesta libre con el polinomio Atil y el
%polinomio B, junto con las salidas pasadas y los deltaU pasados
y_free1(n)=-y_pas1*At{1,1}(2:end)'+deltaU11_pas*Bt{1,1}(1:end)'...
+deltaU12_pas*Bt{1,2}(1:end)';% Armazena yfr
%Actualizo el vector de salidas pasadas
y_pas1=[y_free1(n) y_pas1(1:end-1)];                    % Atualiza yfr
n=n+1;
end
%Ahora realizo mis predicciones hasta el horizonte de predicción
n=1; %Indice para la salida libre
for j=1:N2(2)
if j<=dp(2,1)
deltaU21_pas=[deltaU1(k-dp(2,1)-1+j) deltaU21_pas(1:end-1)];
elseif j>dp(2,1)
deltaU21_pas=[0 deltaU21_pas(1:end-1)];
end
if j<=dp(2,2)
deltaU22_pas=[deltaU2(k-dp(2,2)-1+j) deltaU22_pas(1:end-1)];
elseif j>dp(2,2)
deltaU22_pas=[0 deltaU22_pas(1:end-1)];
end
y_free2(n)=-y_pas2*At{2,2}(2:end)'+deltaU21_pas*Bt{2,1}(1:end)'...
+deltaU22_pas*Bt{2,2}(1:end)';% Armazena yfr
%Actualizo el vector de salidas pasadas
y_pas2=[y_free2(n) y_pas2(1:end-1)];                    % Atualiza yfr
n=n+1;
end
%% Calculo del Incremento de Control   
%Calculo los incrementos de control utilizando unicamente los datos de
%las respuestas libres que se encuentran dentro de los horizontes de
%predicción N1 y N2 y hago el producto con la ganancia calculada de la
%respuesta forzada
deltaU1(k)=K1*[(r1(k)*ones(1,length(y_free1(N1(1):N2(1))))'-y_free1(N1(1):N2(1))');(r2(k)*ones(1,length(y_free2(N1(2):N2(2))))'-y_free2(N1(2):N2(2))')];
deltaU2(k)=K2*[(r1(k)*ones(1,length(y_free1(N1(1):N2(1))))'-y_free1(N1(1):N2(1))');(r2(k)*ones(1,length(y_free2(N1(2):N2(2))))'-y_free2(N1(2):N2(2))')];
%% Calculo la acción de Control
if k==1
u1(k)=deltaU1(k);
u2(k)=deltaU2(k);
else
u1(k)=u1(k-1)+ deltaU1(k);
u2(k)=u2(k-1)+ deltaU2(k);
end
end
%% Grafico los Resultados
nm=nit;
t = 0:Ts:(nm-1)*Ts;
figure
subplot(2,1,1),plot(t,r1,t,r2,t,y1,t,y2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
grid on;
hold
subplot(2,1,2),plot(t,u1,t,u2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('u')
grid on;

[/sociallocker]

Entradas relacionadas

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

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

Comentarios (16)

Hola, me gustaria saber como hacer un sistema mimo 2×2 pero con restricciones.
Muchas gracias

Responder

Hola Sergio. Estoy implementando control MPC en una planta real, con el MPC generando los Set Points de las mallas de control. Los PID ya estan bien sintonizados y funcionando con buenos indices de desempeño. Cuando coloco para funcionar el sistema, la respuesta es un poco oscilante (+- 0.5 del valor del set point), si pongo solo los PID a funcionar la salida alcanza el estado estacionario sin oscilacion ninguna.
Mi pregunta es que serán los MPC que estan mal sintonizados?
Usted tendría como disponibilizar algún link en español o portugués sobre sintonía de MPC?

Responder

Hola Manuel, probablemente si sea la sintonía del MPC la que te este causando eso, esa sintonia no es trivial, y es una de las etapas complicadas del proyecto. Por ahora no tengo videos sobre el tema, pero ciertamente lo voy a hacer en el futuro.
• Usar N y Nu grandes da un mejor desempeño en la predicción y en la mayoría de los casos muestra un buen desempeño tanto dentro del horizonte de predicción cuanto fuera de este (a excepción del sistema inestable)
• Generalmente aumentar mucho el Nu no tiene mucha diferencia porque los pasos de control se vuelven similares
Pasos en el Control
• Un lambda grande implica en incrementos de control pequeños pero consecuentemente seguimientos de referencia mayor.

Responder

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

Responder

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.

Responder

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.

Responder

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.

Responder

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.

Responder

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?

Responder

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.

Responder

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

Responder

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.

Responder

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.

Responder

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.

Responder

Usted ha hecho un trabajo muy bueno.

Parabens.

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

Responder

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

Responder