home

=Somos Rafael Molina López, Pedro Pérez Alcántara y Nikita Popov Popov.= =Pertenecemos al grupo 8 de la asignatura Control Automático, de 3º ETSII.= =En esta wiki presentamos un problema de control en el que desarrollamos múltiples ideas de la asignatura.=

Para este trabajo, se nos ha proporcionado una función de MATLAB que nos asigna una función de transferencia en función del DNI del compañero más joven.



El compañero más joven del grupo es Nikita Popov Popov. Llamamos dniNICK a su DNI separando cada cifra. Entonces esto nos proporciona la función de transferencia.



//**1) Hallar la función de transferencia en z del sistema discretizado con T = 0.1, y obtener la ecuación en diferencias correspondiente, y especificar los polos y ceros del sistema discretizado.**//

Para hacer una buena conversión, vamos a trabajar con matemática exacta, en lugar de los decimales truncados. Así evitamos errores de diferencias cancelativas.



Procedemos a la discretización. Recordemos que tenemos que incluir el efecto del muestreador de HOLD 0.



Para una transformación más cómoda, pasamos la //s// del denominador del muestreador a la función de trasnferencia., así la transformada //Z// del numerador es directa.



Hacemos la transformada //Z// del primer término y descomponemos en factores simples lo que nos falta. Averiguamos las constantes por el método del residuo.



Entonces nos queda, antes de transformar, el siguiente resultado.



Todas las transformaciones son directas. Podemos obtenerlas mediante las tablas de trasnformada //Z//.



Entonces, antes de operar, nos queda el siguiente resultado.



Una vez hacemos la operación completa, al final podemos escribir el resultado decimal y prescindir de la exacta, ya que tendríamos la función de transferencia completa.



Hemos tenido suerte. Tenemos una función de transferencia con todos los polos y ceros de estabilidad. En //s// en el semiplano negativo y en //z// en el círculo unidad.


 * Polos: z = 0.7225; z = 0.7286
 * Ceros: z = 0.2935

A continuación, vamos a obtener la ecuación en diferencias correspondiente al sistema.

La función de transferencia de un sistema en tiempo discreto es la relación entre la transformada //Z// de la salida y la transformada //Z// de la entrada con condiciones iniciales nulas. Para un sistema LTI característico: Aplicando esto a nuestro caso tenemos: //**2) Obtener la descripción en espacio de estado G, H, C y D del sistema.**//

Hemos obtenido en el apartado anterior que la función de transferencia es la siguiente



Vamos a hacer la descripción en espacio de estado utilizando la representación en **Forma Canónica de Control**.



Los términos bn representan los coeficientes que acompañan al término de grado m-n en el numerador. Los términos an funcionan igual, pero para el denominador. Dicho esto, vemos que b0 vale cero, debido a que el numerador tiene un grado menor que el denominador. Entonces, tenemos la descripción completa.




 * //3) Hallar la respuesta escalón del sistema y calcular el valor en régimen permanente.//**

Esta es una tarea sencilla. Si para algo sirve la función de transferencia es para esto. Si multiplicamos la señal de entrada por la función de trasferencia obtenemos:



Para hallar el valor en régimen permanente se calcula el siguiente límite (R(z) denota la señal de salida): Por tanto el valor en régimen permanente es 7.2635.

//**4)Añadir un retardo de t 0 =0,17 en serie con la planta original y determinar la función de transferencia discretizada, hallando la respuesta escalón. **//

Para añadir un retardo dividimos el tiempo de retardo en dos partes, por un lado un múltiplo del tiempo de muestreo y por otro una fracción del mismo, tambien calculamos el valor de m que nos sera necesario a la hora de calcular la transformada Z modificada:

Introducimos el retardo en la función de transferencia y aplicamos la transformada en Z, es muy importante no olvidar el controlador hold 0: Descomponemos el segundo termino de la ecuación en fracciones simples, al igual que en el primer ejercicio y lo obtenido lo multiplicamos por la parte fraccional del retardo obteniendo: El termino e no se tiene en cuenta ya que se incluirá en la transformada modificada.Todas las transformaciones son directas y las obtenemos a partir de la tabla de la transformada en Z modificada: Que es la función de transferencia con el retardo incluido. La salida del sistema será: Para obtener la salida escalón la entrada debe ser: Y la salida que obtenemos es: //**5) Determinar el rango de estabilidad del sistema propuesto (sin retardo) cuando se controla el mismo con una ganancia variable K en configuración de realimentación negativa unitaria, utilizando para ello el criterio de Routh-Hurwitz por aplicación de la trasnformación bilineal.**//

Sabemos que la estabilidad de un sistema está asegurada cuando se cumple la siguiente condición.



