Sistemas Automáticos (SAU) (Parcial 1) 21 Octubre 2015

11407 - Grado Ing. Tec. Industriales Tiempo: 1.5 horas

1. [Modelado y linealización] (Total: 5 pts.)

En el sistema de dos tanques de la figura, el depósito izquierdo tiene una toma de lı́quido elevada:

qi

qs K1
h2
K2
hc h1

A1 u A2

Los parámetros constantes son:

hc =0.5 m, A1 =0.6 m2 , A2 =1.2 m2 , K1 =0.3 m3 /(s · m1/2 ), K2 =0.8 m3 /(s · m1/2 )

K1 y K2 incorporan todos los valores constantes de los principios de conservación de la masa

(fluidos no compresibles) y energı́a de Bernouilli. Se pide:
(a) (1.5 pt.) Obtén el modelo que describe la evolución de los niveles de los depósitos (considerar
h2 > h1 ) en función del caudal de entrada qi y de la apertura de la válvula u, cuyo rango es
[0, 1].
Sol.: El modelo que describe la evolución del nivel h2 es:
dh2 p
A2 · = qi − K2 · u h2 − h1
dt
El modelo que describe la evolución del nivel h1 es:
dh1 p
A1 · = K2 · u · h2 − h1 − qs
dt
Dondel valor qs depende de la relación entre h1 y hc . Ası́
( √
K1 · h1 − hc , h1 ≥ hc
qs =
0, h1 < hc
Con lo que estas tres ecuaciones constituyen el modelo pedido (entradas u, qi ; incógnitas
h1 , h2 , qs ).

(b) (0.5 pt.) En el punto de funcionamiento se establece la posición de la válvula ueq =0.5, ¿Qué
valor del caudal de entrada qi,eq será necesario para dar lugar al nivel h1,eq =0.7? ¿Cuál será
la h2,eq ?
Sol.:
Dado que h1,eq = 0.7 > hc = 0.5, el modelo con el cual se obtiene el punto de equilibrio
solicitado es:
K2 p K1 p
ḣ1 = · u h2 − h1 − · h1 − hc
A1 A1
qi K2 p
ḣ2 = − · u · h2 − h1
A2 A2
Sistemas Automáticos (SAU), GITI-ETSII 21 Octubre 2015, pág. 2

En el equilibrio, las derivadas son cero:

p p
0 = K2 · ueq h2,eq − h1,eq − K1 · h1,eq − hc
p
0 = qi,eq − K2 · ueq · h2,eq − h1,eq
Sustityendo y despejando:
p √
0 = 0.8 · 0.5 h2,eq − 0.7 − 0.3 · 0.7 − 0.5 → h2,eq = 0.8125
√
0 = qi,eq − 0.8 · 0.5 · 0.8125 − 0.7 → qi,eq = 0.1342
Por tanto, el punto de equilibrio es:

(ḣ1,eq , h1,eq , ḣ2,eq , h2,eq , ueq , qi,eq ) = (0, 0.7, 0, 0.8125, 0.5, 0.1342)

(c) (2 pt.) Obtener un modelo linealizado alrededor del punto de equilibrio del apartado anterior.
Escribir el modelo linealizado en forma matricial normalizada.
Sol.:
Llamando

K2 p K1 p
g1 (ḣ1 , h1 , h2 , u) = ḣ1 − · u h2 − h1 + · h1 − hc = 0
A1 A1
qi K2 p
g2 (h1 , ḣ2 , h2 , u, qi ) = ḣ2 − + · u · h2 − h1 = 0
A2 A2
Y aplicando

dg1 (·) dg1 (·) dg1 (·) dg1 (·)

∆ḣ1 + ∆h1 + ∆h2 + ∆u = 0
dḣ1 eq dh1 eq dh2 eq du eq

Se llega a
!
K2 ueq K1
∆ḣ1 + p + p ∆h1
2A1 h2,eq − h1,eq 2A1 h1,eq − hc
Ku K2 p
− p 2 eq ∆h2 − h2,eq − h1,eq ∆u = 0
2A1 h2,eq − h1,eq A1
Sustituyendo (eliminando, por brevedad, la notación ∆ y suponiendo todas las variables en
coordenadas incrementales), se tiene:
ḣ1 = −1.5528h1 + 0.9938h2 + 0.4472u
Análogamente con g2 :
ḣ2 = 0.4969h1 − 0.4969h2 − 0.2236u + 0.8333qi
Con lo que estas dos ecuaciones constituirán el modelo linealizado pedido. En forma matri-
cial normalizada:
       
ḣ1 −1.5528 0.9938 h1 0.4472 0 u
= +
ḣ2 0.4969 −0.4969 h2 −0.2236 0.8333 qi
Sistemas Automáticos (SAU), GITI-ETSII 21 Octubre 2015, pág. 3

(d) (1 pt.) Con el modelo linealizado anterior, cuáles son los nuevos valores de equilibrio de
los niveles de los depósitos si se mantiene el mismo caudal de entrada y la válvula se cierra
ligeramente, posicionándola en ueq2 =0.45.
Sol.:
Si la válvula se cierra pasando del punto de linealización ueq =0.5 a un nuevo ueq2 =0.45,
entones ∆u = 0.45 − 0.5 = −0.05. Dado que el caudal de entrada se mantiene igual, se
tiene ∆qi = 0. Partiendo del modelo linealizado en forma matricial, en el nuevo equilibrio,
se tiene:
       
0 −1.5528 0.9938 ∆h1,eq 0.4472 0 −0.05
= + (1)
0 0.4969 −0.4969 ∆h2,eq −0.2236 0.8333 0
Despejando los niveles, se tiene:
   −1  
∆h1,eq 0.4472 0 −0.0224
=− (2)
∆h2,eq −0.2236 0.8333 0.0112
      
∆h1,eq −1.7889 −3.5778 −0.0224 0
=− = (3)
∆h2,eq −1.7889 −5.5903 0.0112 0.0225
Por lo que, el nivel h1 en el nuevo equilibrio sigue igual, mientral que h2 en el nuevo equi-
librio es h2,eq2 =0.8125+0.0225=0.835.
Sistemas Automáticos (SAU), GITI-ETSII 21 Octubre 2015, pág. 4

2. [Respuesta temporal] (Total: 5 pts.)

El circuito eléctrico de la figura, alimentado por una tensión de entrada V (t), tiene un modelo
dado por la ecuación de estado normalizada en forma matricial:
−( C11R1 + C11R2 ) 1
      1 
d Vc1 C1 R2 Vc1
= 1 1 1 + C1 R1 V
dt Vc2 C2 R2
−( C2 R2
+ C2 R 3
) Vc2 0

Una vez sustituidos unos valores particulares de resistencias y capacidades, el resultado es:
      
d Vc1 −20 10.583 Vc1 10
= + V
dt Vc2 10.583 −11 Vc2 0
Nota: El elemento “10” de la matriz deberı́a ser “9.417” para que realmente le ecuación correspondiera a la representación del circuito eléctrico
expuesto. Dado que en el enunciado del examen aparecı́a el valor 10, se harán los cálculos de la solución con dicho valor.

Se pide:
(a) (0.5 pt.) Escribir la ecuación de salida normalizada del sistema, tomando como variable de
salida y(t) = VR2 , esto es, la caı́da de tensión en la resistencia R2 .
Sol.: Como y = Vc1 − Vc2 , en forma matricial normalizada será:
 

Vc1
y(t) = 1 −1 + 0 · V (t)
Vc2

(b) (1.25 pts.) Mediante la Transformada de Laplace, determinar la función de transferencia que
relaciona la salida Y (s) y la entrada V (s) a partir de las ecuaciones de estado y salida, ası́
como el término de condiciones iniciales de Y (s) en función de las tensiones iniciales de los
condensadores Vc1 (0) y Vc2 (0).
Sol.: Haciendo Transformada de Laplace, tendremos:
        
Vc1 Vc1 (0) −20 10.583 Vc1 10
s − = + V
Vc2 Vc2 (0) 10.583 −11 Vc2 0

con lo que:
      
s + 20 −10.583 Vc1 10 Vc1 (0)
= V +
−10.583 s + 11 Vc2 0 Vc2 (0)

que resulta en:

   −1    
Vc1 s + 20 −10.583 10 Vc1 (0)
= V +
Vc2 −10.583 s + 11 0 Vc2 (0)

y, consecuentemente, la transformada de Laplace de la salida vendrá dada por:

s + 20 −10.583 −1
       

Vc1 10 Vc1 (0)
Y (s) = 1 −1 = 1 −1 V +
Vc2 −10.583 s + 11 0 Vc2 (0)
Sistemas Automáticos (SAU), GITI-ETSII 21 Octubre 2015, pág. 5

Usando la fórmula de la matriz inversa 2 × 2:

     

1 s + 11 10.583 10 Vc1 (0)
Y (s) = 1 −1 · · V +
(s + 20)(s + 11) − 10.5832 10.583 s + 20 0 Vc2 (0)

Operando, resulta en:

10s + 4.17 (s + 0.417)Vc1 (0) − (s + 9.417)Vc2 (0)

Y (s) = ∗V (s) +
s2
| + 31s
{z + 108} | s2 + 31s{z + 108 }
Función de Transferencia Térm. Cond. Iniciales

Una alternativa para obtener de forma “directa” la solución buscada Y (s) resolviendo un
sistema de tres ecuaciones, podrı́a ser plantear y−Vc1 +Vc2 = 0 como una ecuación adicional
a las dos de estado y ponerlas de forma conjunta como:
    
s + 20 −10.583 0 Vc1 (s) 10V (s) + Vc1 (0)
−10.583 s + 11 0 Vc2 (s) =  Vc2 (0)  (4)
−1 1 1 y(s) 0

con lo que la regla de Cramer permitirı́a obtener y(s) sin el paso intermedio de despejar las
tensiones de los condensadores, resolviendo dos determinantes 3 × 3:

s + 20 −10.583 10V (s) + Vc1 (0)

−10.583 s + 11 Vc2 (0)
−1 1 0
y(s) = =
s + 20 −10.583 0
−10.583 s + 11 0
−1 1 1
(10V (s) + Vc1 (0)) ∗ ((s + 11) − 10.583) + Vc2 (0)(−(s + 20) + 10.583)
=
(s + 20)(s + 11) − 10.5832

obteniendose, efectivamente, el mismo resultado como era de esperar.

Nota: Usando Matlab (Control Systems Toolbox), la solución podrı́a haberse obtenido con:
A=[-20 10.583;10.583 -11]; B=[10;0]; C=[1 -1]; D=0;
circuito=ss(A,B,C,D);
disp(’Funcion de transferencia: C*inv(sI-A)*B’)
fdt=tf(circuito)
disp(’Term. Cond iniciales: C*inv(sI-A)’)
tci=tf(ss(A,eye(2),C,D))

(c) (2 pt.) Determinar la respuesta temporal de la salida, y(t), ante una tensión de entrada V (t) =
sin(10t) y condensadores inicialmente descargados.
Sol.: Con los condensadores descargados, tenemos:
10s + 4.17 10s + 4.17 10
Y (s) = ∗V (s) = 2 ∗ 2
s2
| + 31s
{z + 108} |s + 31s
{z + 108} |s +
2
{z10 }
Función de Transferencia Función de Transferencia Entrada

Para realizar la descomposición en fracciones simples, debemos comprobar si las raı́ces del
término s2 + 31s + 108 son reales o complejas. En este caso son reales (roots([1 31 108])):
√
−31 ± 312 − 108 · 4 p
= −17.5 ± 529/4 = {−4, −27}
2
Sistemas Automáticos (SAU), GITI-ETSII 21 Octubre 2015, pág. 6

Además, las raı́ces de s2 + 102 son ±10j (imaginarias puras). Por tanto, la descomposición
en fracciones simples puede ser expresada:

100s + 41.7 M N P · 10 Q·s

Y (s) = 2
= + + 2 2
+ 2
(s + 4)(s + 27)(s + 100) s + 4 s + 27 s + 10 s + 102

Debiendo calcular M , N , P y Q. Pasando a común denominador el lado derecho, tendrı́amos:

M (s + 27)(s2 + 100) + N (s + 4)(s2 + 100) + P 10(s + 4)(s + 27) + Qs(s + 4)(s + 27)
Y (s) =
(s + 4)(s + 27)(s2 + 100)

Con lo que operando, tendrı́amos:

Monomio Num. = Num. descompuesto
s3 : 0 = M +N +Q
s2 : 0 = 27M + 4N + 10P + 31Q
s: 100 = 100M + 100N + 310P + 108Q
1: 41.7 = 2700M + 400N + 1080P
que es un sistema de 4 ecuaciones y 4 incógnitas, que debe ser resuelto.
Nota: para resolver rápidamente el sistema de ecuaciones, se pueden calcular M y N de
forma aislada haciendo s = −4, de modo que igualando numeradores:

100 ∗ (−4) + 41.7 = M ∗ (−4 + 27)((−4)2 + 100) + 0

y haciendo s = −27, se tendrı́a:

100 ∗ (−27) + 41.7 = N ∗ ((−27) + 4)((−27)2 + 100) + 0

pudiéndose despejar P y Q una vez obtenidas M y N . Ası́, tras resolver el sistema de

ecuaciones completamente, el resultado es:

M = −0.134, N = 0.139, P = 0.323, Q = −5.12 · 10−3

que dará la transformada inversa, solución final del apartado:

y(t) = −0.134e−4t + 0.139e−27t + 0.323 sin(10t) + 0.0051 cos(10t)

*Para comprobar con Matlab (symbolic toolbox) el resultado numérico, usaremos los comandos <ilaplace> y <vpa>:
>> syms s
>> vpa(ilaplace( (100*s+41.7)/(s+4)/(s+27)/(sˆ2+10ˆ2) ), 3)
ans =
0.323*sin(10.0*t) - 0.00512*cos(10.0*t) + 0.139*exp(-27.0*t) - 0.134*exp(-4.0*t)
También podemos resolver la EDO y 00 + 31y 0 + 108y = 10u0 + 4.17u con condiciones iniciales nulas usando <dsolve>:
>> syms y(t) t
>> u=sin(10*t); Du=diff(u);
>> Dy=diff(y); D2y=diff(y,2);
>> ysol = dsolve( D2y+31*Dy+108*y==10*Du+4.17*u, y(0)==0, Dy(0)==0);
>> vpa( ysol, 3)
ans =
0.139*exp(-27.0*t) - 0.323*cos(10.0*t + 1.55) - 0.134*exp(-4.0*t)

Nota: Una descomposición en fracciones simples alternativa podrı́a haber sido:

100s + 41.7 M N T W
Y (s) = = + + +
(s + 4)(s + 27)(s2 + 100) s+4 s + 27 s − 10j s + 10j

de modo que los coeficientes T y W hubieran sido el resı́duo en s = 10j y en s = −10j, respectivamente. El resultado,
según Matlab es:
Sistemas Automáticos (SAU), GITI-ETSII 21 Octubre 2015, pág. 7

[k,p]=residue([100 41.7],[1 31 208 3100 10800])

k =
0.139418891278125 + 0.000000000000000i
-0.002561769477143 - 0.161356432760700i
-0.002561769477143 + 0.161356432760700i
-0.134295352323838 + 0.000000000000000i
p =
-26.999999999999993 + 0.000000000000000i
-0.000000000000005 +10.000000000000005i
-0.000000000000005 -10.000000000000005i
-3.999999999999997 + 0.000000000000000i
esto es:
100s + 41.7 0.139 0.134 −0.0256 − 0.161j −0.00256 + 0.161j
Y (s) = = + + +
(s + 4)(s + 27)(s2 + 100) s+4 s + 27 s − 10j s + 10j
que darı́a lugar a

y(t) = −0.134e−4t + 0.139e−27t + (−0.0256 − 0.161j)e10j·t + (−0.0256 + 0.161j)e−10j·t

Para completar el ejercicio, esta solución con números complejos deberı́a ser convertida a expresiones exponenciales y
trigonométricas usando la fórmula de Euler eωj·t = cos(ωt) + j sin(ωt). El resultado serı́a el mismo que antes: en efecto,
T e10j·t + W e−10j·t = T (cos(10t) + j sin(10t)) + W (cos(10t) − j sin(10t)) = (T + W ) cos(10t) + (T − W )j sin(10t)

Se observa que las dos descomposiciones están relacionadas por Q = W + T = 2Re(T ), P = (T − W )j = 2Im(W ).
Obviamente, los términos en T y W ( o en P y Q) serán los que determinen el “régimen estacionario senoidal” una vez las
exponenciales asociadas a M y N hayan desaparecido, que en el caso del circuito eléctrico también podrı́a ser determinado
mediante las adecuadas ecuaciones de impedancias/Kirchoff en corriente alterna estacionaria.

(d) (1.25 pt.) Determinar la respuesta temporal de la salida ante condiciones iniciales nulas y una
tensión de entrada dada por
(
π
5 sin(10t) 0 ≤ t ≤ 10
V (t) = π
0 t > 10

Para este último apartado, se recomienda utilizar los resultados obtenidos en el apartado (c).

Sol.: Nótese que la señal V (t) es una sinusoide activa durante un sólo semiperı́odo (π/10).
Como
sin(10t) = − sin(10t − π), se tiene sin(10t) + sin(10(t − π/10)) = 0, con lo que V (t) se puede
descomponer como:
(
π
0 0 ≤ t ≤ 10
V (t) = 5sin(10t) + π
5sin(10(t − π/10)) t > 10

Denominando y(c) (t) a la solución de y(t) en el apartado (c), por linealidad e invarianza temporal
podemos aplicar superposición y decir que, en el caso que se discute en este apartado, será:
(
π
0 0≤t≤ 10
y(t) = 5y(c) (t) + π
5y(c) (t − π/10) t > 10

Puntuación total: 10 puntos.

Tabla de transformadas:
f (t) F(s) Fórmula matriz inversa 2 × 2
1
1 s

A11 A12

1

A22 −A12

e−at sin(wt) (s+a)w2 +w2 A= ; A−1 =
A21 A22 A11 A22 − A21 A12 −A21 A11
e−at cos(wt) (s+a)
s+a
2 +w 2