Hola a todos:
Soy nuevo en este foro y estoy realizando un proyecto para una materia de la escuela, mm según yo estoy algo adelantado, ya he encontrado mi uk que es la ecuación que según yo debo programar en mi pic. La pregunta es que no encuentro como interpretar uk no se si es tiempo en el que debe estar encendido mi sistema o que?.
Mi sistema consta de un foco ubicado en una caja, la temperatura es censada por un sensor de temperatura. En este momento el sistema es un control on-off, es decir, el foco empieza a calentar cuando hace falta temperatura y se desactiva cuando la temperatura es mayor. MM como verán el sistema es fluctuante, pero eso no me sirve así que necesito implementar un sistema pid para que sea totalmente estable.
A continuación les muestro lo que llevo:
1.- Como primer paso necesito determinar la función de transferencia de mi sistema para eso he tomado las siguientes mediciones:
Tiempo en segundos: 0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210
Temperatura: 15,16,18,20,22,24,26,28,29,31,33,34,36,37,38,40,41,42,43,44,44,45
Con ayuda del siguiente código encontramos f(t) de mi sistema:
clc;
tiempo=[0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210];
temperatura=[15 16 18 20 22 24 26 28 29 31 33 34 36 37 38 40 41 42 43 44 44 45];
plot(tiempo,temperatura)
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.
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 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, despues click en el boton fitting, luego en new fit, 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, seleccionamos 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 paréntesis.
Bueno la ecuación que encontré fue:
f(x)=a*exp(b*x)+c*exp(d*x)
a=112.2
b=-0.0007449
c=-98.25
d=-0.003188
despues:
y(t)=a*exp(-b*t)-c*exp(-d*t)
con los valores:
a=112.2
b=0.0007449
c=98.25
d=0.003188
Aplicando Laplace:
*de tablas:
e^-at ---------------------> 1/s+a
Sustituyendo:
y(s)=a(1/s+b)-c(1/s+d)
y(s)=a/s+b-c/s+d
y(s)=as+ad-cs+cb/s^2+ds+bs+bd
Por lo tanto:
y(s)=13.95s+.430880025/s^2+3.9329x10^-3s+2.3747412x10^-6
Utilizando matlab con el siguiente codigo obtenemos:
Gs=tf([13.95 0.430880025],[1 0.0039329 0.0000023747412])
Gz=C2D[Gs,0.1,'zoh']
13.95s + 0.4309
G(s)= -------------------------------------
s^2 + 0.003933s + 2.375x10^-6
1.397z-1.393
G(z)= -------------------------------------
z^2 - 2z + 0.9996
Hasta este punto he obtenido mi funcion de transferencia.
2.- Como segundo paso me decidi por usar el método de ragazzini para encontrar uk:
K 1 F(z)
F(z) = --------- ; P(z) = (z - z1) (z - z2) ; D(z) = -------- -------
P(z) G(z) 1 - F(z)
Condicion de estabilidad -> F(1) = 1 ; z1 y z2 son polos deseados ; essp=0
Condiciones iniciales
Zita = 0.7 ; tiempo de establecimiento al 5% (Test5%) = 2seg ; essp = 0 ; T = 0.1
1) hallamos polos deseados (z1 y z2)
Wn = 3 / (Zita * Test5%)
|z| = euler^(-Zita * Wn * T)
<z = (T * Wn) * Raiz cuadrada(1 - (Zita^2)) --> resultado en radianes, se debe pasar a grados
Resultado en radianes * 180 / PI ---> nos entrega el resultado en grados
con una calculadora escribimos REC( |z| , <z ) ----> <z debe ser en grados
lo anterior nos entrega los polos deseados z1 y z2
z1,2 = Alfa +/- j Beta
z1,2 = 0.85084 +/- j 0.13105
2) hallamos K
P(z) = (z - z1) (z - z2) ----> simplificando -----> P(z) = z^2 - (2 * Alfa * z) + Alfa^2 + Beta^2
P(z) = z^2 - 1.70168 z + 0.74110
F(z) = K / P(z)
F(1) = 1 ----------> Se reemplaza P(z) y luego todas las z se cambian por un 1
K
F(z) = -----------------------------------
z^2 - 1.70168 z + 0.74110
K = (1)^2 - 1.70168 (1) + .74110 -----------> K = .03942
B(z) A(z) K
G(z) = --------- ; D(z) = -------- ---------------
A(z) B(z) [ P(z) - K ]
mi G(z):
1.397(z-.997)
g(z)= -------------------
(z-.99979998)^2
.02821z^2 -.05643z + .02820
Reemplazando: D(z) = -------------------------------------------
z^3 - 2.4153z^2 + 1.9161z - 0.50076
En este caso el z de mayor potencia es z^3
por lo tanto hay que multiplicar arriba y abajo por z^ -3
y reemplazar D(z) por u(z) / e(z)
u(z) .02821z^-1 - .05643^- 2 + .02820z^-3
------ = -----------------------------------------------------
e(z) 1 - 2.4153z^-1 + 1.9161z^-2 - 0.50076z^-3
Resumiendo, haber..... hay que despejar u(z)
cuando pase eso van a quedar partes como u(z) 0.79 z^ -3
por alla alguna vez nos enseñaron una propiedad que nos salva la patria
entonces seria algo como 0.79(z - 3) y asi con todo ya sea u(z) o e(z)
cambiamos las z por k (minusculas)
uk = ( .02821 ek-1 ) - ( .05643 ek-2 ) + ( .02820 ek-3 ) + ( 2.4153 uk-1 ) - ( 1.9161 uk-2 ) + ( .50076 uk-3)
uk = es la salida en el tiempo presente
uk-1 = es la salida en el tiempo anterior
ek-1 = error en el tiempo anterior
ek-2 = error en el tiempo anterior del anterior
Pues como se abran dado cuenta según yo tengo todo, el problema ahora es como manejar este uk para la programación en el micro
