▷ Control por Mínima Varianza

El control por mínima varianza (MVC, por sus siglas en ingles Minimum Variance Control o en portugues Variancia Mínima) tiene como objetivo reducir el efecto de las perturbaciones sobre la salida. Es de especial utilidad cuando la salida del proceso que se desea controlar se encuentra contaminada por una perturbación estocástica la cual no puede eliminarse completamente, pero si puede reducirse su varianza.

En la entrada de hoy abordaremos este interesante controlador y mostraremos algunos ejemplos de implementación en el software de MATLAB.

Contenidos mostrar

Control por Mínima Varianza – MVC

Como ya explicamos, un proceso cualquiera afectado por un ruido puede verse como:

Control por Mínima Varianza

La estrategia de este controlador consiste en aplicar la señal de control en el instante k o sea u(k) tomando los valores pasados de la entrada ( $u(k-1),u(k-2),...,u(k-n_b)$ ) de la salida( $y(k),y(k-1),y(k-2),...,y(k-n_a)$ ) de forma a minimizar un indice J

J=E\left(y^2(k+d+1/k)\right)\ [1]

Donde E es la esperanza matemática, k es el instante actual y d es el retardo.

Esta ecuación nos muestra que a cada periodo de muestreo k se debe determinar una señal de control que minimice este indice y la inclusión del retardo d+1 es una consideración razonable dado que cualquier acción de control apenas afectara la salida una vez el retardo haya transcurrido. Por eso para el desarrollo de este controlador vamos a suponer tres condiciones:

Todos los ceros del proceso deben estar dentro del circulo unitario (Sistemas de fase mínima)
El retardo del sistema debe ser conocido.
El orden del sistema debe ser conocido, o por lo menos establecido un limite máximo.

Con relación a la primera suposición, si aplicamos esta técnica a sistemas de fase no mínima (ceros fuera del circulo unitario), la estructura de control todavía será capaz de minimizar la varianza, pero sera sumamente sensible a cualquier variación de parámetros del sistema, lo que muestra que esta estrategia no es robusta.

A continuación se presenta el video en Youtube con la explicación que continua más abajo en el post. Al final del post está el código de implementación ejemplo junto con su video de Youtube.

Video en Español

Video en Portugues

Estructura de Control del MVC

Como ya dijimos que el proceso está siendo afectado por perturbaciones estocásticas y fuera de eso estamos trabajando en el dominio discreto, la ecuación que representa mi proceso físico lo representaremos por un modelo ARMAX descrito a continuación:

A(z^{-1})y(k)=z^{-d}B(z^{-1})u(k)+C(z^{-1})v(k)\ [2]

Donde los polinomios A, B y C vienen dados por:

A(z^{-1})=1+a_1z^{-1}+a_2z^{-2}+...+a_{n_a}z^{-n_a}

B(z^{-1})=b_1+b_2z^{-1}+b_3z^{-2}+...+b_{n_b}z^{-n_b}

C(z^{-1})=1+c_1z^{-1}+c_2z^{-2}+...+c_{n_c}z^{-n_c}

Lo que estamos viendo es que el polinomio C está filtrando la perturbación estocástica v(k) (Ruido Blanco) y es sumado con la respuesta del sistema. Dado que dicho polinomio aparece en el modelo del proceso, este termino es empleado como un parámetro de diseño del control por mínima varianza y en particular el grado de este polinomio es un grado menor al polinomio A, $n_c=n_a-1$ . Este polinomio C, esta fuertemente relacionado a un polinomio observador.

Deduciendo que la señal de control calculada u(k) va afectar el sistema una ves pasado el retardo, tal y como lo vimos en la entrada de sistemas con retardo del predictor de Smith, quiere decir que la señal u(k) afecta la salida en y(k+d+1). Esto quiere decir que para efectos de calcular la predición de la salida en el futuro, nos es más conveniente reescribir la ecuación [2] de la siguiente forma:

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+\dfrac{C(z^{-1})}{A(z^{-1})}v(k+d+1)\ [3]

De esta manera, la salida futura del sistema (que equivale a d intervalos de tiempo en el futuro, asumiendo k como el instante actual) es la suma de dos términos:

$\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)$ representa el efecto de la acción de control actual y todos los controles pasados en la planta
$\dfrac{C(z^{-1})}{A(z^{-1})}v(k+d+1)$ representa el futuro, presente y pasado de las perturbaciones.

Esta descomposición de la perturbación puede escribirse de manera explicita si escribimos el polinomio $\dfrac{C(z^{-1})}{A(z^{-1})}$ , utilizando la identidad polinomial $C=A(z^{-1})F(z^{-1})+z^{-(d+1)}G(z^{-1})$ , como:

