ER4.1. Control discreto de asignación de polos por enfoque polinomial

Introducción (planteamiento)

Dado el siguiente modelo lineal de un proceso, diseñar y simular un controlador con los siguientes requerimientos de diseño para una referencia tipo escalón: $M_p\leqslant 10\%, t_s\simeq 0.8 \mathrm{seg}, e_{ss}=0$. Realizar un análisis de sensibilidad e incertidumbre.

$G_p(s)=\frac{1}{s(s+1)}\\ T_s=0.1 \mathrm{seg}$

Método (plan de solución)

Pasos: (i) discretizar el modelo; (ii) dados los requerimientos de diseño especificar la función de transferencia deseada en lazo cerrado de tiempo continuo utilizando las características de un sistema de segundo orden; (iii) discretizar el modelo deseado; (iv) calcular $\mu$ y $\nu$; (v) calcular el regulador discreto por asignación de polos; (vi) calcular el controlador de adelanto; (vii) simular el control con el proceso de tiempo continuo y graficar las variables de salida y control; (viii) realizar el análisis de sensibilidad e incertidumbre on la herramienta GSUA de MATLAB; (ix) analizar los resultados. La operaciones que se realizan a continuación se detallan en el código de MATLAB de este enlace.

Resultados (solución)

Discretización del modelo con la función c2d de MATLAB, donde el período de muestreo $T_s=0.1$ se calculó a partir del tiempo de crecimiento:

 

Se observa que el modelo tiene un integrador dado por el término $s$ en la función de transferencia continua y $(z-1)$ en la función de transferencia discreta, por lo que el error en estado estacionario para una referencia escalón será igual a cero (si no tuviera un integrador se podría adicionar uno "artificialmente" al modelo discreto de la planta y luego pasarlo al controlador). De acuerdo con los requerimientos de diseño y lo explicado en el libro, se tiene:

$M_p=\frac{y_{\max}-y_{ss}}{y_{ss}}100\%=e^{-\pi \zeta /\sqrt{1-\zeta ^2}}100\%$

$\zeta =\frac{1}{\sqrt{1+\left[ \pi /\ln \left( M_p/100 \right) \right] ^2}}=0.5912$

$t_s\simeq \frac{3.2}{\zeta \omega _0}\,\,  (0<\zeta <0.69)$

$\omega _0\simeq \frac{3.2}{\zeta t_s}=6.7659$

La función de transferencia deseada en MATLAB de tiempo continuo y tiempo discreto con las características anteriores es.

    

El orden del regulador es, con $n=2$:

$\mu =\nu =n-1=1$

Como el modelo deseado en lazo cerrado es de orden 2 y no de orden $2n-1=3$, se debe cancelar un polo estable de la planta con un cero del regulador. Siguiendo el procedimiento expuesto en el libro se llega al siguiente regulador, donde es necesario resolver la ecuación para hallar los parámetros $(q_0,q_1,p_1)$ del regulador con un término común $(Q^C)$ con el proceso $(A^c)$ para disminuir el orden del modelo en lazo cerrado:

$G_c(z)=\frac{q_0z+q_1}{z+p_1}$

Para el proceso:

$G_p(z)=\frac{0.004837(z+0.9672)}{(z-1)\color{red}{(z-0.9048)}}$

Para el modelo deseado:

$G_m(z)=\frac{0.172z+0.1314}{z^2-1.146z+0.4493}=\frac{B_m}{A_m}$

Se obtiene la siguiente ecuación:

$\color{red} {A^c}A^{nc}P^cP^{nc}+B\color{red}{Q^c}Q^{nc}=\color{red}{A^c}A_m$

$P^c=1$

De manera más simple,

$Q^{nc}\color{red}{Q^c}=q_0\color{red}{(z-0.9048)}$

$\color{red}{(z-0.9048)}(z-1)(z+p_1)+0.004837(z+0.9672)\color{red}{(z-0.9048)}q_0\\=\color{red} {(z-0.9048)}(z^2-1.146z+0.4493)$

Eliminando el término común de fase mínima se obtiene un sistema de dos ecuaciones con dos incógnitas (se pudo dejar la ecuación con tres incógnitas):

$(z-1)(z+p_1)+0.004837(z+0.9672)q_0=z^2-1.146z+0.4493$

Ecuaciones por resolver, comparando los términos a la izquierda y derecha en la expresión anterior y después de destruir paréntesis:

$z^1:  p_1-1+0.004837q_0=-1.146\\ z^0: -p_1+0.004678q_0=0.4493$

La solución es: 

$p_1=-0.3002\\ q_0=31.88$

El controlador es:

$G_c(z)=\frac{31.88(z-0.9048)}{z-0.3002}$

Si se utiliza la matriz de Sylvester se debe tener en cuenta lo siguiente, según lo explicado en el libro (si se incluye el término común en $A_m$ se tiene el mismo efecto anterior de cancelar un término común):

$A(z)=(z-1)(z-0.9048)=z^2-1.905z+0.9048$

$n=2, a_0=1, a_1=-1.905, a_2=0.9048$

$B(z)=0.004837z+0.004678, b_0=0, b_1=0.004837, b_2=0.004678$

