Ejercicio 1 Función de Transferencia Pulso

Para hallar la función de transferencia pulso de este sistema tenemos que calcular

$\dfrac{Y(z)}{X(z)}$

Necesitamos conocer la salida en los instantes de muestreo por lo que colocaremos un muestreador imaginario ideal en la salida del sistema y supondremos que los periodos de muestreo de todos los muestreadores son el mismo.

 

De esta forma podemos calcular Y*(s) para a partir de ella calcular Y(z).

Ahora vamos a calcular Y*(s):

Tenemos que  $Y(s)=U*(s)H(s)$

Para obtener Y*(s) tenemos que muestrear Y(s), por lo tanto tenemos que:

$Y(s)=[U*(s)H(s)]*$

A su vez, U*(s) lo calculamos a partir de un muestreo sobre U(s), que lo calcularemos a partir de G(s) y X*(s)

$U*(s)=[U(s)]*=[X*(s)G(s)]*$

Como el muestreo de una señal muestreada con el mismo periodo nos da la misma señal, podemos sacarla del muestreo.

$U*(s)=X*(s)[G(s)]*$

Sustituyendo en Y*(s) tenemos:

$Y*(s)=X*(s)[G(s)]*[H(s)]*$

Haciendo el cambio z=e^{aT} tenemos

$Y(z)=X(z)Z\{G(s)\}Z\{H(s)\}$

Ahora calcularemos G(z) y H(z)

$G(z)=\mathbf{Z}\{G(s)\}=Z\lbrace\dfrac{1}{s+a}\rbrace=\dfrac{z}{z-e^{-aT}}$

$H(z)=\mathbf{Z}\{H(s)\}=Z\{\dfrac{1}{s+b}\}=\dfrac{z}{z-e^{-bT}}$

Ahora podemos calcular la función transferencia pulso.

$\dfrac{Y(z)}{X(z)}=G(z)H(z)$

$\dfrac{Y(z)}{X(z)}=\dfrac{z}{z-e^{-aT}}\dfrac{z}{z-e^{-bT}}=\dfrac{z^2}{(z-e^{-aT})(z-e^{-bT})}$

 

Ejercicio Test de Jury

Vamos a estudiar la estabilidad del siguiente sistema teniendo en cuenta que la función de transferencia de G(z) es

$G(z)=\dfrac{0.3679z+0.2642}{(z-0.3679) (z-1)}$

 

 

Vamos a obtener la función de transferencia del sistema en lazo cerrado.

$\dfrac{Y(z)}{R(z)}$

$E(s)=R(s)-Y(s)$  que de forma discreta es $E(z)=R(z)-Y(z)$

$U(z)=K E(z)$

$Y(z)=U(z)G(z)$

De las dos primeras ecuaciones obtenemos
$U(z)=K(R(z)-Y(z))=KR(z)-KY(z)$

Utilizando U(z) en la ecuación de Y(z) obtenemos

$Y(z)=U(z)G(z)=KR(z)G(z)-KY(z)G(z)$

$Y(z)(1+KG(z))=KR(z)G(z)$

Teniendo finalmente la función de transferencia pulso del sistema

$\dfrac{Y(z)}{R(z)}=\dfrac{KG(z)}{1+KG(z)}=$

La ecuación característica de un sistema es el denominador de la función de transferencia en lazo cerrado. Por lo tanto en nuestro ejemplo la ecuación característica es:

$E.C. \equiv 1+KG(z)=0$

Sustituyendo G(z) por el valor que nos indica el problema tenemos que:

$E.C.\equiv(1+K)\dfrac{0.3679z+0.2642}{(z-0.3679) (z-1)}=0$

Para aplicar el test de Jury tenemos que tener la ecuación característica en forma de polinomio, por lo que vamos a quitar el denominador de la fracción.

$(z-0.3679) (z-1)+K(0.3679z+0.2642)=0$

$z^2-1.3679z+0.3679+0.3679Kz+0.2642K=0$

$z^2+z(0.3679K-1.3679)+0.2642K+0.3679=0$

Una vez colocada en forma de polinomio vamos a aplicar el test de Jury. La matriz del test tendrá un número de filas 2n-3, siendo n el grado de la ecuación característica. En este caso 2n-3=2*2-3=1.

 