\dfrac{C(z^{-1})}{A(z^{-1})}=F(z^{-1})+z^{-(d+1)}\dfrac{G(z^{-1})}{A(z^{-1})}\ [4]

Donde el grado de $G(z^{-1})$ es menor que el grado de $A(z^{-1})$ y

F(z^{-1})=1+f_1z^{-1}+f_2z^{-2}+...+f_{d}z^{-d}

G(z^{-1})=g_0+g_1z^{-1}+g_2z^{-2}+...+g_{n_a-1}z^{-n_a-1}

Sustituyendo la ecuación (4) en (3) tenemos que:

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+\left(F(z^{-1})+z^{-(d+1)}\dfrac{G(z^{-1})}{A(z^{-1})}\right) v(k+d+1)

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{A(z^{-1})}v(k)\ [5]

Los dos últimos términos del lado derecho de la ecuación 5 tienen el siguiente significado:

$F(z^{-1})v(k+d+1)$ es una combinación lineal de las perturbaciones producidas entre el instante $k$ y $k+d$ , donde el efecto sobre la salida del proceso $y(k+d)$ no se puede controlar con la ley de control actual $u(k)$ , debido a que la perturbación futura $v(k+d)$ para un retardo mayor a cero $d>0$ es totalmente independiente de los controles y salidas pasadas $u(k-1),u(k-2),...,u(k-n_b)$ , $y(k),y(k-1),y(k-2),...,y(k-n_a)$ .

$\dfrac{G(z^{-1})}{A(z^{-1})}v(k)$ es el efecto de las perturbaciones sobre la salida en instantes anteriores a k.

Si despejamos v(k) de la ecuación (2):

v(k)=\dfrac{A(z^{-1})}{C(z^{-1})}y(k)-z^{-d}\dfrac{B(z^{-1})}{C(z^{-1})}u(k)\ [6]

Y reemplazamos la ecuación (6) en (5)

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{A(z^{-1})}\left(\dfrac{A(z^{-1})}{C(z^{-1})}y(k)-z^{-d}\dfrac{B(z^{-1})}{C(z^{-1})}u(k)\right)

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{C(z^{-1})}y(k)-z^{-d}\dfrac{B(z^{-1})G(z^{-1})}{C(z^{-1})A(z^{-1})}u(k)\ [7]

de la ecuación (4) podemos inferir que:

z^{-(d+1)}\dfrac{G(z^{-1})}{A(z^{-1})}=\dfrac{C(z^{-1})}{A(z^{-1})}-F(z^{-1})

z^{-d}\dfrac{G(z^{-1})}{A(z^{-1})}=\dfrac{zC(z^{-1})}{A(z^{-1})}-zF(z^{-1})\ [8]

entonces reemplazamos (8) en (7)

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{C(z^{-1})}y(k)-\dfrac{B(z^{-1})}{C(z^{-1})}\left(\dfrac{zC(z^{-1})}{A(z^{-1})}-zF(z^{-1})\right) u(k)

y(k+d+1)=\dfrac{zB(z^{-1})}{A(z^{-1})}u(k)+F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{C(z^{-1})}y(k)-\dfrac{zB(z^{-1})}{A(z^{-1})} u(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)

Reduciendo tenemos:

y(k+d+1)=F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)\ [9]

La ecuación (9) representa la señal de salida de nuestro proceso que relaciona la entrada del sistema, el ruido y la salida del proceso, de la expresión observamos que los dos términos a la derecha corresponden a la predicción que se tiene de la salida real del proceso.

\hat{y}(k+d+1)=\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)

donde nuestra mejor predicción de salida viene dada por:

y(k+d+1)=F(z^{-1})v(k+d+1)+\hat{y}(k+d+1)

Ahora como lo que queremos es atenuar la varianza en la salida, procedemos a minimizar el criterio J y aplicando la esperanza matemática:

J=E\{[y(k+d+1)]^{2}\}

J=E\left\lbrace \left[ F(z^{-1})v(k+d+1)+\left(\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)\right) \right]^2\right\rbrace

Resolvemos el binomio:

J=E\{[F(z^{-1})v(k+d+1)]^{2}\}+E\left\lbrace \left[ \dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)\right]^2\right\rbrace\\\\+2E\left\lbrace \left(F(z^{-1})v(k+d+1)\right)\left( \dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)\right) \right\rbrace