El problema que se nos presenta es que para la función discretizada en //z// no podemos aplicar el criterio de Routh como lo hacíamos en continua. Para ello, haremos un traslado al plano //w//, es decir, convertiremos el círculo unidad en el eje imaginario //jv//. Para ello aplicamos la transformación bilineal



En primer lugar, vamos a obtener la función de transferencia del sistema en bucle cerrado, con realimentación negativa unitaria.



Y operando se obtiene

Del lazo cerrado podemos obtener el polinomio característico como el denominador de esta función de transferencia.

Si ahora aplicamos la transformación bilineal para //w// tenemos

Si multiplicamos la ecuación por el denominador del término cuadrático nos queda

Operando para expandir el polinomio característico en //w// obtenemos //P(w)//: Con este polinomio característico podemos empezar a hacer el análisis para el criterio de Routh-Hurwitz. Lo primero es asegurar que los coeficientes son mayor que 0.

El polinomio característico es solo de orden 2. En el momento en el que todos los términos son positivos, si hay raíces complejas, la parte real negativa está asegurada. Entonces, hemos acabado el criterio de Routh-Hurwitz. El sistema será estable si y sólo si //K// está en el intervalo //**6) Calcular el error en régimen permanente en función de K ante entrada escalón.**// Para calcular el error en régimen permanente ante entrada escalón lo primero que hacemos es calcular la K_p: Por tanto la K_p es: Ahora podemos calcular el error en régimen permanente: //**7)Calcular para k=1,5 la situación de los polos del sistema discretizado en bucle cerrado y ubicarlos sobre el plano z, comprobando la estabilidad del mismo.**// Primero comprobamos que k=1,5 está en el rango de estabilidad calculado en el ejercicio 5 y vemos que -0,1377<k<2,97 con lo cual el sistema es estable. Una vez comprobado que el sistema es estable sustituimos k en la ecuación de bucle cerrado calculada anteriormente: Y nos da como resultado:
 * Polos z1= 0,14468+0,4056j, z2=0,14468+0,4056j.
 * Zero z0=0,2934.

La ubicación en el plano Z sería: //**8) Diseñar un compensador aplicando el método de Truxal-Ragazzini ante entrada escalón, y graficar la respuesta.**// Para diseñar un controlador mediante el método de Truxal-Ragazzini, lo primordial es pensar en //G(z)// en potencias de //z// negativas.

Este resultado es fácil de obtener. Tan solo hay que multiplicar y dividir por la //z// de mayor grado, //z// al cuadrado. Truxal-Rugazzini dice que la función de transferencia en lazo cerrado será la siguiente.

Esta //F(z)// debe quedar configurada de la siguiente forma. Así tenemos una función de transferencia rapidísima (todos los polos en el cero) y una serie de ceros que compensan el sistema. La finalidad de todo esto es poder despejar de la fórmula anterior el controlador.

El controlador obtenido aplicando Truxal-Ragazzini debe cumplir 4 condiciones básicas:
 * //Condición de Realización física.//

Esta condición es simple. Lo único que dice es que, si la //G(z)// en desarrollo de potencias de //z//, la //F(z)// debe empezar igual, para no llevar adelanto. Es fácil ver cómo la función de transferencia //G(z)// en desarrollo de potencias empieza en //1/z//. En ese caso obtenemos
 * //Condición de estabilidad//

En realidad esta condición son dos condiciones al mismo tiempo. Ninguna de estas dos condiciones es necesaria, ya que el sistema es estable de por sí.
 * 1) 1 - //F(z)// debe tener por ceros los polos de //G(z)// que estén fuera del círculo unidad.
 * 2) //F(z)// debe tener por ceros los ceros de //G(z)// que estén fuera del círculo unidad.
 * //Error en régimen permanente nulo//

Esto solo se consigue añadiendo integradores a la //F(z)//. El número de integradores depende de la señal de entrada.



En nuestro caso, la entrada es escalón. Así que la condición queda configurada como:



//N(z)// es una función manipulable a nuestro antojo en principio. Y nuestra señal de entrada escalón es:




 * //Condición de rizado nulo//

Esta condición tiene un poco más de trabajo. El objetivo es que, una vez conseguido el valor final, permanezca constante, no oscile.

La condición depende de la existencia de integradores en la //G(z).// Si tenemos integradores, //U(z)// (nuestra señal de control) debe ser truncada. Si no los tiene, no.

Es decir, a partir de la última muestra correcta, en presencia de integradores, las siguientes muestras serán cero. Si no hay integradores, serán la misma cte.

No tenemos integradores en nuestra función de transferencia //G(z)//. Entonces //F(z)// NO DEBE SER DIVISIBLE POR //G(z)//.

//F(z)// no debe ser divisible entre los términos del denominador entre paréntesis, pero conocemos la estructura de //F(z).// Sabiendo esto, //F(z)// se trunca en //1/z.// Esto es así porque solo tenemos un término, el término de la condición de erp = 0.

