```
%matplotlib inline
from matplotlib.pyplot import plot, figure, gca, legend
```


```
from control import tf, rlocus, ss, ss2tf, step_response
from numpy import linspace, matrix, roots
from sympy import init_printing, Matrix, var, eye
init_printing()
```

# Comparación de RL contra Retroalimentación de Estado

## Root Locus

Tenemos el siguiente sistema retroalimentado:

$$
G(s) = \frac{4}{s(s+2)}
$$


```
G = tf([0, 0, 1],[1, 2, 0])
```


```
figure(figsize=(10, 10))
r, k = rlocus(G, linspace(0, 100, 1000))
ax = gca()
ax.set_xlim(-2.5, 0.5)
ax.set_ylim(-3.0, 3.0);
```

Tenemos como parametros deseados $\zeta = 0.5$ y $\omega_n = 4$, para los cuales habiamos considerado anteriormente el siguiente compensador:

$$
G_C(s) = \frac{s + 2.9}{s + 5.4}
$$


```
Gc = tf([0, 1, 2.9], [0, 1, 5.4])
K = tf([0, 0, 4], [0, 0, 1])
```


```
figure(figsize=(10, 10))
r, k = rlocus(Gc*K*G, linspace(0, 100, 1000))
ax = gca()
ax.set_xlim(-6.5, 0.5)
ax.set_ylim(-4.0, 4.0);
```

Y graficando ambos juntos:


```
figure(figsize=(10, 10))
r, k = rlocus(G, linspace(0, 1000, 10000))
r, k = rlocus(Gc*K*G, linspace(0, 1000, 10000))
ax = gca()
ax.set_xlim(-6.5, 0.5)
ax.set_ylim(-4.0, 4.0);
```


```
print G*Gc*K
```

    
          4 s + 11.6
    ----------------------
    s^3 + 7.4 s^2 + 10.8 s
    


Y si quieremos obtener la respuesta ante una señal tipo escalon tendremos:


```
Gf = G.feedback(tf([1],[1]))
Gtf = (G*Gc*K).feedback(tf([1],[1]))
```


```
tf, yf = step_response(Gf, T=linspace(0, 10, 1000))
tt, yt = step_response(Gtf, T=linspace(0, 10, 1000))
```


```
figure(figsize=(10, 10))
p1, = plot(tf, yf)
p3, = plot(tt, yt)
legend([p1, p3], [r"$G_F$", "$k G_c G_F$"]);
```

## Representación en Espacio de Estados

Ahora, para analizar el mismo sistema con su representación de estado, primero lo tenemos que llevar a esta representación:

Primero notamos que el grado de numerador es $m = 0$ y que el grado del denominador es $n = 2$, por lo que, de acuerdo a nuestra forma normal para las Ecuaciones Diferenciales Ordinarias, tenemos que los coeficientes de la EDO son:

$$
\frac{y(s)}{u(s)} =
\frac{b_0 s^m + b_1 s^{m-1} + ... + b_{m-1} s + b_m}{s^n + a_1 s^{n-1} + ... + a_{n-1} + a_n}
$$

Si $G(s) = \frac{4}{s^2 + 2s}$, entonces:

$$
b_0 = 4
$$

$$
a_1 = 2 , \quad a_2 = 0
$$

Por lo que tenemos que nuestra representación de estado será la siguiente:

$$
\frac{d}{dt}
\begin{pmatrix}
    x_1 \\
    x_2
\end{pmatrix} = 
\begin{pmatrix}
    0 & 1 \\
    0 & -2
\end{pmatrix}
\begin{pmatrix}
    x_1 \\
    x_2
\end{pmatrix} +
\begin{pmatrix}
    0 \\
    1
\end{pmatrix} u
$$

$$
y =
\begin{pmatrix}
    4 & 0 
\end{pmatrix}
\begin{pmatrix}
    x_1 \\
    x_2
\end{pmatrix}
$$


```
var("s")
```


```
A = Matrix([[0, 1], [0, -2]])
b = Matrix([[0],[1]])
c = Matrix([[4], [0]])
d = 0
```

Podemos obtener su función de transferencia con la siguiente formula:

$$
F(s) = c^T (sI - A)^{-1} b
$$


```
G = c.T*(s*eye(2) - A).inv()*b
```


```
G[0]
```

Y obtenemos su polinomio caracteristico con:

$$
P(s) = |s I - A|
$$

y sus polos con los eigenvalores de $A$.


```
Pc = (s*eye(2) - A).det()
Pc
```