Notemos que el ruido $v(k+d+1)$ es una señal en el futuro, como no sabemos el ruido que va a suceder en el futuro la mejor estimativa que tenemos es considerar que ese ruido será cero, es decir su esperanza matemática es CERO. Por lo tanto tenemos que:

J=E\{[y(k+d+1)]^{2}\}=E\left\lbrace \left[ \dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)\right]^2\right\rbrace

Es claro que la estimación óptima en el sentido de error cuadrático medio mínimo es dada por la igualación de la predicción de salida a cero:

\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)=0\ [10]

Si resolvemos la ecuación (10) para u(k) tenemos que:

u(k)=-\dfrac{G(z^{-1})}{zB(z^{-1})F(z^{-1})}y(k) [11]

Por lo tanto hemos llegado a una ecuación de la ley de control que minimiza la varianza de la salida del proceso. Conocido como control por mínima varianza, MVC, minimum variance control o variancia mínima. Con esta ecuación estamos garantizando que la varianza en la salida intente llegar a cero, sin embargo este es un controlador que presenta error en estado estable.

minimum variance control

Podemos observar también que este controlador de la ecuación (11) cancela los ceros del sistema en lazo abierto, dado que el polinomio B esta en el denominador. Si por alguna razón estos ceros se encuentran fuera del circulo unitario cualquier variación de los parámetros va a conducir en un comportamiento inestable del sistema en malla cerrada.

Este únicamente es un controlador regulatório (No sigue referencia), como en muchos casos es interesante seguir referencia, en lugar de aplicar la ley de control sobre la salida del proceso, ésta es aplicada sobre la diferencia entre el setpoint y la salida (Error), quedando la ley de control de la siguiente forma:

u(k)=\dfrac{G(z^{-1})}{zB(z^{-1})F(z^{-1})}[w(k)-y(k)]

Control por variancia mínima (portugues)

Sin embargo este controlador presenta error en estado estable, por lo tanto para eliminar el offset o error en estado estable, podemos incluir un integrador en el lazo de control asi:

u(k)=\dfrac{G(z^{-1})}{zB(z^{-1})F(z^{-1})}\left[ \dfrac{z-1+\alpha}{z-1} \right] [w(k)-y(k)]

donde $\alpha$ es un de diseño del integrador con valores entre 0 y 1.

Control por Mínima Varianza con Integrador

Control por Mínima Varianza con Seguimiento de Referencia

En este caso en lugar de minimizar la salida, se debe minimizar la relación que existe entre la salida del proceso y el setpoint (Error), de la siguiente forma:

J=E\{[y(k+d+1)-w(k+d+1)]^{2}\}\ [12]

Como ya sabemos el valor de $y(k+d+1)$ el cual es dado en la ecuación (9), procedemos a reemplazarlo en la ecuación (12):

J=E\left\lbrace \left[[F(z^{-1})v(k+d+1)+\left(\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)\right)-w(k+d+1)\right]^2\right\rbrace

El procedimiento se hace exactamente igual al hecho en el caso anterior, con eso se llega a que para minimizar J debemos igualar la diferencia entre la predicción de salida y el setpoint a cero:

\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)-w(k+d+1)=0

Si resolvemos la ecuación (10) para u(k) tenemos que:

\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)-w(k+d+1)=0

u(k)=\dfrac{1}{zB(z^{-1})F(z^{-1})}\left[C(z^{-1})w(k+d+1)-G(z^{-1})y(k)\right]\ [14]

La ecuación (14) corresponde al control por mínima varianza, MVC, con seguimiento de referencias. En esta ecuación también observamos que ocurre los cancelamientos de ceros del proceso en lazo abierto, lo que nos indica que la estructura tiene problemas a la hora de tratar con sistemas de fase no mínima.

La siguiente figura representa el diagrama de bloques correspondiente al sistema de control de minina varianza con ley de control dada por la ecuación (14)

control por mínima varianza con seguimiento de referencia

Controlador por Asignación de Polos – RST Posicional

Control Predictivo Basado en Modelo – DMC

CONTROLABILIDAD Y OBSERVABILIDAD

Control por Mínima Varianza Generalizado (GMV) Metodo I

En este caso además de incluir el error en la función costo, incluimos una ponderación sobre la acción de control. Así para minimizar la varianza utilizamos el siguiente indice:

J=E\{[y(k+d+1)-w(k+d+1)]^{2}+ru^2(k)\}\ [15]

Voy a mostrar como hacer el procedimiento para encontrar el mínimo de la ecuación (15), el procedimiento es el mismo que se puede aplicar en las dos estrategias anteriores. Se sabe que para encontrar el mínimo del indice con el objetivo de determinar la ley de control, debemos derivar J con relación a u(k) e igualarlo a cero:

