Velocidad del Motor de CC : Modelación

El motor de CC es un actuador común en control sistemas. Provee movimiento rotatorio directamente y, acoplado con ruedas dentadas o poleas y cables, puede proveer movimiento transicional. El circuito eléctrico de la armadura y el diagrama de cuerpo libre del rotor se muestran en la siguiente figura:

Para este ejemplo, asumimos los valores siguientes para los parametros físicos. Estos valores se derivaron experimentalmente de un motor real del laboratorio de control para alumnos de grado del Carnegie Mellon.
* momento de inercia del rotor (J) = 0.01 kg.m^2/s^2
* coeficiente de amortiguamiento del sistema mecanico (b) = 0.1 Nms
* constante de fuerza electromotriz (K=Ke=Kt) = 0.01 Nm/Amp
* resistencia eléctrica (R) = 1 ohm

* inductancia eléctrica (L) = 0.5 H
* entrada (V): Fuente de Tensión
* salida (theta): posición del eje
* el rotor y eje se consideran rígidos
El torque del motor, T, se relaciona con la corriente de armadura, i, por un factor constante Kt. La fuerza contraelectromotriz (emf), e, se relaciona con la velocidad de rotación mediante las siguientes ecuaciones

En unidades del sistema internacional SI (las que usaremos), la Kt (constante de armadura) es igual a Ke (constante del motor).
De la figura de arriba podemos escribir las siguientes ecuaciones basadas en la ley de Newton combinado con la ley de Kirchhoff:
:

1. Función de Transferencia
Usando Transformadas deLaplace, las ecuaciones del modelo de arriba pueden expresarse en términos de s.

Eliminando I(s) podemos obtener la siguiente función de transferencia a lazo abierto , donde la velocidad de rotación es la salida y el voltaje es la entrada.

2. Espacio de Estado
En la forma espacio de estado, las ecuaciones de arriba pueden expresarse escogiendo la velocidad de rotación y corriente eléctrica como las variables de estado y la tensión como una entrada. La salida se elige que sea la velocidad de rotación.

Requerimientos de diseño
Primero, el motor sin compensar puede rotar solo a 0.1 rad/seg. con una tensión de entrada de 1 Volt (esto se demostrara luego cuando se simule la respuesta a lazo abierto). Como el requerimiento mas basico de un motor es que debe rotar a la velocidad deseada, el error de estado estacionario de la velocidad del motor debe ser menor que 1%. El otro requerimiento de performance es que el motor debe acelerarse hasta su velocidad de estado estacionario apenas se encienda. En este caso, queremos tener un tiempo de establecimiento de 2 segundos. Como una velocidad mayor que la referencia podría dañar el equipo, queremos tener un sobrepico menor que 5%.
Si simulamos la entrada de referencia (r) con una entrada escalón unitario, entonces la salida velocidad del motor debería tener:
* Tiempo de establecimiento menor que 2 segundos
* Sobrepico menor que 5%
* error de estado estacionario menor que 1%
Representación en Matlab y respuesta a lazo abierto
1. Función de Transferencia
Podemos representar lafunción de transferencia anterior en Matlab definiendo las matrices numerador y denominador como sigue:

Cree un nuevo archivo-m e ingrese los siguientes comandos:
J=0.01;
b=0.1;
K=0.01;
R=1;
L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];
Ahora veamos qué hace el sistema original a lazo abierto. Agregue los siguientes comandos al final del archivo-m y ejecútelo en la ventana de comandos del Matlab:
step(num,den,0:0.1:3)
title('Respuesta al Escalón del sistema a lazo abierto')
Debería obtenerse la figura siguiente:

De la figura vemos que cuando se aplica 1 volt al sistema, el motor puede lograr solo una velocidad maxima de 0.1 rad/seg., diez veces menor que la velocidad deseada. Ademas, el motor tiene 3 segundos para alcanzar su velocidad de estado estacionario; esto no satisface los 2 segundos del criterio de tiempo de establecimiento.
2. Espacio de Estado
Podemos también representar el sistema usando las ecuaciones de espacio de estado. Pruebe los siguientes comandos en un nuevo archivo-m.
J=0.01;
b=0.1;
K=0.01;
R=1;
L=0.5;
A=[-b/J K/J
-K/L -R/L];
B=[0
1/L];
C=[1 0];
D=0;