```
A.eigenvals()
```

Los ceros del sistema se pueden obtener con el determinante de la matriz sistema.

$$
\Sigma (s) =
\begin{pmatrix}
    sI - A & b \\
    -c^T & d
\end{pmatrix}
$$


```
S =((s*eye(2) - A).row_join(b)).col_join((-c.T).row_join(Matrix([d])))
S
```


```
S.det()
```

### Retroalimentación de Estado

Para obtener el comportamiento deseado (coeficiente de amortiguamiento $\zeta = 0.5$ y frecuencia natural de $\omega_n = 4$), tenemos la siguiente función de transferencia.

$$G_C = \frac{1}{s^2 + 2 \zeta \omega_n s + \omega_n^2}$$

por lo que podemos calcular la función de transferencia deseada facilmente:


```
zeta = 0.5
omega = 4
1/(s**2 + 2*zeta*omega*s + omega**2)
```

Por lo que sabemos, que los coeficientes deseados para $a_1$ y $a_2$ son $4$ y $16$. Pero recordemos que en el sistema tenemos:

$$a_1 = 2 \quad y \quad a_2 = 0$$

por lo que nuestro vector de retroalimentación $f$ será:

$$
f =
\begin{pmatrix}
    0 - 16 \\
    2 - 4
\end{pmatrix} =
\begin{pmatrix}
    -16 \\
    -2
\end{pmatrix}
$$


```
f = Matrix([[-16], [-2]])
f
```

Y ahora tan solo tenemos que aplicar lo siguiente:

$$
A_f = A + b f^T
$$


```
Af = A + b*f.T
Af
```

Y de nuevo podemos obtener sus polos con los eigenvalores de $A_f$:


```
Af.eigenvals()
```

Su función de transferencia:


```
Gc = c.T*(s*eye(2) - Af).inv()*b
Gc
```

Tenemos que calcular tambien $c_f^T = c^T + d f^T$ y $b$ y $d$, permanecen igual.


```
cf = c + d*f
cf
```

Y ahora integramos la matriz sistema y podemos obtener sus ceros:


```
Sf =((s*eye(2) - Af).row_join(b)).col_join((-cf.T).row_join(Matrix([d])))
Sf
```


```
Sf.det()
```

Si queremos gráficar el comportamiento de este sistema con retroalimentación de estado con el escalon unitario, tenemos que declarar el sistema de espacio de estados:


```
# Volvemos a definir las matrices, pero como arreglos de numpy
# debido a que el modulo de control no sabe que hacer con las
# variables simbolicas
A = matrix([[0, 1], [0, -2]])
b = matrix([[0], [1]])
c = matrix([[4], [0]])
d = 0
f = matrix([[-16], [-2]])
Af = A + b*f.T
cf = c + d*f
# Se agrega una ganancia de 4 a la entrada del sistema, debido
# a que crece la norma del denominador, y por la tanto se atenua
# la salida del sistema.
SR = ss(Af, 4*b, cf.T, d)
```


```
tr, yr = step_response(SR, T=linspace(0, 10, 1000))
```


```
figure(figsize=(15, 8))
p1, = plot(tf, yf)
p2, = plot(tr, yr)
p3, = plot(tt, yt)
legend([p1, p2, p3], [r"$G_F$", r"$SS$", "$k G_c G_F$"]);
```

### Observador de estado

No siempre es posible tener acceso al estado del sistema, por lo que es importante saber como estimar el valor del estado. Esto se hace por medio de un observador.

Un observador de estado, basicamente es una copia de la planta del sistema, la cual alimentamos con la salida del sistema original $y(s)$ y nos devuelve el estado del sistema $x(s)$.


```

```

13614599ff2b6b45c6febc6b60d7f982da22a3f7
ipynb
Jupyter Notebook
IPythonNotebooks/Propedeutico/Clase 10.ipynb
chelizalde/DCA
34fd4d500117a9c0a75b979b8b0f121c1992b9dc
MIT
IPythonNotebooks/Propedeutico/Clase 10.ipynb
chelizalde/DCA
34fd4d500117a9c0a75b979b8b0f121c1992b9dc
MIT
IPythonNotebooks/Propedeutico/Clase 10.ipynb
chelizalde/DCA
34fd4d500117a9c0a75b979b8b0f121c1992b9dc
MIT
2021-03-20T12:44:13.000Z
2021-03-20T12:44:13.000Z
Qwen/Qwen-72B
1. YES
2. YES