\dfrac{\partial J}{\partial u(k)}=0

Vemos que en la expresión de J tenemos la ecuación y(k+d+1), la cual ya sabemos por la ecuación (9) que es:

y(k+d+1)=F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)

y que podríamos expresarla como:

C(z^{-1})y(k+d+1)=F(z^{-1})C(z^{-1})v(k+d+1)+G(z^{-1})y(k)+zB(z^{-1})F(z^{-1})u(k)

Sabemos que B, C, F y G son polinomios, habría que sustituir cada uno de los polinomios, para poder despejar únicamente el y(k+d+1). Para demostrar como obtener el mínimo únicamente voy a sustituir el C que acompaña la y(k+d+1) y el zBF que acompaña el u(k).

(1+c_1z^{-1}+...+c_{n_c}z^{-n_c})y(k+d+1)=\\F(z^{-1})C(z^{-1})v(k+d+1)+G(z^{-1})y(k)+z(b_1+b_2z^{-1}...+b_{n_b}z^{-n_b})(1+...+f_{d}z^{-d})u(k)

de la ecuación anterior podemos despejar $y(k+d+1)$ para poderla sustituir en la función objetivo de la ecuación (15), recordemos que como queremos mínimizar con relación a $u(k)$ , si resolvemos el producto de todos los polinomios de la ecuación anterior el único coeficiente que va a acompañar a $u(k)$ sería el primer termino del polinomio B o sea $b_1$ . Así en términos muy generales la ecuación anterior queda más o menos así:

y(k+d+1)=-[C(z^{-1})-1]y(k+d+1)+F(z^{-1})C(z^{-1})v(k+d+1)+G(z^{-1})y(k)\\+b_2u(k)+[b_2u(k)-zB(z^{-1})F(z^{-1})u(k)]

Ya ustedes pueden hacer todo el procedimiento en casa, lo importante notar aquí es que la ecuación (9) y (16) son equivalentes, es decir son las mismas ecuaciones expresadas de forma diferente, por lo pronto, si procedemos a minimizar el indice J sustituyendo la ecuación (9) tenemos que:

J=E\left\lbrace \left[[F(z^{-1})v(k+d+1)+\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)-w(k+d+1)\right]^2+ru^2(k)\right\rbrace

Podemos cancelar de una vez el ruido blanco porque sabemos que su esperanza matemática es cero entonces tendriamos que:

J=E\{[y(k+d+1)-w(k+d+1)]^{2}+ru^2(k)\}=

E\left\lbrace \left[ \dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)-w(k+d+1)\right]^2+ru^2(k)\right\rbrace

y por la ecuación (16) sabemos que (este paso únicamente es para ver como hacer la derivada más fácilmente)

J=E\{[y(k+d+1)-w(k+d+1)]^{2}+ru^2(k)\}=

E\{[ -[C(z^{-1})-1]y(k+d+1)+F(z^{-1})C(z^{-1})v(k+d+1)+G(z^{-1})y(k)\\+b_2u(k)+[b_2u(k)-zB(z^{-1})F(z^{-1})u(k)-w(k+d+1)]^2+ru^2(k)\}

hallamos el mínimo de la ecuación anterior derivando J con relación a $u(k)$ y recordando que la ecuación (16) y la (9) son iguales, solo que la (9) es más compacta, entonces expresamos el resultado de la forma de dicha ecuación:

\dfrac{\partial J}{\partial u} = 2b_2 \left[\dfrac{G(z^{-1})}{C(z^{-1})}y(k)+\dfrac{zB(z^{-1})F(z^{-1})}{C(z^{-1})} u(k)-w(k+d+1)\right]+2ru(k)=0

Si resolvemos la ecuación (17) para $u(k)$ tenemos que:

u(k)=\dfrac{1}{zB(z^{-1})F(z^{-1})+\dfrac{r}{b_1} C(z^{-1})}\left[C(z^{-1})w(k+d+1)-G(z^{-1})y(k)\right]

El siguiente es el diagrama de bloques del control por minima varianza generalizado o ponderado:

Generalised Minimum Variance control

Control por Mínima Varianza Generalizado (GMV) Metodo II

En este caso además de incluir el error en la función costo, incluimos una ponderación sobre la acción de control.

Esta estructura permite solucionar los problemas del controlador de Mínima Varianza estandar, como por ejemplo el controlar procesos de fase no mínima y seguir referencia sin error de offset. Para eso, se introduce en el sistema una Pseudo-Salida $\phi(k+d)$ definido como:

]\phi(k+d)=Py(k+d)+Qu(k)-Rw(k) [15]