Para averiguar el valor del término de //F(z)// efectuamos la división, que debe ser exacta.



Entonces,

Si aplicamos esto a la fórmula donde despejamos con anterioridad el controlador tenemos la función de transferencia del mismo. Sustituimos, operamos y nos queda:



En el proceso de creación del controlador hemos visto que el método está orientado a una entrada particular, con lo que hay una leve pérdida de generalidad.

Ahora debemos hacer la comprobación de la condición de rizado. Es decir, debemos ver si es cierto que //F(z)// no es divisible por //G(z)//.

Vamos a graficar ahora la salida. Hemos metido un escalón, así que multiplicamos la señal escalón por la nueva función de transferencia //F(z)// queda



El resultado es una división. Si operamos, el resultado (el cociente) es la salida. Cada término significa una muestra diferente. Desde la primera muestra, el error es 0.

Como resultado podemos decir que el controlador no sólo controla, sino que lo hace rapidísimo, desde la primera muestra. El dibujo sería el siguiente.




 * PARTE PRÁCTICA**

//**1) Programar el sistema obtenido en SIMULINK utilizando la planta original en continua, obteniendo la respuesta escalón y la respuesta rampa, graficandolas ambas para valores de T=0.1, T=0.5 y T=1, comparando las respuestas.**// Para ver como responde el sistema ante entrada __escalón__ montamos el siguiente esquema en el SIMULINK: Donde ZOH tiene el T=0.1, ZOH1 -> T=0.5 y ZOH2 -> T=1. Las respuestas obtenidas son las siguientes:

Para ver como responde el sistema ante entrada __rampa__ montamos el siguiente esquema en el SIMULINK:

Y obtenemos las siguientes salidas: Podemos observar que todas las gráficas son parecidas y que las respuestas de los sistemas discretizados por medio de Hold 0 (con distintos tiempos de muestreo) se parecen bastante a la respuesta continua.


 * //2) Expresar el sistema discretizado asignado para T=0.1 en formato tf,zpk y ss y especificar el diagrama de polos y ceros en z.//**

Introducimos la función en continua en formato zpk y la discretizamos: El resultado es la función discretizada en formato zpk. Al no especificar ningún método de discretización Matlab ha utilizado 'zoh'. Para pasar al formato tf solo será necesario hacer la conversión: Para pasar a espacio de estados (ss): El diagrama de polos y ceros en z será:



//**3) Aplicar ltiview al sistema discretizado asignado y calcular las especificaciones de respuesta transitoria (t subida, t pico, t establecimiento, Sobreoscilación) ante entrada escalón unitario con T=0.1 y T=1 y compararlas.**//

Primero vamos a utilizar el comando //ltiview// para la función de transferencia discretizada con T = 0.1.



Habiendo definido Gz con MATLAB anteriormente como la función de transferencia discretizada para T = 0.1 es tan fácil como escribir los comandos que se ven.

¿Cómo hemos podido ver todos los datos? Es fácil. Tras escribir //ltiview(Gz)// obtenemos el primer diagrama. Hacemos clic derecho, Characteristics, una opción.



De esta forma nos aparece un circulo en el punto del mapa de la respuesta donde se alcanza la especificación. Con solo poner el ratón encima aparece el valor.

Vemos que para este sistema las especificaciones que cumple son:
 * Sobreoscilación: 0%
 * Tiempo de subida: 1.1
 * Tiempo de establecimiento: 1.8
 * Tiempo de pico > 3

El tiempo de pico es > 3, pero eso no nos da el dato real. Para verlo tenemos que alejar el zoom en el eje horizontal. Hacemos clic derecho, Properties.



Se nos abre el menú opciones. En la pestaña //Limits// podemos cambiar el final del eje de abscisas, de 3 a 8. De esta forma sí podemos ver el tiempo de pico.

Entonces las especificaciones son:
 * Sobreoscilación: 0%
 * Tiempo de subida: 1.1
 * Tiempo de establecimiento: 1.8
 * Tiempo de pico = 4.5

Vamos a ver qué pasa para T = 1. Para ello utilizaremos //Gz3// definida en el apartado 1.



Entonces las especificaciones son:
 * Sobreoscilación: 3.67 e-14 %
 * Tiempo de pico: 15
 * Tiempo de subida: 1
 * Tiempo de establecimiento: 2

El hecho de que aparezca una tan pequeña sobreoscilación puede deberse a un fallo de cálculo.

Se ve la muy leve mejora del tiempo de subida. Lógico. Si el tiempo de muestra es más grande, llegaremos antes a un valor más cercano al aproximado.

El tiempo de establecimiento se ve levemente perjudicado, pero también tiene sentido. Es más fácil también pasarse del tiempo óptimo.

