Tarea3 08 sol+maqpapel

Departamento de Ingeniería de Sistemas y Automática

Universidad Politécnica de Valencia

INGENIERÍA DE CONTROL I 2008-9

___________________________________________________________________________________________ P. Albertos. Dept. Ingeniería de Sistemas y Automática. UPV. Ext.: 79570 e-mail [email protected] J.V. Roig. Dept. Ingeniería de Sistemas y Automática. UPV. Ext.: 75768 e-mail [email protected]

Tarea 3. Diseño control máquina de papel (entrega 4/12/08) La caja de distribución de una máquina de papel se ha modelado experimentalmente, obteniéndose la matriz de transferencia:

( )( )

( )( )⎥⎦

⎤⎢⎣

⎡⋅

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

+++

+−

++=⎥⎦

⎤⎢⎣

⎡−

−

pf

ssse

ssse

vh

s

s

73.0211.0

71.0086.0

92.1129.0

83.021245.0

2.0

1.0

Existe una fuga constante de material de entrada, w, que no se puede medir. Se pide:

1. Considerando los pares nivel/presión y velocidad/flujo, diseña sendos reguladores que consigan, en ambos bucles de control, amortiguamiento 0.5 y test ≤ π seg, sin error estacionario, independientemente, como si no hubiera interacción entre variables.

2. Determina el comportamiento real del sistema acoplado, con estos reguladores.

3. Se discretiza el modelo con un periodo de muestreo de T = 0.1s. Sin tener en cuenta retardos, obtener una realización mínima y una reducida de 4º orden.

4. Incluyendo retardos, diseña un control centralizado, por realimentación del estado, para conseguir los mismos polos que en el diseño hecho en 1) (introduce acción integral).

a. Con el modelo de realización mínima. b. Con el sistema reducido.

5. Diseña un observador del estado (del sistema reducido) y de la perturbación e implementa un control por realimentación de la salida y con prealimentación.

6. Diseña un control óptimo con respecto al índice:

( )∑∞

=

+++++=0

2,2

2,1

2222 02.0101.05.0k

kkkkkk eaeafpvhJ

suponiendo que se mide el estado y la salida. eai son los errores acumulados de cada salida. Compara las respuestas a escalón en referencias de ambas salidas.

7. Realiza el diseño para conseguir las prestaciones previstas en 1), habiendo hecho previamente un desacoplamiento:

a. con prealimentación dinámica. b. por realimentación del estado.





Solución

1. Diseño de control feedback descentralizado. Utilizando la herramienta sisotool de MATLAB podemos diseñar reguladores que cumplan con las especificaciones del control. Así obtenemos:

- Para la fdt ( ) ( )( ) 92.1

129.012 +

−==ssp

shsg el regulador ( )s

ssK 58.36261.41+−= .

- Para la fdt ( ) ( )( ) 71.0

086.021 +

==ssf

svsg el regulador ( )s

ssK 182.18072.202+= .

2. Simulación del control feedback descentralizado.

El esquema de SIMULINK de dicho control es el mostrado en la siguiente Figura 2, y su evolución se muestra en la Figura 1.

Salida v

Salida h

Ref_v

Ref_h

f

p

h

v

MáquinaPapel

20.8072(s+1.182)

sK2

-4.6261(s+3.58)

sK1

Entradas (f&p)

Figura 2: Control feedback descentralizado.

0 5 10 15 20 250

0.5

1

1.5

velo

cida

d

Control feedback descentralizado

0 5 10 15 20 25-1

0

1

2

3

Tiempo (seg)

nive

l

Figura 1: Salidas con control feedback descentralizado.





El sistema equivalente en bucle cerrado (sin tener en cuenta los retardos) de este esquema lo podemos obtener con las siguientes instrucciones:

>> s = tf('s'); >> g11 = 0.1245/((s+2)*(s+0.83)); >> g12 = -0.129/(s+1.92); >> g21 = 0.086/(s+0.71); >> g22 = 0.11/((s+2)*(s+0.73)); >> G = [g11 g12; g21 g22]; %matriz de transferencia sin retardos >> K1 = -4.6261*(s+3.58)/s; >> K2 = 20.8072*(s+1.182)/s; >> K = [K1 0;0 K2]; %regulador desacoplado >> Gs = series(K, G, [2 1], [1 2]); %conexión regulador-proceso >> H1 = eye(2); %realimentación unitaria >> Gc = feedback(Gs, H1); %Cerrar el bucle