Control por mínima varianza generalizado

El cual puede interpretarse como la salida generalizada del sistema, ya que se involucra una accion feedforward, se filtra la acción de salida y la accion de control se aplica sobre el error del sistema.

El proyecto del controlador consiste en minimizar la salida del sistema generalizado

J=E\{[\phi(k+d)]^{2}\}\ [16]

Sabemos que el sistema Original viene dado por el modelo ARMAX

A(z^{-1})y(k)=z^{-d}B(z^{-1})u(k)+C(z^{-1})v(k)

Atrasando la señal [/latex]d[/latex] pasos para al frente y despejando la salida tenemos que:

y(k+d)=\dfrac{B}{A}u(k)+\dfrac{C}{A}v(k+d)

si sustituimos esta expresión en la ecuación (15):

\phi(k+d)=\dfrac{PA+AQ}{A} u(k+d)+\dfrac{PC}{A} v(k+d)-Rw(k)

De la misma forma hecha en MV, el controlador GMV utiliza la siguiente identidad polinomial:

PC=FA+z^{-d}G

F(z^{-1})=1+f_1z^{-1}+f_2z^{-2}+...+f_{d}z^{-(d-1)}

G(z^{-1})=g_0+g_1z^{-1}+g_2z^{-2}+...+g_{n_a-1}z^{-n_g}

n_g=max(n_a-1,n_p+n_c-d)

Multiplicando el sistema original

Ay(k+d)=Bu(k)+Cv(k+d)

por F y reemplazando FA por la ecuación (16) tenemos:

PCy(k+d)=BFu(k)+Gy(k)+CFv(k+d)

Si se adiciona la expresión $QCu(k)-CRw(k)$ a ambos lados de la ecuación de arriba:

PCy(k+d)+CQu(k)-CRw(k)=BFu(k)+Gy(k)+CFv(k+d)+CQu(k)-CRw(k)

C[Py(k+d)+Qu(k)-Rw(k)]=(BF+CQ)u(k)+Gy(k)-Rw(k)+Fv(k+d)

\phi(k+d)=\dfrac{1}{C}\left[ (BF+CQ)u(k)+Gy(k)-Rw(k)\right]+Fv(k+d)\ [17]

La expresión dentro de corchetes de (17) es compuesto por informaciones conocidas en el instante k, y el termino $Fv(k+d)$ es compuesto por informaciones futuras no disponibles en el instante k y no correlacionados con ningún otro termino de la expresión.

Para minimizar la ecuación (16) se hace igualando la primera parcela de (17) a cero. Así la ley de control del GMV viene dado por:

(BF+CQ)u(k)+Gy(k)-Rw(k)=0

u(k)=\dfrac{CRw(k)-Gy(k)}{Bf+CQ}

EJEMPLOS

Diseñar un control por mínima varianza, mínima varianza con referencia, con integrador y con seguimiento de referencia para el siguiente proceso inestable:

$P(z)=\dfrac{2 z^{-1}}{1 - 2 z^{-1}}z^{-1}$

Con una perturbación medible con la siguiente función de transferencia:

$Q(z)=\dfrac{1 + 0.5 z^{-1}}{1 - 2 z^{-1}}$

2. Hacer un controlador por mínima varianza generalizado para el siguiente proceso inestable tipo integrador:

$P(z)=\dfrac{z^{-1} + 1.5 z^{-2}}{1 - 1.7 z^{-1} + 0.7 z^{-2}}z^{-1}$

Con la siguiente perturbación medible:

$Q(z)=\dfrac{1 - 0.9z^{-1}}{1 - 1.7 z^{-1} + 0.7 z^{-2}}$

Los dos Ejemplos están resuletos en MATLAB, para poder ver los códigos simplemente basta con compartir el conténido de este post utilizando los botones de compartir que están a continuación o suscribiendose al canal de youtube. De esa forma apoyas este sitio web y permites que más personas se beneficien con esta información.

Para poder ejecutar los siguientes códigos, es necesario descargar, descomprimir y guardar la siguiente función llamada «PolFG.m» en la misma carpeta donde van a guardar los códigos o en la carpeta de MATLAB que se encuentra en «Mis Documentos«. Esta función permite hacer el calculo de los polinomios F y G.

Descargar Función PolFG

Control Por Mínima Varianza

% Control Por Mínima Varianza (Extraido del Libro Self Tuning System)
% Sergio Andres Castaño Giraldo
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - Medellín - 2017
% https://controlautomaticoeducacion.com
% ------------------------------------------------------
clc
clear all
close all