$A_m(z)=(z^2-1.146z+0.4493)\color{red}{(z-0.9048)}=z^3-2.051z^2+1.486z-0.4065$

$a_{m,0}=1, a_{m,1}=-2.051, a_{m,2}=1.486, a_{m,3}=-0.4065$

$\mathbf{A}_m=\left[ \begin{array}{c} 1\\ -2.051\\ 1.486\\ -0.4065\\\end{array} \right]$

$\mathbf{E}=\left[ \begin{array}{cc:cc} b_0&  0&  a_0&  0\\ b_1&  b_0&  a_1&  a_0\\ \hdashline b_2&  b_1&  a_2&  a_1\\ 0&  b_2&  0&  a_2\\\end{array} \right] =\left[ \begin{array}{cc:cc} 0&  0&  1&  0\\ 0.004837&  0&  -1.905&  1\\ \hdashline 0.004678&  0.004837&  0.9048&  -1.905\\ 0&  0.004678&  0&  0.9048\\\end{array} \right]$

El controlador con $2n$ parámetros, el cual coincide con el resultado anterior dentro de los errores numéricos, es:

$\mathbf{M}=\left[\begin{array}{c} q_0\\ q_1\\ p_0\\ p_1\\\end{array} \right] = \mathbf{E}^{-1} \mathbf{A}_m=\left[ \begin{array}{c} 31.87\\ -28.84\\ 1\\ -0.3002\\\end{array} \right]$

$G_c(z)=\frac{31.87(z-0.9048)}{z-0.3002}$

El número de condición de la matriz de Sylvester es igual a 1052, un número aceptable. El control de adelanto es:

$G_f(z) = \frac{B_m(z)}{B^{nc}(z)Q^{nc}(z)} = \frac{0.172z+0.1314}{0.004837(z+0.9672) 31.87} = \frac{1.116(z+0.764)}{z+0.9672}$

Diagrama de simulación híbrido (proceso continuo y controlador discreto) en Simulink:

Los resultados, donde se observa un buen ajuste al modelo deseado y una acción de control que supondremos aceptables, son:

El problema también se puede resolver incluyendo el cero de la planta en lugar del polo en el polinomio deseado, en cuyo caso el diagrama de simulación es el siguiente, con idénticos resultados al simular (se deja al lector el cálculo):

 

El análisis de sensibilidad e incertidumbre con un 30% de incertidumbre en dos parámetros del modelo da los siguientes resultados.

Análisis de incertidumbre:

Análisis de sensibilidad:

Condiciones de método de análisis de sensibilidad e incertidumbre:

 

El controlador diseñado se puede comparar con un controlador que ubique todos los polos en $z=0$, lo cual corresponde a un controlador de oscilaciones muertas (deadbeat control). El diseño con MATLAB y los resultados son los siguientes, donde se observa un esfuerzo de control muy grande:

 

Diagrama de simulación con dicho controlador:

Resultados de la simulación, donde se observa un aumento considerable de la acción de control:

Discusión y verificación

El control entrega buenos resultados y los cálculos coinciden por los dos métodos utilizados para asegurar una solución única: utilizando un cero del controlador para cancelar un polo de la planta y por el método de la matriz de Sylvester. El control se ajusta bien al modelo de referencia, con lo cual se cumple con los requerimientos de diseño. Una ventaja del diseño es que el proceso tiene un integrador que elimina el error en estado estacionario a una entrada escalón; si no lo tuviera, se podía adicionar uno al modelo discreto del proceso y luego pasarlo al controlador, una manipulación matemática bastante interesante.

La acción de control toma unos valores máximos en valor absoluto que pueden ser o no aceptables y depende de si el sistema real admite tales amplitudes; en el caso de una simulación se pueden probar otros diseños que disminuyan dichos valores y, en general, al relajar las especificaciones (menor velocidad o mayor sobreimpulso máximo) se puede bajar el esfuerzo de control.

El análisis de incertidumbre muestra que el sistema de control funciona bien a pesar de una incertidumbre del 30 % en dos parámetros del modelo, en cuyo caso hay un aumento del sobreimpulso máximo y del tiempo de establecimiento, pero el sistema se estabiliza. El análisis de sensibilidad indica que el parámetro que más incide en el ajuste a la curva de referencia es la ganancia del modelo, por lo cual debe prestarse especial atención a su medición más precisa.

Para comparar los resultados, se diseña un controlador que asigna todos los polos en el origen, lo cual corresponde al controlador más rápido posible (controlador de oscilaciones muertas), lo cual efectivamente se observa en la respuesta temporal del sistema en lazo cerrado (el sistema se estabiliza en 0.3 seg, menor que los 0.8 seg del primer diseño), pero con una acción de control con valores mucho más grandes del orden de 300, contra valores de no más de 40. En general, al mover más los polos de la planta continua a la izquierda en el semiplano izquierdo (hacia el centro del círculo unitario en el caso discreto) se incurre en un mayor esfuerzo de control. Aunque todavía se pueden ubicar los ceros para reducir el sobreimpulso máximo, se deja esa tarea al lector.