Para determinar el comportamiento del sistema real con estos reguladores podemos calcular los polos de la realización mínimo del sistema equivalente ya calculado con la siguiente instrucción:

>> pc = pole(minreal(ss(Gc))) 12 states removed. pc = -0.4099 + 0.9406i -0.4099 - 0.9406i -1.7798 + 1.2124i -1.7798 - 1.2124i -0.7112 -1.1534 -2.4084 -1.9239

3. Modelo discretizado con T = 0.1s. Sin tener en cuenta retardos, obtener una realización

mínima y una reducida de 4º orden. Con las siguientes instrucciones de MATLAB, podemos discretizar la matriz de transferencia sin retardos y obtener la realización mínima y reducida:

>> T = 0.1; >> Gd = c2d(G, T, 'zoh'); %Matriz de transferencia discreta Transfer function from input 1 to output... 0.0005669 z + 0.0005159 #1: -------------------------- z^3 - 1.739 z^2 + 0.7535 z 0.008302 #2: ---------------- z^3 - 0.9315 z^2 Transfer function from input 2 to output... -0.01174 #1: ---------- z - 0.8253 0.0005026 z + 0.0004589 #2: ----------------------- z^2 - 1.748 z + 0.7611 >> SYSm = minreal(ss(Gd)); %Realización mínima 6 v.e. a = x1 x2 x3 x4 x5 x6 x1 1.739 -0.7535 0 0 0 0





x2 1 0 0 0 0 0 x3 0 0 0.9315 0 0 0 x4 0 0 0 0.8253 0 0 x5 0 0 0 0 1.748 -0.7611 x6 0 0 0 0 1 0 b = u1 u2 x1 0.03125 0 x2 0 0 x3 0.125 0 x4 0 0.125 x5 0 0.03125 x6 0 0 c = x1 x2 x3 x4 x5 x6 y1 0.01814 0.01651 0 -0.0939 0 0 y2 0 0 0.06641 0 0.01608 0.01468 d = u1 u2 y1 0 0 y2 0 0 Sampling time: 0.1 Discrete-time model. >> SYSr = balred(SYSm,4); %Realización reducida de 4º orden a = x1 x2 x3 x4 x1 0.9281 -0.008598 0.001795 0.01681 x2 -0.01835 0.8604 -0.03425 -0.004856 x3 -0.0349 -0.07803 0.8871 -0.008276 x4 -0.04047 0.04588 -0.0278 0.8159 b = u1 u2 x1 -0.1383 -0.02999 x2 0.02833 -0.1432 x3 -0.02983 -0.04511 x4 -0.04884 0.02865 c = x1 x2 x3 x4 y1 -0.01416 0.07964 0.04624 0.0464 y2 -0.07034 -0.01638 0.07233 -0.02404 d = u1 u2 y1 0 0 y2 0 0 Sampling time: 0.1 Discrete-time model.

4. Incluyendo retardos, diseñar un control por realimentación del estado, para conseguir

los mismos polos que en el diseño hecho en 1) (introduce acción integral). En primer lugar discretizamos los polos continuos obtenidos en el apartado 2) para conocer los polos que debemos asignar:

>> pd=exp(T*pc) pd = 0.9556 + 0.0901i





0.9556 - 0.0901i 0.8308 + 0.1012i 0.8308 - 0.1012i 0.9313 0.8911 0.7860 0.8250

a. Con el modelo de la realización mínima.

A la matriz de transferencia discreta Gd(z) le añadimos los retardos, obtenemos su realización mínima (que tendrá 8 variables de estado: las 6 que tenía la realización mínima más los dos retardos que hemos añadidos).

>> Gd(1,1) = Gd(1,1)*tf(1,[1 0],T); %Se añade 1 retardo a g11 >> Gd(2,1) = Gd(2,1)*tf(1,[1 0 0],T); %Se añaden 2 retardos a g21 >> SYSmr = minreal(ss(Gd)); %realización mínima 8 v.e.

Para incluir la acción integral en la realimentación del estado, añadimos el error acumulado como variable de estado (ea1 y ea2):

( ) ( )( )

( )( ) ( ) ( )kr

Iku

DB

keakx

ICA

keakx

kxa ⋅⎥⎦

⎤⎢⎣

⎡+⋅⎥

⎦

⎤⎢⎣

⎡−

+⎥⎦

⎤⎢⎣

⎡⋅⎥

⎦

⎤⎢⎣

⎡−

=⎥⎦