El tiempo de pico sí que se ve muy perjudicado. Ha tenido una variación muy grande, casi 10 unidades. Esto es malísimo para el sistema.

//**4) Trazar el lugar de las raices discreto con rltool y situar los polos en bucle cerrado para K = 2. Repetir el lugar de las raíces para valores de T = 0.5 y T = 1 y comentar los resultados**// Primero vamos a trazar el lugar de las raices discreto para la función de transferencia G(z) con T=0.1: Al establecer K=2, en la herramienta rtool se obtiene el siguiente lugar de las raices: Si cambiamos el tiempo de muestreo a T=0.5 cambia también la función de transferencia discreta, que ahora es: Trazamos de nuevo el lugar de las raices con rltool (situando la ganancia K=2): Al igual que en el caso anterior al cambiar el tiempo de muestreo a T=1 la función de transferencia pasa a ser: Y el lugar de las raices correspondiente es ahora tal que así: Como conclusión de los resultados obtenidos podemos señalar que a medida de que vamos aumentando el tiempo de muestreo los polos del bucle cerrado del sistema tienden a desplazarse a la izquierda haciendo el sistema mas inestable (e incluso oscilante e inestable).

//**5) Verificar si es posible compensar el sistema con un controlador proporcional en serie (ganancia K) para que el sistema en bucle cerrado y realimentación negativa unitaria cumpla las especificaciones de Sobreoscilación <15% y tiempo de establecimiento t est < 3s y error escalón nulo, graficando la respuesta escalón y dibujando asimismo ambas zonas en el plano z.**//

Utilizando el comando rltool visualizamos el lugar de las raices y añadimos las constricciones:

Obtenemos así la zona donde deben estar los polos para que el sistema cumpla las especificaciones: Ahora solo resta modificar la ganancia y ver si es posible, o no, cumplir las especificaciones: Comprobamos que cuando la salida se aproxima a 1 la sobreoscilacion está muy lejos del máximo permitido. Al reducir la ganancia podemos cumplir la especificaciones de sobreoscilación y tiempo de establecimiento pero el error escalón se hace más grande, con lo cual comprobamos que con un controlador proporcional no se pueden llegar a cumplir todas las especificaciones. //**6) En caso contrario verificar si con la adición de polos y ceros reales es posible cumplir las especificaciones requeridas, expresando el controlador obtenido, graficando el nuevo lugar de las raíces discreto, y comprobando a través de ltiview que se cumplen ciertamente las especificaciones requeridas.**//

Ya teníamos el lugar de las raíces visualizado, con las constricciones añadidas.



Para poder cumplir la condición del error en régimen permanente nulo ante entrada escalón debemos introducir al menos un integrador. ¿Cómo lo hacemos?

Encima del lugar de las raíces vemos un recuadro titulado como //Current Compensator//. Si hacemos clic izquierdo en él se nos abre un cuadro de diálogo.

Queremos probar si añadiendo ceros y polos reales la cosa funciona. Para añadir el integrador añadimos un polo en //z = 1//. Pulsamos Add Real Pole, escribimos 1 y OK.



Este compensador no es suficiente, ya que, aunque cumplamos una de las especificaciones, la otra no podrá cumplirse. La SO se cumple, pero el T.est es mayor.



Casi todo el lugar de las raíces se nos ha desplazado fuera del círculo unidad. Para desplazarlo otra vez dentro vamos a poner un cero cerca del polo de bucle abierto.



Los polos están en //0.72...// así que con que pongamos un cero en 0.775 tendremos un buen resultado. Entonces, los polos, menos uno, quedan dentro del circulo unidad.



El polo que queda fuera puede arrastrarse hacia dentro del círculo unidad, incluso de la restricción de la sobreoscilación, sin que esto produzca inestabilidad al sistema.

Abrimos la ventana de respuesta escalón y vemos cómo ha ido la compensación. Pinchamos en Analysis, Response to Step Command. Activamos las especificaciones.



Las especificaciones pedidas son:
 * Sobreoscilación: 4.48 %
 * Tiempo de establecimiento: 1.7
 * Error en régimen permanente: 0

//**7) Verificar el funcionamiento del controlador obtenido mediante el método de Truxal-Ragazzini en SIMULINK, graficando respuesta escalón.**//

Para verificar el funcionamiento del controlador utilizamos Gd(z): Con la ayuda de SIMULINK hacemos el siguiente esquema: Simulando el sistema (Simulation -> Start) y luego haciendo doble click en el icono de scope podemos ver la salida que es la siguiente: Podemos observar que el sistema sigue a la perfección la entrada escalón unitario (erp=0) y por tanto queda demostrado el funcionamiento del controlador obtenido por método de Truxal-Ragazzini.