Ts=1; %Tiempo de Muestreo 

% Función de Trasnferencia del Proceso
P=filt([0 2],[1 -2],Ts);
d=1; %Retardo
P.iodelay=d
% Función de Transferencia de la Perturbación
Q=filt([1 0.5],[1 -2],Ts)

%Obtengo los polinomios del modelo ARMAX
[B,A]=tfdata(P,'v');
C=tfdata(Q,'v');

A=[A zeros(1,d)];

% d=10;
% B=[0 1 0.5];
% A=[1 -1.7 0.7 zeros(1,d)];
% C=[1 -0.9 0 0];


m=length(A)-1-d; %Orden del Sistema

%% calculo de poloinomio F,G y polinomio FB
[F,G]=PolFG(A,C,d,m);
%Producto entre B y F
zBF=conv(B(2:end),F);

%Simula el proceso en Simulink
% sim('MVSim');

%% Simula el proceso
nit=400;  %Numero de Interacciones
Var=0.5; %Varianza del Ruido
e = sqrt(Var)*randn(nit,1); %Ruido Blanco
y(1:nit)=0; %Salida del Proceso
u(1:nit)=0; %Acción de Control

y1(1:nit)=0; %Salida del Proceso MA
u1(1:nit)=0; %Acción de Control MA

for k=4+d:nit
    %Minima Varianza
   y(k)=-A(2)*y(k-1)+B(2)*u(k-1-d)+C(1)*e(k)+C(2)*e(k-1);
   u(k)=(-G*y(k)-zBF(2)*u(k-1))/zBF(1);
   %Malla Abierta
   y1(k)=-A(2)*y(k-1)+B(2)*u1(k-1-d)+C(1)*e(k)+C(2)*e(k-1);
end
t=0:Ts:(nit-1)/Ts;
figure
plot(t,y,t,y1),grid
legend('Minima Varianza','Lazo Abieto')

Control Por Mínima Varianza Sobre el Error Estandar Vs Incremental

% Control Por Mínima Varianza Estandar Vs Incremental
%(Extraido del Libro Self Tuning System)
% Sergio Andres Castaño Giraldo
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - Medellín - 2017
% https://controlautomaticoeducacion.com
% ------------------------------------------------------
clc
clear all
close all

Ts=1; %Tiempo de Muestreo 

% % Función de Trasnferencia del Proceso
% P=filt([0 1],[1 -0.8],Ts);
% d=2; %Retardo
% P.iodelay=d
% % Función de Transferencia de la Perturbación
% Q=filt([1 0.98],[1 -0.8],Ts)

% % Función de Trasnferencia del Proceso
% P=filt([0 2],[1 -0.32],Ts);
% d=1; %Retardo
% P.iodelay=d
% % Función de Transferencia de la Perturbación
% Q=filt([1 0.5],[1 -0.32],Ts)

% Función de Trasnferencia del Proceso
P=filt([0 2],[1 -2],Ts);
d=1; %Retardo
P.iodelay=d
% Función de Transferencia de la Perturbación
Q=filt([1 0.5],[1 -2],Ts)

%Obtengo los polinomios del modelo ARMAX
[B,Ap]=tfdata(P,'v');
C=tfdata(Q,'v');

A=[Ap zeros(1,d)]; %Caso Sin Integrador
%Adición del Integrador en el polinomio A
Ai=[conv(Ap,[1 -1]) zeros(1,d)];

% d=10;
% B=[0 1 0.5];
% A=[1 -1.7 0.7 zeros(1,d)];
% C=[1 -0.9 0 0];


m=length(A)-1-d; %Orden del Sistema
mi=length(Ai)-1-d; %Orden del Sistema Integrador


%% calculo de poloinomio F,G y polinomio FB
[F,G]=PolFG(A,C,d,m); 
[Fi,Gi]=PolFG(Ai,C,d,mi); % Caso Integrador


%Producto entre B y F
zBF=conv(B(2:end),F);
zBFi=conv(B(2:end),Fi);


%% Simula el proceso
nit=400;  %Numero de Interacciones
Var=0.5; %Varianza del Ruido
e = sqrt(Var)*randn(nit,1); %Ruido Blanco
er(1:nit)=0; %Error entre Referencia y Salida
y(1:nit)=0; %Salida del Proceso
u(1:nit)=0; %Acción de Control

eri(1:nit)=0; %Error entre Referencia y Salida Integrador
yi(1:nit)=0; %Salida del Proceso Caso Integrador
ui(1:nit)=0; %Acción de Control Caso Integrador
du(1:nit)=0; %Incremento de Control
%Referencia
w(1:100)=4;w(100:200)=-4;w(200:300)=4;w(300:400)=-4;