⎤⎢⎣

⎡++

=+00

11

1

>> Aa =[SYSmr.A zeros(8,2); -SYSmr.C eye(2)]; %Matriz A ampliada >> Ba = [SYSmr.B; -SYSmr.D]; %Matriz B ampliada

El sistema ampliado tendrá 10 variables de estado; podemos asignar los 8 polos que tendría el sistema discreto equivalente del apartado 1), y por ejemplo 2 polos en 0.

>> Ka = place(Aa, Ba, [pd' 0 0]); El esquema de SIMULINK de dicho control se muestra en la Figura 3 y su evolución en la Figura 4.

ea1

ea2

xax

y1y2

ea

ref_y2

ref_y1 ModeloMínimo

K*u

Ka

y(n)=Cx(n)+Du(n)x(n+1)=Ax(n)+Bu(n)

Discrete State-SpaceModelo mínimo

z -1

1-z -1

Discrete Filter1

z -1

1-z -1

Discrete Filter

em

Figura 3: Realimentación del estado con control integral.





b. Con el modelo del sistema reducido.

De la realización mínima con retardos obtenemos una realización reducida de 6º orden (4 variables de estado de la realización reducida inicial más 2 de los retardos).

>> SYS4b = balred(SYSmr, 6);

De forma similar al apartado anterior, incluimos el error acumulado de cada salida como variable de estado adicional en un sistema ampliado para el cual calculamos la realimentación del estado que asigne 6 de los polos que tiene el sistema discreto equivalente del apartado 1):

>> Aa2 = [SYS4b.A zeros(6,2); -SYS4b.C eye(2)]; % A ampliada >> Ba2 = [SYS4b.B; -SYS4b.D]; % B ampliada >> Ka2 = place(Aa2, Ba2 , pd);

Con un esquema de SIMULINK similar al de la Figura 3 anterior pero realimentando sólo 6 variables de estado, se obtiene la evolución de las salidas mostrada en la Figura 5.

0 5 10 15 20 250

0.5

1

1.5

velo

cida

d

Control centralizado+control integral (Sistema mínimo)

0 5 10 15 20 25-1

0

1

2

3

Tiempo (seg)

nive

l

Figura 4: Salidas con realimentación del estado + control integral.





5. Diseño de observador del estado (del sistema reducido) y de la perturbación, y diseño

del control por realimentación de la salida y con prealimentación. Al sistema discreto reducido con retardos incluidos (SYS4b) le añadimos la perturbación, w; (consiste en cambiar la entrada f por f – w).

( ) ( ) [ ] ( )( )⎥⎦

⎤⎢⎣

⎡⋅+⋅=+

kwku

FBkxAkx 1

>> Ap = SYS4b.A; >> F = -SYS4b.B(:,1); >> Bp = [SYS4b.B F]; %añadimos la entrada: w = -u1 >> Cp = SYS4b.C; >> Dp = [SYS4b.D -SYS4b.D(:,1)] >> SYS5 = ss(Ap, Bp, Cp, Dp);

Para estimar la perturbación, la añadimos como una variable de estado más al sistema a observar, y diseñamos un observador de orden completo para dicho sistema:

( ) ( )( )

( )( ) ( )ku

BkwkxFA

kwkx

kx p ⋅⎥⎦

⎤⎢⎣

⎡+⎥

⎦

⎤⎢⎣

⎡⋅⎥

⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡++

=+0101

11

>> A5p = [SYS5.A SYS5.B(:,3); zeros(1,6) 1]; >> B5p = [SYS5.B(:,1:2); zeros(1,2)]; >> C5p = [SYS5.C zeros(2,1)]; >> D5p = SYS5.D(:,1:2) >> S5p = ss(A5p, B5p, C5p, D5p); >> K0 = place(S5p.A',(S5p.C*S5p.A)',[0 0 0.1 0.2 0.3 0.25 0.21])';

Si se conoce la perturbación, se puede compensar su efecto sobre la salida mediante una prealimentación estática:

( ) ( )FBPre

kwPreku⋅−=

⋅=−1

ˆ

Figura 5: Salidas con realimentación del estado + control integral.

0 5 10 15 20 250

1

2

3

4

velo

cida

d

Control centralizado+control integral (Sistema reducido)

0 5 10 15 20 25-1

0

1

2

3

Tiempo (seg)

nive

l





La acción de control, por tanto estará formada por una prealimentación de la perturbación estimada más la realimentación (Ka2) del estado observado que ya hemos diseñado en el apartado 4b) y que incluía también la acción integral. El esquema de SIMULINK de este control lo podemos ver en la Figura 7, y la evolución de las salidas en la Figura 6.