step(A, B, C, D)

Método de Diseño PID for Velocidad del Motor de CCControl
Control proporcional
control PID
Sintonización de ganancias
A partir del problema principal, las ecuaciones dinamicas y la función transferencia a lazo abierto del Motor de CC son:

y el esquema del sistema se ve:

Para las condiciones originales del problema y la derivación de las ecuaciones de arriba, refiérase por favor a Modelación de un Motor de CC .
Los criterios de diseño con una entrada escalón de 1 rad/seg. son:
* Tiempo de establecimiento menor que 2 segundos
* Sobrepico menor que 5%
* Steady-stage error menor que 1%
Diseñemos ahora un controlador PID adicionémoslo al sistema. Cree primero un archivo-m nuevo y tipee los siguientes comandos (pinche en Modelación del motor de cc por detalles sobre estos comandos).
J=0.01;
b=0.1;
K=0.01;
R=1;
L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];
Recordemos que la función de transferencia para un controlador PID es:

Control proporcional
Tratemos de usar primero un controlador proporcional con una ganancia de 100. Agregue el siguiente código al final de su archivo-m:
Kp=100;
numa=Kp*num;
dena=den;
Para hallar la función de transferencia a lazo cerrado , usamos el comando cloop . Agregue la siguiente línea a su archivo-m:
[numac,denac]=cloop(numa,dena);
Note que numac y denac son el numerador y el denominador de la función de transferencia a lazo cerrado general.Ahora veamos cómo se ve la respuesta al escalón, agregue lo siguiente al final de su archivo-m, y ejecútelo en la ventana de comandos:
t=0:0.01:5;
step(numac,denac,t)
title('Respuesta al escalón con Control Proporcional')
Debería obtenerse la figura siguiente:

control PID
De la figura de arriba vemos que tanto el error de estado estacionario cuanto el sobrepico son muy grandes. Recordemos del tutorial PID que incorporando un término integral se eliminara el error de estado estacionario y un término derivativo reducira el sobrepico. Probemos un controlador PID con Ki y Kd pequeños. Modifique su archivo-m de manera que se sea ve lo siguiente. Luego de correr este nuevo archivo-m se da la figura siguiente.

J=0.01;
b=0.1;
K=0.01;
R=1;
L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];

Kp=100;
Ki=1;
Kd=1;
numc=[Kd, Kp, Ki];
denc=[1 0];
numa=conv(num,numc);
dena=conv(den,denc);
[numac,denac]=cloop(numa,dena);
step(numac,denac)
title('Control PID con pequeño Ki and Kd')

Sintonización de ganancias
Ahora el tiempo de establecimiento es muy largo. Incrementemos Ki para reducir el tiempo de establecimiento. Regrese a su archivo-m y cambie Ki a 200. Vuelva a ejecutar elarchivo y debería obtenerse un grafico com éste:

Ahora vemos que la respuesta es mucho mas rapida que antes, pero el Ki grande ha empeorado la respuesta transitoria (gran sobrepico). Incrementemos Kd para reducir el sobrepico. Vuelva al archivo-m y cambie Kd a 10. Vuelva a correrlo para obtener este grafico:

Entonces sabemos que si usamos un controlador PID con
Kp=100,
Ki=200,
Kd=10,
todos nuestros requerimientos de diseño seran satisfechos.
Método de Diseño por Lugar de Raíces para Control de Velocidad del Motor de CC
Dibujo del lugar de raíces a lazo abierto
Calculo de la ganancia mediante el comando rlocfind
Adición de un controlador en atraso
Dibujo de la respuesta a lazo cerrado
A partir del problema principal, las ecuaciones dinamicas y la función de transferencia a lazo abierto de la Velocidad del Motor de CC son:

y el esquema del sistema se ve:

Para las condiciones originales del problema y la derivación de las ecuaciones de arriba, refiérase por favor a Modelación de un Motor de CC.
Con una referencia escalón de 1 rad/seg., los criterios de diseño son:
* Tiempo de establecimiento menor que 2 segundos
* Sobrepico menor que 5%
* error de estado estacionario menor que 1%
Diseñemos ahora un controlador mediante el método del lugar de raíces.
Cree un nuevo archivo-m y tipee los siguientes comandos (refiérase a problema original para detalles de cómo obtener estos comandos).
J=0.01;
b=0.1;
K=0.01;
R=1;L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];
Dibujo del lugar de raíces a lazo abierto
La idea principal del diseño por lugar de raíces es para hallar la respuesta a lazo cerrado a partir del diagrama del lugar de raíces a lazo abierto. Entonces agregando ceros y/o polos a la planta original, puede modificarse la respuesta a lazo cerrado. Veamos primero el lugar de raíces de la planta. Agregue los siguientes comandos al final de su archivo-m.
rlocus(num,den)
sgrid(.8,0)
sigrid(2.3)
title('Lugar de Raíces sin a controlador')
El comando sigrid es la función definida por el usuario. Debe copiar el archivo sigrid.m a su directorio para poder usarlo. Para mas información sobre cómo usar las funciones, refiérase a funciones.
Los dos argumentos en el comando sgrid son el coeficiente de amortiguamiento (zeta) (0.8 corresponde a un sobrepico del 5%), y la frecuencia natural (Wn) (= 0 corresponde a no tiempo de subida ) respectivamente. El único argumento en el comando sigrid es el término sigma (4.6/2 segundos = 2.3). Luego de salvar el archivo sigma.m en su directorio, corra el archivo-m de arriba en la ventana de comandos. Debería obtener el diagrama del lugar de raíces que se muestra abajo:

Calculo de la ganancia mediante el comando rlocfind
Si recuerda, necesitamos que el tiempo de establecimiento y el sobrepico sean lo mas pequeño posible. Amortiguamientos grandes se corresponden con puntos cercanos al eje real en el lugar de raíces.Una respuesta rapida se corresponde con puntos muy a la izquierda del eje imaginario en el lugar de raíces. Para hallar la ganancia correspondiente a un punto en el lugar de raíces, podemos usar el comando rlocfind . Podemos hallar la ganancia y dibujar la respuesta al escalón usando esta ganancia todo al mismo tiempo. Para hacerlo, entre los siguientes comandos al final de su archivo-m y ejecútelo nuevamente.
[k,poles] = rlocfind(num,den)
[numc,denc]=cloop(k*num,den,-1);
t=0:0.01:3;
step(numc,denc,t)
title('Respuesta al escalón con ganancia')
Vaya a la figura y elija un punto en el lugar de raíces a medio camino entre el eje real y el requerimiento de amortiguamiento , digamos en -6+2.5i. Matlab debería devolver la salida similar la siguiente.

selected_point =
-5.9596 + 2.0513i

k =
10.0934

poles =
-6.0000 + 2.0511i
-6.0000 - 2.0511i
Note que los valores devueltos en la ventana de comandos del Matlab podrían no ser exactamente los mismos, pero debería al menos tener el mismo orden de magnitud. Debería obtener también la figura siguiente:

como puede ver, el sistema esta sobreamortiguado y el tiempo de establecimiento es alrededor de un segundo, Así que se satisfacen los requerimientos de sobrepico y tiempo de establecimiento. El único problema que podemos ver en esta figuraes el error de estado estacionario de acerca de 50%. Si incrementamos la ganancia para reducir el error de estado estacionario, el sobrepico se vuelve muy grande (pruébelo). Necesitamos agregar un controlador en atraso para reducir el error de estado estacionario.
Adición de un controlador en atraso
De la figura vemos que tenemos un lugar de raíces muy simple. Los criterios de amortiguamiento y tiempo de establecimiento se han logrado con el controlador proporcional. El error de estado estacionario es el único criterio no satisfecho con el controlador proporcional. Un compensador en atraso puede reducir el error de estado estacionario. Mediante esto, podríamos sin embargo incrementar el tiempo de establecimiento. Primero pruebe el siguiente controlador en atraso:

Esto puede hacerse cambiando su archivo-m al siguiente:
J=0.01;
b=0.1;
K=0.01;
R=1;
L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];

z1=1;
p1=0.01;
numa = [1 z1];
dena = [1 p1];
numb=conv(num,numa);
denb=conv(den,dena);

rlocus(numb,denb)
sgrid(.8,0)
sigrid(2.3)
title('Lugar de Raíces con un controlador en atraso')
numa y dena son el numerador y denominador del controlador, y numb y denb son el numerador y denominador de la función de transferencia a lazoabierto entrada-salida .
Debería obtenerse el siguiente lugar de raíces, el cual luce muy parecido al original:

Dibujo de la respuesta a lazo cerrado
Cerremos ahora el lazo y veamos la respuesta al escalón a lazo cerrado. Entre el siguiente código al final de su archivo-m:
[k,poles]=rlocfind(numb,denb)
[numc,denc]=cloop(k*numb,denb,-1);
t=0:0.01:3;
step(numc,denc,t)
title('Respuesta al escalón con un controlador en atraso')
Vuelva a ejecutar este archivo-m en la ventana de comandos del Matlab. Cuando se le pida elegir un punto, escoja uno que esté cerca del requerimiento de amortiguamiento (línea puntedada diagonal). Debería obtenerse una figura similar a la siguiente:

Su ganancia debería ser cercana a 20. Como puede ver la respuesta no es muy satisfactoria. Debería notarse también que a pesar que se seleccionó la ganancia correlacionada con una posición cerca del criterio de amortiguamiento, el sobrepico no se acerca al 5%. Esto se debe al efecto del controlador en atraso . (su polo es mas lento). Lo que significa que podemos ir fuera de las líneas punteadas que representan el límite, y obtener las mas altas ganancias sin preocuparnos por el sobrepico . Vuelva a ejecutar su archivo-m, ponga la ganancia apenas por encima de la línea blanca punteada. Siga intentando hasta que obtenga una respuesta satisfactoria. Debería verse similar al siguiente (usamos una ganancia de aproximadamente 50):

El error de estado estacionario es menor que el 1%, ylos requerimientos de tiempo de establecimiento y sobrepico se han satisfecho. Como puede ver, los procesos de diseño por lugar de raíces tienen mucho de procedimientos de prueba y error. Esto es porque es facil dibujar el lugar de raíces, escoger la ganancia, y dibujar la respuesta en un sólo paso. Si no fuéramos capaces de obtener una respuesta satisfactoria eligiendo las ganancias, hubiésemos podido probar un controlador en atraso diferente, o aún adicionado a un controlador en adelanto.
Diseño Frecuencial para el Control de Velocidad del Motor de CC
Dibujo del diagrama de Bode original
Agregado de ganancia proporcional
Dibujo de la respuesta a lazo cerrado
Adición de un controlador en atraso
A partir del problema principal, las ecuaciones dinamicas y la función de transferencia a lazo abierto de la Velocidad del Motor de CC son:

y el esquema del sistema se ve:

Para las condiciones originales del problema y la derivación de las ecuaciones anteriores, refiérase por favor a Modelación de un Motor de CC .
Los criterios de diseño con una entrada escalón de 1 rad/seg., son:
* Tiempo de establecimiento menor que 2 segundos
* Sobrepico menor que 5%
* Error de estado estacionario menor que 1%
Cree un nuevo archivo-m y tipee los siguientes comandos (pinche en problema original para detalles de cómo obtener estos comandos).
J=0.01;
b=0.1;
K=0.01;
R=1;
L=0.5;
num=K;
den=[(J*L) ((J*R)+(L*b))((b*R)+K^2)];
Dibujo del diagrama de Bode original
La idea principal del diseño basado en la frecuencia es usar el diagrama de Bode del función de transferencia a lazo abierto para estimar la respuesta a lazo cerrado. Agregando un controlador al sistema se cambia el diagrama de Bode a lazo abierto, por lo tanto también se cambia la respuesta a lazo cerrado. Primero grafiquemos el Diagrama de Bode para la función de transferencia a lazo abierto original . Agregue el siguiente código al final de su archivo-m, y entonces ejecútelo en la ventana de comandos del Matlab.
bode(num,den)
Debería obtenerse el siguiente diagrama de Bode:

Agregado de ganancia proporcional
Del diagrama de Bode anterior, vemos que el margen de fase puede ser mayor que aproximadamente 60 grados si w es menor que 10 rad/seg.. Agreguemos ganancia al sistema de modo que el ancho de banda sea 10 rad/seg., lo que nos dara un margen de fase de acerca de 60 grados. Para hallar la ganancia a 10 rad/seg., puede intentar de extraerlo del diagrama de Bode (parece ser apenas mas que -40 dB, o 0.01 en magnitud). El comando bode, invocado con argumentos a izquierda, puede también ser usado para proporcionar la magnitud exacta:
[mag,fase,w] = bode(num,den,10)
mag =

0.0139
para tener una ganancia de 1 a 10 rad/seg., multiplique el numerador por 1/0.0139 o aproximadamente 72.
num = 70*num
y re-ejecute su archivo-m. Debería tener el siguiente diagrama de Bode:

Dibujo de larespuesta a lazo cerrado
De la figura de arriba vemos que el margen de fase ahora es bastante grande. Veamos cómo luce la respuesta a lazo cerrado. Agregue un signo % al fina de los comandos bode y también el siguiente código al final de su archivo-m:
[numc,denc]=cloop(num, den, -1);
t=0:0.01:10;
step(numc,denc,t)
Vera la figura siguiente:

El tiempo de establecimiento es bastante rapido, pero el sobrepico y el error de estado estacionario son muy altos. El sobrepico puede ser reducido bajando un poco la ganancia para obtener un margen de fase mas alto, pero esto puede causar que el error de estado estacionario se incrementara. Probablemente se necesite un controlador en atraso.
Adición de un controlador en atraso
Podemos agregar un controlador en atraso para reducir el error de estado estacionario. Al mismo tiempo, debería tratarse de achicar el sobrepico reduciendo la ganancia. Reduzcamos la ganancia a 50, y probemos un controlador en atraso de

el cual debe reducir el error de estado estacionario en un factor de 1/0.01 = 100 (pero pudiera incrementar el tiempo de establecimiento). Regrese y cambie su archivo-m de manera que se sea ve lo siguiente:
num=K;
den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];
num=50*K;

z=1;
p=0.1;
numa=[1 z];
dena=[1 p];
numb=conv(num,numa);
denb=conv(den,dena);bode(numb,denb)
Vuelva a ejecutar el archivo y obtendra este grafico:

El margen de fase aparenta bueno. El error de estado estacionario predicho estaría cerca de 1/40dB o 1%, como era lo deseado. Cierre el lazo y observe la respuesta al escalón . Agregue las siguientes líneas de código al final de su archivo-m y vuelva a ejecutarlo.
[numc,denc]=cloop(numb, denb, -1);
t=0:0.01:10;
step(numc,denc,t)

Ahora tiene a respuesta al escalón que satisface los requerimientos de diseño. El error de estado estacionario es menor que 1%, el sobrepico es alrededor de 5%, y el tiempo de establecimiento es alrededor de 2 segundos.
Control Digital de Velocidad
del Motor de CC con Control PID
Conversión de Continuo a Discreto
Controlador PID
En esta pagina, consideremos la versión digital del problema de control de velocidad de un motor de CC . Puede obtenerse un modelo digital del motor de CC a partir de la conversión del modelo analógico , como describiremos. El controlador para este ejemplo se diseñara por un método PID.
De la pagina Modelación: un Motor de CC, la función de transferencia a lazo abierto para la velocidad del motor se derivó como:

Donde:
*resistencia eléctrica (R) = 1 ohm
*inductancia eléctrica (L) = 0.5 H
*constante de fuerza electromotriz (Ke=Kt) = 0.01 Nm/Amp
*momento de inercia del rotor (J) = 0.01 kg*m^2/s^2
*coeficiente de amortiguamiento del sistema mecanico (b) = 0.1 Nms
*entrada (V): Fuente de Tensión
*salida (theta punto): Velocidad de rotación
*Se asume queel rotor y el eje son rígidos
Los requerimientos de diseño para entrada escalón de 1 rad/seg. son
* Tiempo de establecimiento: Menor que 2 segundos
* Sobrepico: Menor que el 5%
* Error de estado estacionario: Menor que el 1%
Conversión de Continuo a Discreto
El primer paso en el diseño de un sistema de control discreto es convertir la función de transferencia continua a una función de transferencia discreta. El comando Matlab c2dm lo hara para ud.. c2dm requiere los siguiente cuatro argumentos: el polinomio numerador (num), el polinomio denominador (den), el tiempo de muestreo (Ts) y el tipo de mantenedor. En este ejemplo, usaremos el mantenedor de orden cero ('zoh').
Del requerimiento de diseño, hagamos el tiempo de muestreo, Ts igual a 0.12 segundos, lo cual es 1/10 de la constante de tiempo de un sistema con un tiempo de establecimiento de 2 segundos. Creemos un nuevo archivo-m e ingrese los siguientes comandos:
R=1;
L=0.5;
Kt=0.01;
J=0.01;
b=0.1;

num = Kt;
den = [(J*L) (J*R)+(L*b) (R*b)+(Kt^2)];

Ts = 0.12;
[numz,denz] = c2dm(num,den,Ts,'zoh')
Luego de ejecutado este archivo-m debe devolver lo siguiente:

numz =

0 0.0092 0.0057


denz =

1.0000 -1.0877 0.2369De estas matrices, la función de transferencia discreta puede escribirse como:

Primero, nos gustaría vea cómo se ve la respuesta del sistema a lazo cerrado sin ningún control. Si aprecia la matriz numz arriba, ésta tiene un cero adicional al principio, tenemos que librarnos de él antes de cerrar el lazo con el comando Matlab cloop. Agregue el siguiente código al final de su archivo-m:
numz = [numz(2) numz(3)];
[numz_cl,denz_cl] = cloop(numz,denz);
Luego de hacerlo, observemos cómo se ve la respuesta a lazo cerrado al escalón . El comando dstep generara el vector de señales de salida discreto y el comando stairs las conectara . Pinche aquí para mas información. Agregue el siguiente código de Matlab al final del archivo-m anterior y ejecútelo nuevamente.
[x1] = dstep(numz_cl,denz_cl,101);
t=0:0.12:12;
stairs(t,x1)
xlabel('Tiempo (segundos)')
ylabel('Velocidad (rad/s)')
title('Respuesta tipo escalera :Original')
Debería verse la figura siguiente:

Controlador PID
Recordemos que la función de transferencia de tiempo continuo para un controlador PID es:

Existen varias maneras para el mapeo del plano s al plano z . La mas precisa es . No podemos obtener la función de transferencia del filtro PID de este modo porque la función de transferencia de tiempo discreto tendría mas ceros que polos, lo cual no es realizable. En su lugar usaremos la transformación bilineal , definida como sigue:

Por lo quepodemos derivar el controlador PID discreto con el mapeo por la transformación bilineal. Para una derivación mas detallada de un controlador PID discreto , vea Controlador PID Discreto . Equivalentemente, el comando c2dm del Matlab le ayudara para convertir el compensador PID de tiempo continuo al compensador PID de tiempo discreto usando el método 'tustin' en este caso. El método 'tustin' usara la aproximación bilineal para la conversión a tiempo discreto de la derivada. De acuerdo con la pagina Método de Diseño del PID para el Motor de CC , Kp = 100, Ki = 200 y Kd = 10 satisfacen el requerimiento de diseño. Se usara todas las ganancias en este ejemplo. Ahora agregue los siguientes comandos Matlab a su archivo-m anterior y ejecútelo nuevamente en la ventana del Matlab.
% Controlador PID Discreto con aproximación bilineal
Kp = 100;
Ki = 200;
Kd = 10;