for k=4:nit
   %Salida del Proceso Mínima Varianza Estandar
   y(k)=-A(2)*y(k-1)+B(2)*u(k-1-d)+C(1)*e(k)+C(2)*e(k-1);
   %Salida del Proceso Mínima Varianza Incremental
   yi(k)=-A(2)*yi(k-1)+B(2)*ui(k-1-d)+C(1)*e(k)+C(2)*e(k-1);
   
   %Error de Referencia
   er(k)=w(k)-y(k);
   eri(k)=w(k)-yi(k);
   
   %Control de Mínima Varianza Estandar
   u(k)=(G*er(k)-zBF(2)*u(k-1))/zBF(1);
   %Control de Mínima Varianza Incremental
   du(k)=(Gi(1)*eri(k)+Gi(2)*eri(k-1)-zBFi(2)*du(k-1))/zBFi(1);
   ui(k)=du(k)+ui(k-1);
end
t=0:Ts:(nit-1)/Ts;
figure
%Mínima Varianza Estandar
subplot(2,1,1)
plot(t,y,t,w),grid
title('Control de Mínima Varianza Estandar');
xlabel('Tiempo(Muestras)');
%Mínima Varianza Incremental
subplot(2,1,2)
plot(t,yi,t,w),grid
title('Control de Mínima Varianza Incremental');
xlabel('Tiempo(Muestras)');
axis([0 nit -40 40])

Control Por Mínima Varianza con Seguimiento de Referencia Vs Incremental

% Control Por Mínima Varianza (Extraido del Libro Self Tuning System)
% Con Seguimiento de Referencia Vs Incremental
% Sergio Andres Castaño Giraldo
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - Medellín - 2017
% https://controlautomaticoeducacion.com
% ------------------------------------------------------
clc
clear all
close all

Ts=1; %Tiempo de Muestreo 

% Función de Trasnferencia del Proceso
P=filt([0 2],[1 -2],Ts);
d=1; %Retardo
P.iodelay=d
% Función de Transferencia de la Perturbación
Q=filt([1 0.5],[1 -2],Ts)

%Obtengo los polinomios del modelo ARMAX
[B,Ap]=tfdata(P,'v');
C=tfdata(Q,'v');

A=[Ap zeros(1,d)]; %Caso Sin Integrador
%Adición del Integrador en el polinomio A
Ai=[conv(Ap,[1 -1]) zeros(1,d)];

% d=10;
% B=[0 1 0.5];
% A=[1 -1.7 0.7 zeros(1,d)];
% C=[1 -0.9 0 0];


m=length(A)-1-d; %Orden del Sistema
mi=length(Ai)-1-d; %Orden del Sistema Integrador


%% calculo de poloinomio F,G y polinomio FB
[F,G]=PolFG(A,C,d,m); 
[Fi,Gi]=PolFG(Ai,C,d,mi); % Caso Integrador


%Producto entre B y F
zBF=conv(B(2:end),F);
zBFi=conv(B(2:end),Fi);


%% Simula el proceso
nit=400;  %Numero de Interacciones
Var=0.5; %Varianza del Ruido
e = sqrt(Var)*randn(nit,1); %Ruido Blanco
er(1:nit)=0; %Error entre Referencia y Salida
y(1:nit)=0; %Salida del Proceso
u(1:nit)=0; %Acción de Control

eri(1:nit)=0; %Error entre Referencia y Salida Integrador
yi(1:nit)=0; %Salida del Proceso Caso Integrador
ui(1:nit)=0; %Acción de Control Caso Integrador
du(1:nit)=0; %Incremento de Control
%Referencia
r(1:100)=4;r(100:200)=-4;r(200:300)=4;r(300:400)=-4;

for k=4:nit
   %Salida del Proceso Mínima Varianza con Seguimiento de Referencia
   y(k)=-A(2)*y(k-1)+B(2)*u(k-1-d)+C(1)*e(k)+C(2)*e(k-1);
   %Salida del Proceso Mínima Varianza Incremental
   yi(k)=-A(2)*yi(k-1)+B(2)*ui(k-1-d)+C(1)*e(k)+C(2)*e(k-1);
   
   %Error de Referencia
   er(k)=C(1)*r(k)+C(2)*r(k-1)-G*y(k);
   eri(k)=r(k)-yi(k);
   
   %Control de Mínima Varianza con Seguimiento de Referencia
   u(k)=(er(k)-zBF(2)*u(k-1))/zBF(1);
   %Control de Mínima Varianza Incremental
   du(k)=(Gi(1)*eri(k)+Gi(2)*eri(k-1)-zBFi(2)*du(k-1))/zBFi(1);
   ui(k)=du(k)+ui(k-1);