Condiciones del test de Jury

Ahora vamos a comprobar las condiciones del test para comprobar la estabilidad según los valores de K y de la E.C.

1) D(1)>0  Tenemos que comprobar que el valor de la E.C. evaluado en z=1 es mayor que 0.

$1+1(0.3679K-1.3679)+0.2642K+0.3679>0$

$(0.3679+0.2642)K>0$

Para que se cumpla esta inecuación, forzosamente K tiene que ser mayor que 0

 

2) $(-1)^n D(-1) >0$  Tenemos que comprobar que -1 elevado al grado de la ecuación multiplicado por la ecuación evaluada en z=-1 es mayor que 0. Como n=2, -1 elevado a 2 es 1, por lo que solo vamos a evaluar E.C. en z=-1.

$-1+(-1)(0.3679K-1.3679)+0.2642K+0.3679>0$

$1-0.3679K+1.3679+0.2642K+0.3679>0$

$2.7358-0.1037K>0$

$-K>\dfrac{-2.7358}{0.1037};K<26.382$

 

3) | a0 | < an  Esta condición la comprobamos utilizando la tabla del test de Jury. Como solo esta formada por una fila, tenemos que evaluar  únicamente una condición.

$0.2642K+0.3679<1$

$K<\dfrac{1-0.2642}{0.3679};K<2.3925$

 

Finalmente acotamos el valor de K para que cumpla las 3 inecuaciones obtenidas.

 

$1) K>0$

$2) K<26.382$

$3) K<2.3925$

 

Las que nos van a condicionar el valor de K son las inecuaciones 1 y 3, que nos van a limitar el valor de K para que sea estable. Por lo tanto K queda acotada de la siguiente forma:

$0< K<2.39$

nota: redondeamos en este caso a la baja para asegurarnos un margen a la hora de asignar un valor concreto a K.

De esta forma el sistema obtenido será estable.

Síntesis Analítica de Controladores Discretos

Enunciado

Se desea controlar la posición angular de un motor de corriente continua que servirá como articulación para un brazo robot. Para ello se ha estudiado la dinámica de dicho motor junto con su etapa de potencia, la reductora empleada para mover el brazo y el potenciómetro que se empleará como sensor de la posición angular del mismo. El modelo obtenido es:

$$Gp(s) = \dfrac{{V_{pot} (s)}}{{U(s)}} = \dfrac{{10}}{{s(s + 1)}},$$

siendo $V_{pot} (s)$ la tensión del potenciómetro de salida (función directa del ángulo en el que está situado el motor), y U(s) la tensión de entrada a la etapa de potencia del motor. El esquema de control empleado es el mostrado en la siguiente figura:

 

siendo $V_{in} (s)$ la entrada de referencia (tensión) del sistema, $C(z)$ el controlador digital y$G_p (s)$ la función de transferencia de la planta a controlar. El conversor analógico-digital se considerará como un muestreador impulsivo, y el digital analógico como un retenedor de orden 0.

Procedimiento de diseño:

1. Obtener la función de transferencia pulso en lazo cerrado: $G_{LC}(z)=\dfrac{{V_{pot}(z)}}{{V_{in} (z)}}.$
2. Simular el comportamiento del sistema con $C(z)=1$, para una entrada de tipo escalón y para una entrada de tipo rampa. Comente los resultados de la simulación teniendo en cuenta el significado de las señales que se representan.
3. Diseñar el controlador $C(z)$ mediante el Lugar Geométrico de las raíces para que el tiempo de asentamiento sea, según el criterio del 2%, inferior a 2 segundos, y el sobreimpulso máximo no supere el 5%. Para ello obtenga en primer lugar los polos que debería tener un sistema en tiempo continuo para cumplir las especificaciones, y calcule a partir de éstos los de un sistema en tiempo discreto. En último lugar, utilice el L.G.R. para situar los polos de lazo cerrado del sistema aproximadamente en el punto calculado.
4. Simular el comportamiento del sistema con el controlador diseñado para una entrada de tipo escalón y para una entrada de tipo rampa. Comente los resultados de la simulación.
5. Obtenga el algoritmo de control mediante el método directo para el controlador diseñado y programe mediante pseudo-código dicho algoritmo.

