CONTROL DE PROCESOS QUÍMICOS - Prof. Juan … · Funciones elementales Matlab dispone de las...
Transcript of CONTROL DE PROCESOS QUÍMICOS - Prof. Juan … · Funciones elementales Matlab dispone de las...
UNIVERSIDAD NACIONAL EXPERIMENTAL POLITECNICA
“ANTONIO JOSÉ DE SUCRE”
VICERRECTORADO BARQUISIMETO
DEPARTAMENTO DE INGENIERÍA QUÍMICA
CONTROL DE PROCESOS QUÍMICOS
Prof: Ing. (MSc).
Juan Enrique Rodríguez C.
Octubre, 2013
Matlab es un programa interprete de comandos. Esto quiere
decir que es capaz de procesar de modo secuencial una serie de
comandos previamente definidos, obteniendo de forma
inmediata los resultados. Los comandos pueden estar ya
definidos en el propio Matlab y pueden también ser definidos
por el usuario.
Operadores
Hay operadores para números (reales o complejos) y para matrices.
• Para números: + - * / ^
• Números complejos: Esta definida la unidad imaginaria, √-1, que se denota indistintamente por los
símbolos i y j
• Para matrices: + - * / ^
• Para matrices elemento por elemento: .+ .- .* ./ .^
Funciones elementales
Matlab dispone de las funciones elementales mas comunes (las que tienen las calculadoras de bolsillo) y
otras especiales, propias. Realizan una operación sobre un argumento numerico dado de tipo matriz y
operan elemento por elemento. Las mas usuales son:
• Trigonometricas: sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh.
• Lógicas: any, all, find, exist, isnan, finite, isempty, isstr, strcomp.
• Otras: abs, angle, sqrt, real, imag, conj, round, x, floor, ceil, sign, rem, exp, log, log10.
• Especiales: bessel, gamma, rat, ert, invertf, ellipk, ellipj. 3
Operaciones numéricas
>> a = 4
a =
4
>> b = 5 + a
b =
9
>> c = a^2 + b^2
c =
97
>> sin (30*pi/180)
ans =
1/2
Si el usuario no ha asignado el
resultado a una variable, Matlab lo
hace utilizando la variable \ans"
Operaciones con vectores
>> a = [1 2 3 4 6 4 3 4 5]
a =
1 2 3 4 6 4 3 4 5
>> b = a + 2
b =
3 4 5 6 8 6 5 6 7
>> c = a + b
c =
4 6 8 10 14 10 8 10 12
>> d = a .* b
c =
3 8 15 24 48 24 15 24 35
5
Operaciones con polinomios
6
Los polinomios se representan en Matlab como vectores fila. Por ejemplo, el polinomio
3s3 - 5s2 + 7s + 3 se representa por
>> p = [ 3 -5 7 3 ]
Las raíces de un polinomio se hallan mediante la función roots:
>> r = roots(p)
La función polyval sirve para hallar el valor de un polinomio. Si el parámetro que le pasamos
es un vector, calcula otro vector con los valores del polinomio para cada uno de los del vector.
Veamos como se efectúan algunas de las operaciones mas comunes con matrices:
Introducir una matriz A:
Operaciones con matrices
>> A = [1 2 3 ; 4 5 6 ; 7 8 0]
• Calculo de la transpuesta:
>> B = A'
• Producto matricial:
>> C = A * B
• Determinante:
>> det(A)
• Matriz inversa:
>> inv(A)
7
• Valores propios y vectores propios:
>> [val,vec]=eig(A)
• Valores singulares:
>> svd(A)
• Exponenciación matricial (eA):
>> expm(A)
• Polinomio característico:
>> p = poly(A)
La librera Symbolic Math Toolbox da acceso a Matlab a algunas funciones del núcleo de Maple
que permiten operar con expresiones simbólicas. las variables simbólicas se hayan declarado
previamente con sym() o syms. Por ejemplo:
>> syms a b p x
>> a = x^3 + 3*x^2 - 2*x + 7;
>> b = x^2 + x + 3;
>> p = a * b
p =
(x^3 + 3*x^2 - 2*x + 7)*(x^2 + x + 3)
>> expand(p)
ans =
x^5 + 4*x^4 + 4*x^3 + 14*x^2 + x + 21
Operaciones simbólicas
8
La sustitución de un símbolo por otro en una expresión simbólica se puede realizar con la orden
subs. La forma de hacerlo es subs(expr, old, new)
Sustitución de variables
>> syms x a b c
>> f = a*x^2 + b*x + c;
>> g = subs(f,x,-1)
g =
a - b + c
Las ordenes poly2sym y sym2poly sirven, respectivamente, para convertir un polinomio
expresado en forma numérica (vector de coeficientes) en su expresión simbólica, y viceversa.
El siguiente ejemplo ilustrara su utilización.
>> syms x
>> p = [1 2 3 4 5]
>> px = poly2sym(p,x)
px = x^4 + 2*x^3 + 3*x^2 + 4*x + 5
>> sym2poly(px)
ans =
1 2 3 4 5
Conversión de polinomios
Se pueden también sustituir varias variables a la vez.
>> syms x a b c k
>> f = a*x^2 + b*x + c;
>> g = subs(f,[a,b,c],[1,2,k])
g =
x^2+2*x+k
9
El comando diff() de Matlab permite calcular derivadas, totales y parciales, de una expresión
algebraica, función de una o varias variables y parámetros, respecto de una de ellas (o de ellos).
Supongamos que nos dan una expresión f(x), por ejemplo el polinomio
f(x) = a3x3 + a2x
2 + a1x + a0
Derivadas
>> syms a3 a2 a1 a0 x
>> f=a3*x^3+a2*x^2+a1*x+a0;
>> diff(f)
ans =
3*a3*x^2 + 2*a2*x + a1
El comando int() de Matlab permite resolver integrales, tanto indefinidas como definidas.
Integrales
>> syms a b x
>> f=log(x)/x;
>> int(f)
ans =
log(x)^2/2
Para obtener la expresión de una integral definida, pondremos
>> int(f,a,b)
ans =
1/2*log(b)^2-1/2*log(a)^2
10
En las aplicaciones interesa a veces conocer el valor numérico de una función y = f(x) para uno
o varios valores de la variable. Por ejemplo, dada la función y = 10(1-e-x/3sin(10x)), definida en
el intervalo [0; 15], una posible representación en Matlab, seguida de su representación grafica,
será:
>> x=[0:0.1:15];
>> y=10*(1-exp(-x/3).*sin(10*x));
>> plot(x,y),title('Grafica de una función')
Gráficos en 2D
11
Las funciones de dos variables, de la forma f(x;y) se pueden representar gráficamente con
Matlab en 3D. Veámoslo con un ejemplo.
Gráficos en 3D
yyyxxxxz 22623 2224
Supongamos que queremos calcular los valores de z en una región rectangular del plano (x; y)
definida en el rango -3 < x < 3 y -3 < y < 13 y representarla gráficamente. Para ello
escribiremos:
>> x1 = linspace(-3, 3);
>> y1 = linspace(-3, 13);
>> [x,y] = meshgrid(x1, y1);
>> z = x.^4+3*x.^2-2*x+6-2*y.*x.^2+y.^2-2*y;
>> mesh(z)
>> surf(z)
12
Descomposición en fracciones parciales
Considere la función de transferencia
611s6ss
63s5s2s
sA
sBsG
23
23
Para esta función:
num=[2 5 3 6]
den=[1 6 11 6]
Escribimos:
[r,p,k]=residue(num,den)
Por lo que queda expresado:
21s
3
2s
4
3s
6
sA
sBsG
13
Matlab esta dotado de un mecanismo que le permite interpretar ficheros de texto, con la
condición de que su nombre termine por .m. Se utilizan principalmente para crear funciones (en
el sentido matemático), programas y funciones (ordenes) de Matlab.
Ficheros-m
14
t=0:0.05:10;
y=sin(t);
z=cos(t);
plot(t,y,'o',t,z,'x')
grid
title('Graficas del seno y coseno')
xlabel('seg')
ylabel('y=seno(t); z=coseno(t)')
text(3,0.45,'sen(t)')
text(0.8,-0.3,'cos(t)')
15
Se desarrolla una reacción termoquímica en donde el reaccionante A se convierte en un producto B.
Velocidad de reacción: r(t)=k*CA(t)
Concentración de la entrada: CAi(t)
En la condición inicial
t=0; CAi(0)=12,5 mol/L
Demás constantes
Constante de velocidad de reacción: k=0,2 min-1
Volumen de la masa reaccionante: V=5 Litros
Flujo de entrada: Fi= 1 L/min
Ecuaciones diferenciales de primer orden
Balance para el componente A (E-C+F-S=A) en estado no estacionario:
tCFV*k
FV*ktC*
FV*k
F
dt
tdC*
FV*k
V
tC*FV*kF*tCdt
tdC*V
dt
tdC*VF*tCV*tC*kF*tC
dt
tdC*V
dt
V*Cd
dt
ndF*tCV*rF*tC
AAiiA
AiAiA
AAAiAi
AAAAiAi
16
tC*KtC
dt
tdC*τ
tC*FV*k
FtC
dt
tdC*
FV*k
V
AieAA
Aii
AA
5,0
L/min 1L 5*min 0,2
L/min 1
FV*k
FK
min5,2L/min 1L 5*min 0,2
L 5 τentonces ,
FV*k
V τSea
1-e
1-
mol/L 6,25(0)C
valoresdoreemplazan
00C*F0C*V*k0C*F
iónconcentrac la de inicialCondición
A
AAAi
.m matlaben programaun en problema este Realicemos
20
Respuestas transitorias: se refiriere a la que va del estado inicial al estado final, tales como respuesta a
un salto o entrada escalón, se utilizan frecuentemente para investigar las características en el dominio
temporal de los sistemas de control.
Respuesta a una entrada escalón
Sea la siguiente función de transferencia 254ss
25
sR
sC2
23
Sea la respuesta a un impulso unitario del sistema de segundo orden
12ζs
s
sR
sCsG
2
s
Para obtener la respuesta a un impulso unitario, realicemos el siguiente procedimiento
num=[1 0]
den=[1 2*(zeta) 1]
Consideremos cinco valores diferentes de zeta: ζ = 1, 0.7, 0.5, 0.3, y 0.1
26
Sea la siguiente función de transferencia
2010ss
1
sR
sCsG
2
El objetivo de este problema es mostrar como Kp, Ki y Kd contribuyen a obtener
• Un tiempo de subida más rápido
• Un sobreimpulso mínimo
• Un error en régimen permanente nulo
Veamos, en primer lugar, la respuesta en bucle abierto del sistema ante una entrada escalón. Cree
un nuevo archivo de instrucciones con el siguiente código:
num=[1];
den=[1 10 20];
step(num,den)
27
Ahora diseños un controlador que reduzca el tiempo de subida y el tiempo de establecimiento y
elimine el error en régimen permanente.
Control Proporcional (P)
De acuerdo a lo estudiado en las clases anteriores, el controlador proporcional (Kp) reduce el
tiempo de subida, incrementa el sobreimpulso y reduce el error en régimen permanente. Por lo
que la función de transferencia de lazo cerrado del sistema anterior con un controlador
proporcional sería:
Kp20s*10s
Kp
sX
sY2
Establezcamos 300 como valor de la ganancia proporcional (Kp) y modifiquemos el archivo de
instrucciones de la siguiente manera:
Kp=300;
num=[Kp];
den=[1 10 20+Kp];
t=0:0.01:2;
step(num,den,t)
Al ejecutar este archivo en la ventana de instrucciones de matlab se deberá obtener la siguiente
figura:
28
La función de Matlab denominada cloop puede usarse para obtener directamente la función de
transferencia de lazo cerrado a partir de la función de transferencia de lazo abierto (en vez de
obtener la función de lazo cerrado a mano). El siguiente archivo de instrucciones en el que se usa
el comando cloop generará una gráfica idéntica a la mostrada anteriormente.
num=[1];
den=[1 10 20];
Kp=300;
[numCL,denCL]=cloop(Kp*num,den);
t=0:0.01:2;
step(numCL,denCL,t)
29
El gráfico anterior muestra como el controlador proporcional ha reducido tanto el tiempo de
subida como el error en régimen permanente, ha incrementado el sobreimpulso y disminuido, en
una pequeña cantidad el tiempo de establecimiento.
Control Proporcional y Derivativo (PD)
Veamos ahora el control PD. Este tipo de controlador reduce el sobreimpulso como el tiempo de
establecimiento. La función de transferencia de lazo cerrado del sistema propuesto con un
controlador PD es:
Kp20s*K10s
Kps*K
sX
sY
D
2
D
Establezcamos 300 como valor de la ganancia proporcional (Kp) y en 10 el de KD. Introduzca las
siguientes instrucciones y ejecútelo desde la ventana de Matlab:
Kp=300;
KD=10;
num=[KD Kp];
den=[1 10+KD 20+Kp];
t=0:0.01:2;
step(num,den,t)
Al ejecutar este archivo en la ventana de instrucciones de matlab se deberá obtener la siguiente
figura:
30
Esta gráfica muestra como el controlador derivativo reduce el sobreimpulso y el tiempo de
establecimiento y tiene poco efecto sobre el tiempo de subida y el error en régimen permanente
Control Proporcional e Integral (PI)
El control PI, disminuye el tiempo de subida, incrementa tanto el sobreimpulso como el tiempo
de establecimiento y elimina el error en régimen permanente. Para el sistema propuesto, la
función en lazo cerrado para el control PI es:
I
23
I
Ks*Kp20s*10s
Ks*Kp
sX
sY
Reduzcamos a 30 el valor de Kp y establezcamos en 70 el de KI. Introduzca las siguientes
instrucciones y ejecútelo desde la ventana de Matlab:
31
Kp=30;
KI=70;
num=[Kp KI];
den=[1 10 20+Kp KI];
t=0:0.01:2;
step(num,den,t)
Al ejecutar este archivo en la ventana de instrucciones de matlab se deberá obtener la siguiente
figura:
Se ha reducido el valor de la ganancia proporcional (Kp) porque el controlador integral también
reduce el tiempo de subida e incrementa el sobreimpulso tal y como hace el controlador
proporcional. La respuesta anterior muestra como el controlador integral elimina el error en
régimen permanente.
32
Control Proporcional, Integral y Derivativo (PID)
Veamos ahora el controlador PID, la función de transferencia de lazo cerrado del sistema dado
con un controlador PID es:
I
2
D
3
IP
2
D
Ks*Kp20s*K10s
Ks*Ks*K
sX
sY
Tras varias pruebas se ha comprobado que los valores de las ganancias KP=350; KI=300 y
KD=50, proporcionan una respuesta adecuada. Para confirmarlo, introduzca las siguientes
instrucciones en matlab.
KP=350;
KI=300;
KD=50;
num=[KD KP KI];
den=[1 10+KD 20+KP KI];
t=0:0.01:2;
step(num,den,t)
Al ejecutar este archivo en la ventana de instrucciones de matlab se deberá obtener la siguiente
figura:
33
Consejos generales para el diseño de un controlador PID
1.- Obtenga la respuesta en lazo abierto y determine los parámetros que deben ser mejorados
2.- Añada un control proporcional para mejorar el tiempo de subida
3.- Añada un control derivativo para mejorar el sobreimpulso
4.- Añada un control integral para eliminar el error en régimen permanente
5.- Ajuste los valores de KP, KI y KD para obtener la respuesta deseada.
6.- Tenga en cuenta que, si no es necesario, no tiene porque implantar los tres controladores en un
único sistema
35
Como el lugar de las raíces es realmente las ubicaciones de todos los posibles polos de lazo
cerrado, del lugar de las raíces se puede seleccionar una ganancia(para un controlador
proporcional puro) tal que el sistema de lazo cerrado se comporte de manera adecuada.
Trazando el lugar de las raíces de una función de transferencia
Considere un sistema en lazo abierto, cuya función de transferencia sea:
20s*15s*5s*s
7s
sX
sY
Primero, debemos desarrollar el polinomio que se encuentra en el denominador, quedando:
num=[1 7];
den=[1 40 475 1500 0];
rlocus(num,den)
Para lograr determinar el lugar de las raíces, introduzca las siguientes instrucciones en matlab.
s*1500s*475s*40s
7s
sX
sY234
36
¿Cómo se puede diseñar un controlador para este sistema usando el método del lugar de las raíces?
Ejemplo: Digamos que nuestro criterios de diseño para el controlador sea 5% de sobreimpulso y
un (1) segundo de tiempo de subida.
Usar la función del Matlab denominada sgrid sirve para hallar una región aceptable del lugar
de raíces. Esta función sgrid requiere dos argumentos: Frecuencia natural (Wn) y coeficiente
de amortiguamiento (zeta). Estos dos argumentos pueden determinarse a partir de los
requerimientos de tiempo de elevación, tiempo de establecimiento, sobrepico y las 3
ecuaciones de abajo.
2
p
2
p
r
n
s
n
π
MLn1
π
MLn
ζ
T
1,8ω
T
4,6ω*ζ
Donde:
Wn=Frecuencia natural
ζ=coeficiente de amortiguamiento
Ts=Tiempo de establecimiento
Tr=Tiempo de Subida
Mp=Máximo sobrepico
37
05,0esoSobreimpul2ζ1
*ζπ
Recordemos que existe una relación entre el sobreimpulso y la razón de amortiguamiento, la cual es:
Operando y resolviendo en la ecuación, se tiene que la razón de amortiguamiento (ζ, zeta) vale
o tiene que ser mayor a: 0,6901
Estudiando el sistema a lazo cerrado, tenemos:
num=[1 7];
den=[1 40 475 1500 0];
step(num,den)
[numCL,denCL]=cloop(num,den)
Obviamente, no todos los polos satisfacen los
criterios de diseño. Para determinar que parte del
lugar es aceptable usamos el comando:
sgrid(zeta,Wn) que traza las líneas de razón de
amortiguamiento (zeta) y frecuencia natural (Wn)
Además, de la ecuación que relaciona la frecuencia natural y la razón de amortiguamiento, se tiene:
8,11
8,1
T
1,8ω
r
n
38
Ahora, introduzca las siguientes instrucciones en matlab.
rlocus(num,den)
zeta=0.7;
Wn=1.8;
sgrid(zeta,Wn)
En la grafica anterior, las
dos líneas discontinuas
de un ángulo de
aproximadamente 45º
indican la ubicación de
los polos con zeta=0.7; el
semicírculo indica la
ubicación de los polos
con una frecuencia
natural Wn=1.8
39
Para hacer el sobreimpulso < al 5%, los polos deben estar entre las dos líneas discontinuas
Para hacer el tiempo de subida sea < al 1 seg, los polos deben estar fuera del semicírculo.
Ahora conocemos que sólo es aceptable la parte del lugar de las raíces comprendida entre las dos
líneas y fuera del semicírculo. Como todos los polos comprendidos en esta zona están en el
semiplano izquierdo, el sistema será estable en lazo cerrado.
Para ello, se puede usar el comando de matlab: rlocfind para elegir la ubicación de los polos
deseados, como se muestra a continuación:
[kd,poles]=rlocfind(num,den)
Select a point in the graphics window
selected_point =
-3.2611 - 0.0074i
kd =
298.0271
poles =
-21.7782
-12.6360
-3.2610
-2.3247
40
Continuemos y usemos la ganancia (kd) seleccionada en nuestro controlador proporcional.
Respuesta de lazo cerrado
Para obtener la respuesta ante una entrada en escalón se necesita conocer la función de
transferencia de lazo cerrado. Se puede hallar utilizando las reglas del algebra de bloques o dejar
que Matlab lo haga por nosotros, como se muestra:
[numCL,denCL]=cloop((kd)*num,den)
numCL =
1.0e+03 *
0 0 0 0.2980 2.0862
denCL =
1.0e+03 *
0.0010 0.0400 0.4750 1.7980 2.0862
Los dos argumentos de la función cloop son el numerado y el denominador del sistema en lazo
abierto. Se debe incluir la ganancia proporcional que previamente se ha elegido, obteniendo:
Comprobemos la respuesta del sistema en lazo cerrado ante un escalón unitario.
step(numCL,denCL)