Blackat tu última funcion G(s) no existe fisicamente, debería estar invertida.
saludos,
Que tal, oigan esta muy interesante y muy buena la informacion que dan aqui, la vdd me esta ayudando mucho a darme idea de como es esto del PID, solo que mi gran y principal problema es que en la escuela me pidieron hacer esto del PID discreto con un microcontrolador (ya sea pic o avr), pero sin enseñarme absolutamente nada al respecto de que es y como se hace un PID, todo me lo dejaron a investigacion, y sii ya he buscado por toodos lados y he encontrado varias cosas, pero lo que sigo sin encontrar o sin entender, es exactamente cuales son las variables que necesito conocer de mi "caja negra" conformada en mi caso por una hielera de unicel, un foco de 40watts y un lm35.. y bueno por el momento esa es mi principal duda, espero pueda alguien auxiliarme porque me lo encargaron la semana pasada y lo tengo que entregar a finales de esta semana :(.. estoy practicamente desesperado. de ante mano gracias.
Y ¿Como se implementa una ecuacion de un controlador en C?Aquí (http://www.todopic.com.ar/foros/index.php?topic=26780.0) hay uno implementado en C.
Un saludo.
Ese promediador parece que no retrasa la señal, la voy a probar para filtrar la velocidad de mi motor ya que uso un periodo de muestreo de 1mseg.
Nocturno te sale un bajo porcentaje de precisión, porque tu curva calculada no se aproxima bien en muchas zonas a la curva real, segun se ve tu tiempo de establecimiento es 0.6seg.
saludos.
Una manera más simple de sacar la función de tu planta es hacerla al ojo :?, la respuesta de velocidad de un motor dc tiene la forma de una funcion de primer orden es decir sin sobreimpulso. Entonces tendrías una función:
Gs=A/(tau*s+1);
donde A=w(rad/seg)/V(voltios)
y tau=Ts/5; donde T es el tiempo de establecimiento que en tu gráfica parece ser 0.6seg. Debes estar usando un motor grandecito o conectado a carga para que te de este valor alto seguramente.
luego lo discretizad con c2d:
Gz=c2d(Gs,0.01,'zoh')
saludos.
Hola a todos.
Amigos tengo un proyecto para una materia y me gustaria que me ayudaran, pero la verdad es que no quiero nada hecho, lo que deseo es que me ayuden a hacerlo paso a paso, uds me guian y yo lo hago, quiza algunas cosas ya esten hechas, pero si no hago las cosas no aprendo y esa es mi mas grande ilusion.
Objetivo: Crear un control PID digital usando pic16f877a para controlar la temperatura de un horno.
Datos: La planta sera un cajon de madera y en su interior estara una bombilla de 60W la cual hara las veces de calefactor, esta bombilla se controlara usando pwm, el cual con la ayuda de un circuito adicional controlara la intensidad luminica de la bombilla y a su vez el calor generado. Para sensar la temperatura se usara un LM35 el cual da 10mV/ºC
Nota: Ire editando este post a medida como vaya realizando nuevas cosas y recibiendo sugerencias y correcciones de su parte.
PASO 1:
* Armamos el circuito sensor usando el LM35, probamos de que sense bien la temperatura y lo introducimos en nuestra caja de madera, previamente elaborada; hay que tener en cuenta de que la caja quede bien sellada una vez se cierre, o habra problemas al momento de controlar y sensar la temperatura.
PASO 2:
* Graficamos los datos y hallamos la ecuacion caracteristica de la planta, esta ecuacion es distinta para cada planta, podemos usar un cicruito/programa sensor de temperatura que nos ayude a tomar los datos de la temperatura cada cierto intervalo de tiempo.... o ..... usamos lapiz, papel y cronometro para hacer todo manual... yo lo hare manual.
-- Con estos datos pasamos a elaborar la grafica, manualmente, usando excel o usando matlab, como aqui lo que queremos es aprender, pues usaremos matlab, para eso ingresamos los datos en matlab en forma de dos matrices (tiempo/temperatura) y que nos la grafique.
Codigo para graficar los datos (esto es un ejemplo, los datos reales son demasiados):
clc;
tiempo[0 1 2 3 4 5 6 7 8 9 10];
temperatura[25 25 26 27 28 29 30 30 31 31 31];
plot(tiempo,temperatura)
Listo ahora tenemos graficada la respuesta de la planta ante una entrada escalon unitario.
Lo siguiente lo pueden hacer como lo explico a continuacion o como lo explica blackcat en un post mas adelante, los dos metodos son validos, ya cada quien escogera.
Para hallar la ecuacion, nos vamos a matlab (donde esta el workspace), alli veremos un boton que dice START (con el logo de matlab), pulsamos ese boton, luego click en toolboxes, luego en curve fitting y por ultimo en curve fitting tool: (figura 1)
(http://img19.imageshack.us/img19/1746/piddigitalf1yy2.th.jpg) (http://img19.imageshack.us/my.php?image=piddigitalf1yy2.jpg)
Figura 1
Se nos abrira una ventana nueva y alli daremos click en el boton Data... en la ventana que se abrio en X data ponemos la variable tiempo y en Y data la variable temperatura (o sus los nombres que uds les tengan asignados) .. si no les aparece nada es por que tienen que haber corrido su codigo de graficacion de sus datos reales para que las variables del tiempo y temperatura aparezcan en el workspace...
luego daremos click en Create date set. y click en close (figura 2)
(http://img209.imageshack.us/img209/7956/piddigitalf2jn7.th.jpg) (http://img209.imageshack.us/my.php?image=piddigitalf2jn7.jpg)
Figura 2
La grafica mostrada sera muy similar a la de su figure 1, pero en forma de puntos, lo cual no es ningun problema, por que ahora vamos a hallar nuestra ecuacion, para eso hacemos lo siguiente:
click en el boton fitting, luego en new fit .... aqui es donde les voy a recomendar algo a titulo personal...
Recomendacion: traten de la ecuacion quede en terminos de eulers y no en forma de polinomio, ya que si es asi todos los polos les quedaran en el mismo lugar (en cero) y no es muy buena idea para los calculos, pero en fin ya sabran uds como hacerlo...
Vamos a modificar unicamente donde dice Tye of fitt. ahi podran ver todo el tipo de ecuaciones con las cuales se puede representar nuestra curva caracteristica de la planta... si siguen mi recomendacion entonces buscar donde dice exponential y escogemos el tipo de ecuacion ( a*exp(bx) + c*exp(dx) ) y hacemos click en apply... luego de eso nos aparecera nuestra ecuacion (la que nos describe nuestros datos) ... si tienen dudas de que valores para las constantes escoger, escojan los promedios, osea los valores que estan fuera de los parentesis ... (figura 3)
(http://img187.imageshack.us/img187/1171/piddigitalf3ig7.th.jpg) (http://img187.imageshack.us/my.php?image=piddigitalf3ig7.jpg)
Figura 3
Ahora damos click en close y confirmamos que la ecuacion reproduce lo mas fielmente posible la grafica de nuestros datos si es asi ya tienemos nuestro y(t) = a*exp(bt) + c*exp(dt) y tambien los valores de las constantes... y si no s asi buscamos la que mejor lo haga.. pero estoy un 90% seguro de que la funcion exponential lo hara de maravilla.
Por lo tanto obtenemos una ecuacion y(t) = a*exp(b*t) + c*exp(d*t) aplicando laplace obtenemos:
32.98 s + 0.06097
Y(S) = ----------------------------
s^2 + 0.00141 s - 9.773e-009
Para mi caso ... a uds les dara una ecuacion similar .. matlab tambien lo hace usando el codigo
syms t
a = 42.98;
b = -6.9x10^-6;
c = -10.15;
d = 0.001417;
yt = a*exp(b*t) + c*exp(d*t)
YS = laplace (yt)
pero la verdad es que se ve muy raro y desordenado, asi que recomiendo lapiz y papel que hacerlo no es dificil y ademas sale directo, luego se aplica un poco de algebra, se reemplazan valores y listo.
Ahora sabemos que dentro de Y(S) tenemos implicito nuestra entrada escalon unitario, po lo tanto nuestra G(S) sera:
32.93 s^2 + 0.06097 s
G(S) = ----------------------------
s^2 + 0.00141 s - 9.773e-009
Ahora tenemos que determinar el tiempo de muestreo (T) para eso podemos usar el metodo que describe blackcat mas adelante o podemos asumirlo, ya cada quien decidira, por mi parte lo sumire como T=0.1
Pasamos a la parte en que tenemos que aplicar la transformada Z para poder discretizar nuestro G(S) y obtener nuestro G(Z) .. pero recuerden que G(Z) no es igual a la transformada Z de G(S) si no que lleva incorporado un retenedor de un orden cualquiera, para este caso usaremos uno de orden cero... y lo llamaremos Bo(S) asi que:
G(Z) = Z/ Bo(S) G(S) donde Z/ sera la tranformada Z de todo eso... bueno para los amantes del lapiz y papel, podeis hacerlo, por que por mi parte no lo voy a hacer, para eso voy a usar matlab y el siguiente codigo.
n1=32.98;
n2=0.06097;
d1=1;
d2=0.0014101;
d3=-(9.773/1000000000);
GS=tf( [n1 n2 0], [d1 d2 d3])
T=0.1;
GZ=C2D(GS,T,'zoh')
donde zoh indica a matlab que se va a usar un retenedor de orden cero. osea nuestro Bo(S)
y el resultado que nos entrega matlab sera el siguiente:
3.293 z - 3.292
G(Z) = ------------------
z^2 - 2 z + 0.9999
Sampling time: 0.1
Bueno hasta aqui todos deberiamos llevar algo similar, ahora es donde yo me he decidido por hacer un controlar por el metodo de ragazzini .. asi que por ahi es donde lo haremos....
lo primero es hallar los polos deseados teniendo en cuenta las siguientes condiciones de lazo cerrado.
essp = 0 ; Test 5% = 10 seg ; T = 0.1 ; Zcita = 0.5
usando formulas ya conocidas en nuestros estudios realizados obtenemos los dos polos deseados Z(1,2):
Z(1,2) = 0.9691 +/- j 0.0504 lo cual llamaremos (alfa +/- j beta) para comodidad
El metodo ragazzini nos dice que:
F(Z) = K / P(Z)
ahora P(Z) = (Z - Z1) (Z - Z2) como ya sabemos los valores de Z1 y Z2 entonces reemplazamos... pero podemos hacer algo mejor y es reemplazar por las letras y despejar para que nos quede algo asi:
P(Z) = Z^2 - 2 alfa Z + (alfa^2 + beta^2) ahora se reemplazan los valores y listo.
para que ragazzini funcione se dice que F(1) = 1 .. lo cual es de gran ayuda.. sabemos por formula que:
F(Z) = K / P(Z) --- reemplazamos P(Z) y nos quedara asi:
F(Z) = K / Z^2 - 1.9382 Z + 0.94169 --- hacemos F(1) = 1 --- quedando: 1 = K / 1^2 - 1.9382 + 0.94169 -- despejando K
K = 0.00349 ahora pasamos a la parte interesante...
D(Z) = [ 1/G(Z) ] [ F(Z) / 1 - F(Z) ]
pero como todo se puede simplificar entonces lo hacemos.
D(Z) = [ 1 / G(Z) ] [ K / P(Z) - K ]
como ya conocemos todas las variables y ecuaciones, simplemente hallamos D(Z)
D(Z) = [ 0.001058Z^2 - 0.002116Z + 0.001058 ] / [ Z^3 - 2.93769Z^2 - 0.9379 ]
Bueno la parte final es hallar la ley de control, para eso cambiamos D(Z) por u(Z)/ e(Z) y multiplicamos al otro lado por Z^-3 / Z^-3
quedandomos asi:
u(Z) 0.001058Z^-1 - 0.002116Z^-2 + 0.001058Z^-3
----- = ---------------------------------------------------------
e(z) 1 - 2.93769Z^-1 + 0.9382Z^-2 - 0.9379Z^-3
ahora el denominador pasa a multimplicar a u(Z) y e(Z) a multiplicar al numerador. luego de lo cual despejamos u(Z) y de paso usamos una un teorema que dice en resumen que u(Z)Z^-1 = u(Z-1) .. lo mismo con e(Z) .... ahora cambiamos Z por k quedando finalmente la ley de control asi:
u(k) = 0.001058 e(k-1) - 0.002116 e(k-2) + 0.001058 e(k-3) + 2.93769 u(k-1) - 0.9382 u(k-2) + 0.9379 u(k-3)
esta ley de control u(k) es la que programamos para poder controlar nuestro sistema
SE EDITARA A MEDIDA QUE SE AGREGUEN NUEVAS COSAS....
Hola ...
Para poder estimar un sistema con matlab se utiliza la herramienta IDENT .. esta es muy facil de usar, lo que se necesita es lo siguiente:
-> un vector de entradas (x) ... la señal de estimulo que aplicaste al sistema
-> un vector de salidas (y) ... la señal de respuesta que obtuviste del sistema al aplicar la señal de entrada
-> un vector de tiempo (t) ... los tiempos de cada muestra, deben ser constantes, es decir, 1, 2, 3, 4 segundos etc ... y no 1, 1.5, 2, 3, 5, 8 .. etc..
Con eso IDENT identificará el sistema que se encuentra adentro.
(http://img218.imageshack.us/img218/551/sist1zd9.jpg)
Lo que necesitas es una señal de prende o apaga el bombillo .. o bien una señal de volaje RMS que le aplicas al bombillo, luego debes medir en tiempo real el comportamiento de la temperatura ante la señal. Un escalon como lo estas haciendo, es posible, sin embargo solo tenes el comportamiento de calentar y no de enfriar.
Para poner a prueba el IDENT me inventé un sistema tipico y muy facil de entender .. un filtro pasobajos RC de primer orden, representado en Laplace como:\begin{equation}
G(s) = \frac{1}{1 + 0.001s}
\end{equation}
Es decir R = 1k y C = 1uF ... siendo tau = 1ms;
>> R = 1000
>> C = 1e-6
>> G = tf( 1, [R*C 1] )
Lo que voy hacer es inventar un vector de entradas, estimulare el sistema con esas entradas y utilizaré IDENT para identificar el sistema; lo ideal es que obtenga un resultado como en de la ecuacion anterior. En el caso practico, no conozco la ecuación del sistema, pero lo que se acostumbra es estimular el sistema con diferentes vectores de entrada y estimar varios modelos y comparar el resultado con los otros vectores.
Como mencioné, es necesario que tenga un vector de datos de entrada y ver la respuesta del sistema (salida) ante esa entrada, el mejor vector de entrada es aquel que tenga un buen contenido de frecuencia, ruido por ejemplo; sin embargo, dependiendo del sistema esto se complica, por ejemplo, en un motor; lo recomendable es estimar el modelo matemáticamente y ver que frecuencias son las aptas para crear una señal de estímulo; por ejemplo, el modelo de un motor es parecido a un filtro pasobajos, entonces utilizamos una señal que contenga componentes de frecuencia dentro del ancho de banda del motor.
Como en este caso, el ejemplo es simulado, utilizo la función RAND para generar un vector de entrada aleatorio de 100 muestras:
>> x = rand([100, 1])
(http://img220.imageshack.us/img220/5128/entradaxs3.jpg)
Esta es mi señal de entrada, ahora, necesito el vector de salida, estimularé el sistema con la función LSIM, entre los datos de esta funcion está el vector de tiempo. Aqui voy a utilizar tiempo discreto y elejiré que el tiempo de muestreo sea un décimo de la constante de tiempo del sistema (tau = RC); en este caso, ts = 0.1ms. Genero un vector de tiempo de 100 datos desde 0 hasta 99*0.1e-3.
>> n = (0:99)'
>> ts = 0.1e-3
>> t = n*ts
Ahora utilizo LSIM ....
>> y = lsim(G, x, t)
La grafica que obtengo en el tiempo es:
(http://img220.imageshack.us/img220/9757/salidaspr6.jpg)
donde la linea azul es la entrada y la linea verde es la salida, veran que esta filtrada.
Ya tengo listo mi vector de entrada, salida y tiempo; esto es lo que tendria en la practica, una señal de entrada, el comportamiento de la salida y el tiempo en que fue tomado cada dato. En el caso de la temperatura, el vector de entrada seria pulsos rectangulares, que me dicen cuando se prendio el bombillo y cuando permanecio apagado; la salida seria el aumento de temperatura y la disminucion de temperatura, el tiempo seria cada cuanto tome la muestra, como es temperatura puede ser de 1 o 3seg, constantes y no interrumpidos ... en este caso, tomar los datos a reloj en mano no es conveniente.
Yo conozco el sistema; sin embargo, la idea es estimar el sistema en tiempo discreto utilizando el vector de entrada, salida y tiempo; luego transformar el sistema a tiempo continuo y comparar.
Ejecuto IDENT
>>ident
Y aparece una pantalla como la siguiente:
(http://img220.imageshack.us/img220/2124/ident1xz8.jpg)
Lo que hacemos es en IMPORT DATA seleccionamos TIME DOMAIN DATA ... esto porque nuestros datos de entrada, salida y tiempo estan en el tiempo. Aparece un cuadro como asi:
(http://img220.imageshack.us/img220/1276/ident2wy0.jpg)
Y lo lleno con la informacion que tengo, entrada (x), salida (y), tiempo de inicio (0), tiempo entre muestras (ts) como lo hice en la figura. Damos IMPORT ..
(http://img220.imageshack.us/img220/1257/iden3xn4.jpg)
En la columna de cuadros izquierda me aparece el bloque de datos que cargue en IDENT; ahi puedo cargar la cantidad de datos que desee y puedo combinar datos como de estimulo y verificacion, esto es que con un conjunto de datos puedo estimar el sistema y con otro conjunto puedo verificar si la estimacion está correcta.
Ahora, sobre ESTIMATE damos LINER PARAMETRIC MODELS ... aqui escogemos el método de estimación que mas nos guste o que mejor de resultados, sobre el tipo de estimador se puede encontrar información en internet:
http://www.ie.itcr.ac.cr/einteriano/control2/Laboratorio/3.Models.pdf (http://www.ie.itcr.ac.cr/einteriano/control2/Laboratorio/3.Models.pdf)
Yo elejí ARX tanto de orden 1 como de orden 2 ... entre mas orden es mejor la estimación; sin embargo, hay un compromiso, pues estimar un compensador para un sistema que tiene un orden alto es mas dificil que para uno bajo, entonces escogemos el modelo de menor orden. Seleccionando la casilla MODEL OUTPUT aparecerá una gráfica que nos dice cual modelo es el que mejor se ajusta a los datos de salida. Vemos que con el modelo de 1º orden tenemos un ajuste de 100.
(http://img294.imageshack.us/img294/4262/ident5kn6.jpg)
Ahora arrastramos el cuadro que dice ARX111 de la columna derecha sobre el cuadro TO WORKSPACE ... eso hara que los datos del modelo lo tengamos en la linea de comando. Para ver el modelo ejecutamos:
>> H = tf(arx111)
Transfer function from input "u1" to output "y1":
0.09516
----------
z - 0.9048
Transfer function from input "v@y1" to output "y1":
1.665e-016 z
------------
z - 0.9048
Input groups:
Name Channels
Measured 1
Noise 2
Sampling time: 0.0001
Entonces nuestro modelo estimado en tiempo discreto es:\begin{equation}
G(z) = \frac{0.09516}{z - 0.9048}
\end{equation}
Usamos:
>> G2 = d2c(H)
y obtenemos:\begin{equation}
G(s) = \frac{1000}{s + 1000}
\end{equation}
Vemos que IDENT identificó el sistema correctamente ... Utilizando este modelo podemos usar SISOTOOL y calcular de manera sencilla un compensador :)
Saludos y espero que les sirva de algo!
PD: Espero sus comentarios y estamos listos para la implementación PID en un PIC!
tiempo=[0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210];
%son voltios de un divisor, ntc 10k al positivo, 3k3 a masa.
%medido en la de 3k3 a masa.
temperatura=[4.77 4.69 4.61 4.51 4.40 4.29 4.14 4.04 3.94 3.86 3.78 3.70 3.64 3.60 3.55 3.51 3.49 3.46 3.40 3.36 3.34 3.34];
plot(tiempo,temperatura)
%en eje x tiempo. Eje y voltios o temperatura.
syms t
a = 3.406;
b = -0.1962;
c = 0.3603;
d = 0.5446;
yt = a*exp(b*t) + c*exp(d*t)
YS = laplace (yt)
%yt =
%(1703*exp(-(981*t)/5000))/500 + (3603*exp((2723*t)/5000))/10000
%YS =
%1703/(500*(s + 981/5000)) + 3603/(10000*(s - 2723/5000))
g1=tf([1703 0],[500 500* 981/5000])
g2=tf([3603 0],[10000 -(10000*( 2723/5000))])
g=g1+g2
%rltool(g)
step(g)
K =0.94116;
Tp1 = 6.2583e-05;
Tp2 =1.9424;
num = K;
den = conv([Tp1 1],[Tp2 1]);
sys = tf(num,den)
%sys =
% 0.9412
% ---------------------------
% 0.0001216 s^2 + 1.942 s + 1
hola , no he leido el post completo, pero cuando tuve que identificar una planta de nivel de agua usando matlab tuve que hacer lo siguiente:
-tomar como T0(tiempo cero) el tanque vacio y medir el voltaje de salida de mi sensor
-empezar a llenar hasta cierto nivel, y tomar el tiempo que tardó en llegar a ese nivel y medir el voltaje de salida del sensor
-así mismo para las siguientes muestras
Matlab modela con parámetros en el dominio del tiempo, asi que usé el comando "Ident" y le pase los parametros en forma de vector del tiempo y el voltaje
¿No puedes ajustar el PID a mano? Creo que será más sencillo.
Un problema que te vas a encontrar en el ventilador, es que su comportamiento no es lineal. ¿Lo estás controlando con tensión o con PWM?
Saludos.
ojo, que con los sistemas por aire, en ocasiones hay que aplicar soluciones anti WIND-UP.
un ejemplo.
Click Click (http://www.elai.upm.es/webantigua/spain/Asignaturas/ControlProcesos/archivos/Practicas/Practica_1.pdf)
Yo probé a hacer un control PID para controlar la temperatura de un disipador de pc y acabé volviendome loco
Yo tengo como salida el array de temperatura, y como tiempo desde 0 en pasos de 10seg. ¿que pongo como entrada? sin array de entrada no se puede ejecutar "process models".
La idea es aprender a utilizar ident para encontrar la función de transferencia de una planta y rltool para calcular el controlador.
Yo tengo como salida el array de temperatura, y como tiempo desde 0 en pasos de 10seg. ¿que pongo como entrada? sin array de entrada no se puede ejecutar "process models".
La idea es aprender a utilizar ident para encontrar la función de transferencia de una planta y rltool para calcular el controlador.
La planta será muy parecida a un sistema con un polo.
Puedes calentar las resistencias hasta que la temperatura sea estable sin tensión en el ventilador (Vin = 0v)
Luego pones en marcha el ventilador con 6 voltios y tomas medidas cada segundo:
Vin = [6 6 6 6 6 6 6 6 6 6 6 6 6 ....
Vout = [4 3.8 3.5 3.2 2.9 2.6 2.3 2.1 1.9 1.8 1.7 1.6 1.6 ....
El array de entrada es la tensión del ventilador.
Saludos.
Yo tengo como salida el array de temperatura, y como tiempo desde 0 en pasos de 10seg. ¿que pongo como entrada? sin array de entrada no se puede ejecutar "process models".
La idea es aprender a utilizar ident para encontrar la función de transferencia de una planta y rltool para calcular el controlador.
La planta será muy parecida a un sistema con un polo.
Puedes calentar las resistencias hasta que la temperatura sea estable sin tensión en el ventilador (Vin = 0v)
Luego pones en marcha el ventilador con 6 voltios y tomas medidas cada segundo:
Vin = [6 6 6 6 6 6 6 6 6 6 6 6 6 ....
Vout = [4 3.8 3.5 3.2 2.9 2.6 2.3 2.1 1.9 1.8 1.7 1.6 1.6 ....
El array de entrada es la tensión del ventilador.
Saludos.
Muchas gracias a todos!!
Con ese dato he conseguido llegar a buen puerto, o eso creo, el invento mantiene la temperatura al nivel indicado.
Dejo todo esto documentado por si le sirve a alguien ;-)
Saludos!!
http://www.4shared.com/rar/sAJauyd8ce/PID_Inverso_enfriador.html