0 5 10 15 20 25-1

0

1

2

3

4

velo

cida

d

Prealimentación+Realimentacion estado observado+Control integral (Sistema reducido)

0 5 10 15 20 25-1

0

1

2

3

Tiempo (seg)

nive

l

Figura 6: Salidas con prealimentación+realimentación salida con control integral.

Figura 7: Prealimentación más realimentación salida con control integral.

xp_ex_e

w_e

eaxa

w

ref_y2

ref_y1 Salidas

K*u

Pre


ObservadorActualizado

K*u

Ka2


Discrete State-SpaceModelo reducidocon perturbación

z -1

1-z -1

Discrete Filter1

z -1

1-z -1

Discrete Filter

em

em





6. Control óptimo. En primer lugar expresaremos el índice en forma matricial, y lo adaptaremos para poder aplicar la ecuación de Ricatti, teniendo en cuenta que en el estado del sistema ampliado del apartado 4) hemos incluido el error acumulado como variable de estado.

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( )[ ] ( )( ) ( ) ( )

( ) ( ) ( ) ( )∑

∑

∑

∑

∑

∞

=

∞

=

∞

=

∞

=

∞

=

⋅⋅+⋅⋅=

⋅⋅+⎥⎦

⎤⎢⎣

⎡⋅⎥

⎦

⎤⎢⎣

⎡⋅=

⋅⋅+⋅⋅+⋅⋅⋅⋅=

⋅⋅+⋅⋅+⋅⋅=

⋅⎥⎦

⎤⎢⎣

⎡⋅+⋅⎥

⎦

⎤⎢⎣

⎡⋅+⋅⎥

⎦

⎤⎢⎣

⎡⋅=

0

0

0

0

0

00

02.00002.0

1.00010

5.0001

k

Taa

Ta

k

T

ea

kea

TTy

TT

kea

TTy

T

k

TTT

kuRkukxQkxJ

kuRkukeakx

QQ

keakxJ

keaQkeakuRkukxCQCkxJ

keaQkeakuRkukyQkyJ

keakeakukukykyJ

Para el sistema reducido y ampliado del apartado 4b) diseñamos una realimentación del estado óptima para el sistema ampliado con el comando dlqr de MATLAB:

Qy = [1 0; 0 0.5]; %poderación de las salidas. R = [10 0; 0 0.1]; %ponderación de las entradas Qea = 0.02*[1 0; 0 1]; %ponderación del error acumulado Q = C'*Qy*C; %ponderación del estado x Qa = [Q zeros(6,2); zeros(2,6), Qea]; %ponderación del est ampliado Kopt = dlqr(Aa2, Ba2, Qa, R);

El esquema de SIMULINK es el mismo utilizado en el apartado 4b) de la Figura 3, pero sustituyendo la matriz de realimentación del estado por Kopt. La evolución de las salidas con este control es la mostrada en la Figura 8.

0 10 20 30 40 50 60 700

0.2

0.4

0.6

0.8

1

velo

cida

d

Control óptimo (Sistema reducido)

0 10 20 30 40 50 60 70-0.5

0

0.5

1

1.5

2

Tiempo (seg)

nive

l

Figura 8: Salidas con control óptimo.





7. Control desacoplado.

El desacoplamiento lo podemos realizar de dos formas: - Mediante una prealimentación dinámica, a partir de un modelo de matriz de

transferencia. - Mediante una realimentación del estado, a partir de un modelo en espacio de

estados.

a. Desacoplamiento con prealimentación dinámica. Para el modelo de matriz de transferencia continua sin retardos G(s), diseñamos una prealimentación dinámica que consiga desacoplar el sistema con los pares nivel/presión (y1/u2) y velocidad/flujo(y2/u1). Una posibilidad es descomponer la matriz de transferencia del sistema en una matriz diagonal D(s) (formado por el denominador común de cada fila), y otra matriz H(s) cuya inversa sea realizable y que utilizaremos para desacoplar:

( ) ( )( )

( )( )( ) ( )sHsD

sss

ssssG ⋅=

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

+++

+−

++=

73.0211.0

71.0086.0

92.1129.0

83.021245.0

( ) ( )( )( )

( )( )( ) ⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

+++

+++=0

71.073.021

92.183.0210

sss