Solución

1) Obtener la función de transferencia pulso del sistema en lazo cerrado
Para poder trabajar con la transformada Z, colocamos un muestreador virtual a la salida del sistema, quedando de la siguiente manera:

El conversor analógico-digital es un muestreador impulsivo y el conversor digital-analógico es un retenedor de orden cero con la siguiente función de transferencia:

$G_{h0} (s) = \dfrac{1 - e^{Ts}}{s}$

Para el cálculo de la función de transferencia del sistema, agrupamos en un solo bloque el retenedor de orden cero y la planta, quedando el diagrama de la siguiente manera:

Realizamos los cálculos para despejar la función de transferencia:

$\begin{array}{l} V_{pot} (s) = G_1 (s)U^* (s) \\
U^* (s) = C(z)E^* (s) \\
E^* (s) = \left[ {V_{in} (s) - V_{pot} (s)} \right]^* \\
\end{array}$

Despejamos:

$\begin{array}{l} E^* (s) = \left[ {V_{in} (s) - V_{pot} (s)} \right]^* = V_{in} ^* (s) - V_{pot} ^* (s) \\
V_{pot} ^* (s) = \left[ {G_1 (s)U^* (s)} \right]^* = G_1 ^* (s)U^* (s) = G_1 ^* (s)C(z)E^* (s) = \\
= G_1 ^* (s)C(z)\left[ {V_{in} ^* (s) - V_{pot} ^* (s)} \right] \\
\end{array}$

Finalmente, la última expresión queda de la siguiente manera:

$\begin{array}{l}
V_{pot} ^* (s)\left[ {1 + G_1 ^* (s)C(z)} \right] = G_1 ^* (s)C(z)V_{in} ^* (s) \\
\dfrac{{V_{pot} ^* (s)}}{{V_{in} ^* (s)}} = \dfrac{{G_1 ^* (s)C(z)}}{{1 + G_1 ^* (s)C(z)}} \\
\end{array}$

Expresando todo en Z, la función de transferencia queda:

$G_{LC} (z) = \dfrac{{V_{pot} (z)}}{{V_{in} (z)}} = \dfrac{{G_1 (z)C(z)}}{{1 + G_1 (z)C(z)}}$

2) Simular el comportamiento del sistema con $C(z) = 1$, para una entrada de tipo escalón y para una entrada de tipo rampa. Comente los resultados de la simulación teniendo en cuenta el significado de las señales que se representan.

Como el controlador no tiene ningún efecto, omitimos los conversores A/D.

Simulamos la entrada escalón: (Cuidado, esta simulación puede no ser correcta)

H = tf(10,[1 1 10]);
T = 0:0.1:15;
Y = step(H,T);
X = ones(1,length(T));
hold on;
plot(T,Y,'b');
plot(T,X,'r');

Se observa que el sistema es subamortiguado, y que tarda en estabilizarse aproximadamente unos 7,5s. El sobreimpulso es excesivo, ya que es aproximadamente del 60%. También se observa que el error de estado estacionario de posición es nulo, puesto que es sistema es de tipo 1.

Simulamos la entrada rampa de pendiente 1: (Cuidado, esta simulación puede no ser correcta)

H = tf(10,[1 1 10]);
T = 0:0.1:10;
Y = lsim(H,T,T);
hold on;
plot(T,Y,'b');
plot(T,T,'r');

Es esta respuesta se obsera obviamente la misma dinámica que en la respuesta anterior, es decir, la de un sistema subamortiguado con un gran sobreimpulso. Sin embargo, el error de velocidad no es nulo como lo era el de posición, ya que es aproximadamente del 10%.

Para medir el error mediante MATLAB, puede teclearse:

error_velocidad=Y(end)-T(end)

3) Diseñar el controlador $C(z)$ mediante el Lugar Geométrico de las Raíces, para que el tiempo de asentamiento sea, según el criterio del 2%, inferior a 2s, y el sobreimpulso máximo no supere el 5%. Para ello obtenga en primer lugar los polos que debería tener un sistema en tiempo continuo para cumplir las especificaciones, y calcule a partir de éstos los de un sistema discreto. En último lugar, utilice el L.G.R. para situar los polos de lazo cerrado del sistema aprocimadamente en el punto calculado.