[dencz,numcz]=c2dm([1 0],[Kd Kp Ki],Ts,'tustin');
Note que el numerador y el denominador en el comando c2dm se invirtieron arriba. La razón es porque función de transferencia PID no es propia. Matlab no permite esto. Intercambiando el numerador y el denominador el comando c2dm puede ser engañado para devolver la respuesta correcta. Veamos si la performance de la respuesta a lazo cerrado con el compensador PID satisface los requerimientos de diseño. Agregue ahora el código siguiente al final de su archivo-m y ejecútelo nuevamente. Debería obtenerse la siguiente respuesta tipo escalera a lazocerrado.
numaz = conv(numz,numcz);
denaz = conv(denz,dencz);
[numaz_cl,denaz_cl] = cloop(numaz,denaz);

[x2] = dstep(numaz_cl,denaz_cl,101);
t=0:0.12:12;
stairs(t,x2)
xlabel('Tiempo (segundos)')
ylabel('Velocidad (rad/s)')
title('Respuesta tipo escalera: con controlador PID ')

como puede ver de la figura de arriba, la respuesta del sistema a lazo cerrado es inestable. Por lo tanto debe haber algo mal en el sistema compensado. Así que echemos un vistazo al root locus del sistema compensado. Agreguemos los siguientes comandos Matlab al final de su archivo-m y ejecútelo nuevamente.
rlocus(numaz,denaz)
title('Lugar de Raíces del Sistema Compensado')

De este grafico del lugar de raíces, vemos que el denominador del controlador PID tiene un polo en -1 en el plano z . Sabemos que si un polo de un sistema esta fuera del círculo unitario, el sistema sera inestable. Este sistema compensado siempre sera instable para cualquier ganancia positiva porque habra un número par de polos y ceros a la derecha del polo en -1. Por lo tanto ese polo siempre se movera a la izquierda y afuera del círculo unitario. El polo en -1 viene del compensador, y podemos cambiar su ubicación cambiando el diseño del compensador. Lo elegimos para cancelar el cero en -0.62. Esto hara el sistema estable para al menos algunas ganancias. Ademas podemos elegir una ganancia apropiadadel diagrama del lugar de raíces para satisfacer los requerimientos de diseño usando rlocfind.Ingrese el siguiente código de Matlab a su archivo-m.
dencz = conv([1 -1],[1.6 1])
numaz = conv(numz,numcz);
denaz = conv(denz,dencz);

rlocus(numaz,denaz)
title('Lugar de Raíces del Sistema Compensado');
[K,polos] = rlocfind(numaz,denaz)
[numaz_cl,denaz_cl] = cloop(K*numaz,denaz);

[x3] = dstep(numaz_cl,denaz_cl,101);
t=0:0.12:12;
stairs(t,x3)
xlabel('Tiempo (segundos)')
ylabel('Velocidad (rad/s)')
title('Respuesta tipo escalera:con controlador PID ')
El nuevo dencz tendra un polo en -0.625 en lugar de -1, que casi cancela al cero del sistema no compensado. En la ventana del Matlab, debería ver el comando que le solicita elegir el punto en el grafico de lugar de raíces. Debe pinchar en el grafico como sigue:

Entonces Matlab devolvera la ganancia apropiada y los polos compensados correspondientes, y se ploteara la respuesta a lazo cerrado compensada como sigue.

La figura muestra que el tiempo de establecimiento es menor que 2 segundos y el sobrepico porcentual es alrededor del 3%. Ademas, el error de estado estacionario es cero. Y también, la ganancia, K, del lugar de raíces es 0.2425 lo cual es razonable. Por lo tanto esta respuesta satisface todos los requerimientos de diseño.