ssssD

( ) ( )( ) ( )( ) ( )( )⎥⎦

⎤⎢⎣

⎡++−+

+++=

83.02129.092.11245.071.011.073.02086.0sss

ssssH

El sistema desacoplado ( ) ( ) ( )sDsHsG =⋅ −1 , está compuesto por dos f.d.t. independientes, para las cuales utilizando la herramienta sisotool de MATLAB podemos diseñar sus respectivos reguladores que cumplan con las especificaciones del control. Así obtenemos:

- Para la fdt ( )sd12 el regulador: ( ) ( )( )( )96.2991.10

83.0131.4717.373.46 2

2

12 +++++=

sssssssKd

- Para la fdt ( )sd 21 el regulador: ( ) ( )( )( )11.6841.16

71.0696.2217.347.114 2

2

21 +++++=

sssssssKd

El sistema equivalente en bucle cerrado, de este control por desacoplamiento por prealimentación dinámica y su simulación lo podemos obtener con las siguientes instrucciones:

>> Ds = [0 1/((s+2)*(s+0.83)*(s+1.92)); 1/((s+0.71)*(s+2)*(s+0.73)) 0]; >> Hs = [0.086*(s+2)*(s+0.73) 0.11*(s+0.71); 0.1245*(s+1.92) -0.129*(s+2)*(s+0.83)]; >> Grd = inv(Hs); >> Kd12 = 46.73*(s^2+3.71*s+4.13)*(s+0.83)/(s*(s^2+10.91*s+29.96)); >> Kd21 = 114.47*(s^2+3.21*s+2.7)*(s+0.71)/(s*(s^2+16.41*s+68.11));





>> Kd = [0 Kd12; Kd21 0]; >> Gs = zpk(minreal(G*Grd*Kd), 0.01); >> H1 = eye(2); %realimentación unitaria >> Gc = feedback(Gs, H1); %Cerrar el bucle >> step(Gc);

El esquema de SIMULINK de este control se puede ver en la y la evolución de las salidas en la.

b. Desacoplamiento por realimentación del estado.

El desacoplamiento por realimentación del estado lo aplicaremos al modelo discreto reducido del apartado 4b (SYS4b). En primer lugar calculamos la condición de Gilbert para cada salida:

>> J1 = SYS4b.C(1,:)*SYS4b.A^0*SYS4b.B = [-0.0000 -0.0117] >> J2 = SYS4b.C(2,:)*SYS4b.A^0*SYS4b.B = 1e-3 *[-0.0000 0.5024]

Ambos filas son no nulas y la matriz J es invertible:

Salida v

Salida h

Ref_v

Ref_hf

p

h

v

MáquinaPapel

mat(Kd21.nu

2mat(Kd21.deKd21

mat(Kd12.nu

2mat(Kd12.deKd1

r1

r2

u1

u2

H(s)^-1

Figura 10: Control por desacoplamiento por prealimentación dinámica.

0 5 10 15 20 25-0.5

0

0.5

1

1.5

2

velo

cida

d

Control Desacoplado por Prealimentación

0 5 10 15 20 250

0.5

1

1.5

2

2.5

Tiempo (seg)

nive

l

Figura 9: Salidas con control por desacoplamiento por prealimentación.





>> J = [J1; J2]; >> rank(J) = 2

Por tanto el sistema es desacoplable mediante una realimentación del estado y una prealimentación: ( ) ( ) ( )krFkxKdesku ⋅+⋅= , que se obtienen de la siguiente forma:

>> At = [SYS4b.C(1,:)*SYS4b.A^1; SYS4b.C(2,:)*SYS4b.A^1]; >> Kdes = -inv(J)*At; Kdes = 1.0e+008 * 2.6295 -0.8365 -3.2218 0.3190 5.4074 -0.8875 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 >> F = inv(J); F = 1.0e+009 * -0.1852 -4.3256 -0.0000 0.0000

El sistema equivalente en bucle cerrado tendrá un polo en 0 (1 retardo) en cada uno de los bucles desacoplados:

>> Ad = SYS4b.A+SYS4b.B*Kdes; >> Bd = SYS4b.B*F; >> SYS7b = ss(Ad, Bd, SYS4b.C, SYS4b.D, T); >> Gd7b = minreal(tf(SYS7b), 0.01);

( )⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

=

z

zsbGd 10

01

7

Tarea3 08 sol+maqpapel

Technology

Transcript of Tarea3 08 sol+maqpapel