end
t=0:Ts:(nit-1)/Ts;
figure
%Mínima Varianza con Seguimiento de Referencia
subplot(2,1,1)
plot(t,y,t,r),grid
title('Control de Mínima Varianza con Seguimiento de Referencia');
xlabel('Tiempo(Muestras)');
%Mínima Varianza Incremental
subplot(2,1,2)
plot(t,yi,t,r),grid
title('Control de Mínima Varianza Incremental');
xlabel('Tiempo(Muestras)');
axis([0 nit -40 40])

Control Por Mínima Varianza Generalizado

% Control Por Mínima Varianza (Extraido del Libro Self Tuning System)
% Generalizado
% Sergio Andres Castaño Giraldo
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - Medellín - 2017
% https://controlautomaticoeducacion.com
% ------------------------------------------------------
clc
clear all
close all

Ts=0.1; %Tiempo de Muestreo 
% % Función de Trasnferencia del Proceso
% P=filt([0 2],[1 -2],Ts);
% d=1; %Retardo
% P.iodelay=d
% % Función de Transferencia de la Perturbación
% Q=filt([1 0.5],[1 -2],Ts)

% Función de Trasnferencia del Proceso
P=filt([0 1 1.5],[1 -1.7 0.7],Ts);
d=1; %Retardo
P.iodelay=d
% Función de Transferencia de la Perturbación
Q=filt([1 -0.9 0 0],[1 -1.7 0.7],Ts)

%Obtengo los polinomios del modelo ARMAX
[B,Ap]=tfdata(P,'v');
C=tfdata(Q,'v');

A=[Ap zeros(1,d)]; %Caso Sin Integrador

m=length(A)-1-d; %Orden del Sistema

%% calculo de poloinomio F,G y polinomio FB
[F,G]=PolFG(A,C,d,m); 


%Producto entre B y F
zBF=conv(B(2:end),F);


%% Simula el proceso
nit=400;  %Numero de Interacciones
Var=0.5; %Varianza del Ruido
e = sqrt(Var)*randn(nit,1); %Ruido Blanco
er(1:nit)=0; %Error entre Referencia y Salida
y(1:nit)=0; %Salida del Proceso
u(1:nit)=0; %Acción de Control
du(1:nit)=0; %Incremento de Control

%Referencia
w(1:100)=4;w(100:200)=-4;w(200:300)=4;w(300:400)=-4;
r=0.5;
%dCon=zBF+r*C(1);
 dCon=zBF+r/B(2)*C;
for k=4:nit
   %Salida del Proceso Mínima Varianza Generalizado
    y(k)=-A(2)*y(k-1)-A(3)*y(k-2)+B(2)*u(k-1-d)+B(3)*u(k-2-d)+C(1)*e(k)+C(2)*e(k-1);
    
   %Error de Referencia
   er(k)=C(1)*w(k)+C(2)*w(k-1)-G(1)*y(k)-G(2)*y(k-1);
   
   %Control de Mínima Varianza Generalizado
   u(k)=(er(k)-dCon(2)*u(k-1)-dCon(3)*u(k-2))/dCon(1);

end
t=0:Ts:(nit-1)*Ts;
figure
%Mínima Varianza Generalizado
plot(t,y,t,w),grid
title('Control de Mínima Varianza Generalizado');
xlabel('Tiempo(Muestras)');

Eso es todo por la entrada del dia de hoy, espero les haya gustado y hayan aprendido algo nuevo. Si te ha servido el contenido de esta entrada, de los videos y los códigos de implementación y deseas apoyar mi trabajo invitandome a un café super barato, puedes hacerlo en el siguiente link:

👉 Invitar a Sergio a un Café ☕️

Que esten muy bien, nos vemos en la siguiente entrada.

Volver al Curso de Control de Procesos Avanzados

Sergio A. Castaño Giraldo

Mi nombre es Sergio Andres Castaño Giraldo, y en este sitio web voy a compartir una de las cosas que mas me gusta en la vida y es sobre la Ingeniería de Control y Automatización. El sitio web estará en constante crecimiento, voy a ir publicando material sobre el asunto desde temas básicos hasta temas un poco más complejos. Suscríbete al sitio web, dale me gusta a la página en Facebook y únete al canal de youtube. Espero de corazón que la información que comparto en este sitio, te pueda ser de utilidad. Y nuevamente te doy las gracias y la bienvenida a control automático educación.