Calculamos los polos teóricos de tiempo continuo que cumplan las especificaciones:

$\delta = \sqrt {\dfrac{{\ln ^2 (M_p )}}{{\ln ^2 (M_p ) + \pi }}}$

$\omega _n = \dfrac{4}{{\delta T_s }}$

$p_{1,2} (s) = - \delta \omega _n \pm j\omega _n \sqrt {1 - \delta ^2 }$

Resolviendo, para T=0.1:
delta = sqrt((log(Mp)^2) / (pi^2 + log(Mp)^2))
wn = 4 / (delta * ts)
p1 = -delta * wn + j * wn * sqrt(1 - delta^2)
p2 = -delta * wn - j * wn * sqrt(1 - delta^2)

Obtenemos los polos en tiempo continuo que cumplen las especificaciones:
p1 = -2.0000 + 2.0974i
p2 = -2.0000 - 2.0974i
Para obtener los valores en tiempo discreto solo tenemos que realizar:
p1z = exp(T * p1)
p2z = exp(T * p2)
Obtenemos los polos en tiempo discreto:
p1 = 0.8008 + 0.1705i
p2 = 0.8008 - 0.1705i

Observando en la herramienta sisotool la estabilidad del sistema nos damos cuenta que a partir de K=0.16 el sistema se vuelve inestable:

 

 

Obtenemos el siguiente controlador:
$C(z) = 0.16 \cdot \dfrac{1}{{1 + 0.58z^{ - 1} + (0.38z^{ - 1} )^2 }}$
4) Simular el comportamiento del sistema con el controlador diseñado para una entrada de tipo escalón y para una entrada de tipo rampa. Comente los resultados de la simulación.

Simulamos la entrada escalón:

clear all, close all, clc;
T = 0.1;
% Constantes
G = tf(10, [1 1 0]);
C = tf([0.16 0 0], [1 0.58 0.1444],T);
Gz = c2d(G,T);
GC = Gz*C;
T = 0:0.1:15;
GLC1 = GC/(1+GC);
Y = step(GLC1,T);
X = ones(1,length(T));
hold on;
plot(T,Y,'b');
plot(T,X,'r');

 

 

Vemos que el sistema se va estabilizando en torno al valor 1, con alguna pequeña oscilación, pero respetando los valores indicados.

Para una entrada tipo rampa de pendiente 1:

clear all, close all, clc;
T = 0.1;
% Constantes
G = tf(10, [1 1 0]);
C = tf([0.16 0 0], [1 0.58 0.1444],T);
Gz = c2d(G,T);
GC = Gz*C;
T = 0:0.1:10;
GLC1 = GC/(1+GC);
Y = lsim(GLC1,T,T);
hold on;
plot(T,Y,'b');
plot(T,T,'r');

 

 

 

Se observa un error de estado estacionario de aproximadamente valor 1. Este error es debido a que el sistema con el controlador es de tipo 1.

5) Obtenga el algoritmo de control mediante el método directo para el controlador diseñado y programe mediante pseudo-código dicho algoritmo.
$C(z) = \dfrac{{Y(z)}}{{X(z)}} = \dfrac{1}{{1 + 0.58z^{ - 1} + 0.1444z^{ - 2} }}$

$Y(z) \cdot (1 + 0.58z^{ - 1} + 0.1444z^{ - 2} ) = X(z)$

$Y(z) + 0.58z^{ - 1} Y(z) + 0.1444z^{ - 2} Y(z) = X(z)$

Pasamos a ecuaciones diferenciales:

$Y[n] + 0.58 \cdot Y[n - 1] + 0.1444 \cdot Y[n - 2] = X[n]$

$Y[n] = X[n] - 0.58 \cdot Y[n - 1] - 0.1444 \cdot Y[n - 2]$


Este pseudo-código está mal, hay que programar el algoritmo de control, es decir C(z):

pseudo-código:
Funcion Y=Y(n)
Si n=0 entonces retorna Y=0.144
Si n=1 entonces retorna Y=0.58
Si n>1 entonces retorna Y=X(n)-0.58*Y(n-1)-0.144*Y(n-2)