```python
import sympy 
sympy.init_printing()
%matplotlib inline

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import control


import scipy.optimize
from tbcontrol.responses import fopdt, sopdt
from ipywidgets import interact

def resultplot1(K, tau, theta, y0):
    plt.figure(figsize=(21, 7))
    plt.scatter(ts, ym)
    plt.plot(ts, fopdt(ts, K, tau, theta, y0), color='red')
    plt.grid(True)

def resultplot2(K, tau, zeta, theta, y0):
    plt.figure(figsize=(21, 7))
    plt.scatter(ts, ym)
    plt.plot(ts, sopdt(ts, K, tau, zeta,theta, y0), color='red')
    plt.grid(True)
    

```

# PID

El controlador PID se compone de:
- Una acción Proporcional para brindar velocidad de respuesta y disminuir el error de estado estable
- Una acción Integral para eliminar el error de estado estable
- Una acción Derivativa para encargarse de errores cuando hay un cambio en el valor absoluto del error teniendo en cuenta que si el error es constante, solo actúan los modos proporcional e integral. 

La función de transferencia del PID es:

$$K_C\left ( 1 + \frac{1}{\tau_I s}+ \tau_D s \right )$$

La forma anterior no es realizable debido a la acción derivativa. Por esta razon, es común usar la forma:

$$K_C\left ( 1 + \frac{1}{\tau_I s}+ \frac{\tau_D s}{\alpha \tau_I s+1} \right )$$

donde se escoge $\alpha$ de manera que la parte derivativa sea una "buena" aproximación de la derivada.



## Diseño directo

Esta técnica parte del supuesto de definir el comportamiento deseado para el sistema en lazo cerrado, escribir una función de transferencia de lazo cerrado de una forma apropiada y despejar el compensador



Partiendo del diagrama anterior, se define la relación deseada para $\frac{Y}{Y_{SP}}$ considerando $D=0$. Después se despeja $G_C$.


```python
s, G_C = sympy.symbols('s, G_C')
G_C
```

Se parte de un modelo de primer orden más tiempo muerto como expresión deseada de $\frac{Y}{Y_{SP}}$.


```python
tau_c, phi = sympy.symbols('tau_c, phi', positive=True, nonzero=True)
desired_Y_over_Y_sp = sympy.exp(phi*s)/(tau_c*s + 1)
desired_Y_over_Y_sp
```

Recuerde la forma de la respuesta del sistema ante un escalón.


```python
from ipywidgets import interact
t = sympy.Symbol('t', positive=True)
def plotresponse(theta=(0, 3.), tau_c_in=(1., 5.)):
    desired_response = sympy.inverse_laplace_transform(desired_Y_over_Y_sp.subs({phi: -theta, tau_c: tau_c_in})/s, s, t)
    p = sympy.plot(desired_response, (t, 0, 10), show=False)
    p2 = sympy.plot(1, (t, 0, 10), show=False)
    p.append(p2[0])
    p.show()
interact(plotresponse);
```


    interactive(children=(FloatSlider(value=1.5, description='theta', max=3.0), FloatSlider(value=3.0, description…


A partir de un modelo aproximado para el proceso $\widetilde{G_P}$, puede calcularse la función de transferencia de lazo cerrado.


```python
GPtilde = sympy.Symbol(r'\widetilde{G_P}')
actual_Y_over_Y_sp = GPtilde*G_C/(1 + GPtilde*G_C)
actual_Y_over_Y_sp
```

En este punto, podría encontrarse el compensador al resolver la igualdad con la función de transferencia  deseada para el sistema en lazo cerrado.


```python
G_C_solved, = sympy.solve(desired_Y_over_Y_sp - actual_Y_over_Y_sp, G_C)
G_C_solved
```

Esta función de transferencia tiene $e^{\phi s}$ en el numerador y en el denominador.  
- El significado de la exponencial en el numerador es de desplazamiento temporal (tiempo muerto).
- El denominador se aproxima mediante una expansión de serie de Taylor de primer orden.


```python
denom = sympy.denom(G_C_solved)
G_C_rational = G_C_solved*denom/denom.subs(sympy.exp(phi*s), 1 + phi*s)
G_C_rational
```

Se define la forma para el compensador PID

$$K_C\left ( 1 + \frac{1}{\tau_I s}  + \tau_D s      \right )$$


```python
K_C, tau_I, tau_D = sympy.symbols('K_C, tau_I, tau_D', positive=True, nonzero=True)
PID = K_C*(1 + 1/(tau_I*s) + tau_D*s)
PID.expand().together()
```

Aunque la forma anterior no es realizable debido a la acción derivativa. Por esta razon, es común usar la forma:

$$K_C\left ( 1 + \frac{1}{\tau_I s}+ \frac{\tau_D s}{\alpha \tau_I s+1} \right )$$


```python
alpha = sympy.symbols('alpha')
ISA = K_C*(1 + 1/(tau_I*s) + tau_D*s/(alpha*tau_D*s + 1))
ISA.expand().together().simplify()
```

El objetivo es hallar los parámetros del compensador para satisfacer la expresión racional $G_C$


```python
G_C_rational
```

Se requiere un modelo de proceso $G$. En este ejemplo se usará un modelo de segundo orden con tiempo muerto.


```python
K, tau, zeta, phi,= sympy.symbols('K, tau, zeta, phi', positive=True)
G2 = K*sympy.exp(phi*s)/((tau*s)**2 + 2*tau*zeta*s + 1)
G2
```

Se incorpora el modelo del sistema a la expresión del compensador.


```python
target_G_C = G_C_rational.subs(GPtilde, G2).expand().together()
target_G_C
```

Se buscan los parámetros del PID que conformen el compensador despejado.


```python
# diferencia entre compensadores determinados
zeroeq = (target_G_C - PID).simplify()
zeroeq
```


```python
numer, denom = zeroeq.as_numer_denom()
eq = sympy.poly(numer, s)
eqs = eq.coeffs()
sympy.solve(eqs, [K_C, tau_D, tau_I])
```

La solución obtenida muestra:

$$K_C = -\frac{2\tau \zeta}{K(\phi-\tau_C)}$$
$$\tau_D = \frac{\tau}{2\zeta}$$
$$\tau_I = 2\tau\zeta$$


Los parámetros del compensador quedan completamente definidos a partir de un modelo aproximado del proceso

$$\frac{K}{(\tau s)^2 + 2\tau \zeta s + 1} \cdot e^{s\phi}$$

Y una respuesta deseada

$$\frac{e^{\phi s}}{\tau_c s + 1}$$

**Ejemplo**

Ajuste un controlador PID para un sistema que fue excitado con un escalón unitario y su respuesta fue registrada.


```python
GP = control.tf([1, 2], [2, 3, 4, 1])
GP
df = pd.read_csv('Datos\DatoEscalon.csv')
df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Unnamed: 0</th>
      <th>tiempo</th>
      <th>U</th>
      <th>Y</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0.000000</td>
      <td>1.0</td>
      <td>5.011953</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>0.196392</td>
      <td>1.0</td>
      <td>5.050839</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2</td>
      <td>0.392785</td>
      <td>1.0</td>
      <td>5.101922</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3</td>
      <td>0.589177</td>
      <td>1.0</td>
      <td>4.996570</td>
    </tr>
    <tr>
      <th>4</th>
      <td>4</td>
      <td>0.785569</td>
      <td>1.0</td>
      <td>5.060810</td>
    </tr>
  </tbody>
</table>
</div>




```python
ym = df.Y
ts = df.tiempo
interact(resultplot2, K=(0.0, 10.0), tau=(0., 10.), 
         theta=(0., 10.), zeta = (0, 10.0), y0=(0., 10.));
```


    interactive(children=(FloatSlider(value=5.0, description='K', max=10.0), FloatSlider(value=5.0, description='t…



```python
p_inicial = [2, 1, 1.1, 0.5, 4.9]
[K_1, tau_1, zeta_1, theta_1, y0_1], pcov = scipy.optimize.curve_fit(sopdt, ts, ym, p_inicial)
[K_1, tau_1, zeta_1, theta_1, y0_1]
```

Se define el comportamiento deseado.


```python
tauC = 1.2
```

Se calculan los parámetros del PID 


```python
sympy.solve(eqs, [K_C, tau_D, tau_I])
```


```python
KC1 = -2*tau_1*zeta_1/(K_1*(-theta_1 - tauC))
TD1 = tau_1/(2*zeta_1)
TI1 = 2*tau_1*zeta_1
[KC1,TI1,TD1]
```


```python
PID1 = KC1*(1 + 1/(TI1*s) + TD1*s)
PID1.expand().together()
```


```python
GC_PID1 = KC1*(1 + control.tf(1,[TI1,0]) + control.tf([TD1,0],1))
GC_PID1
```




$$\frac{0.02943 s^2 + 1.931 s + 0.7193}{2.685 s}$$




```python
_, yLA = control.step_response(GP,ts)
_, yLC1 = control.step_response(control.feedback(GP,1),ts)
_, yLCPID1 = control.step_response(control.feedback(GC_PID1*GP,1),ts)

plt.plot(ts,yLA,
         ts,yLC1,
         ts,yLCPID1)
plt.legend(('Lazo abierto','Realimentado','con PID directo'))
plt.grid()

```

Observe que se presentan oscilaciones aunque la respuesta deseada se modeló como un sistema de primer orden con tiempo muerto. 

- ¿A qué se debe esta discrepancia?
- ¿Qué puntos débiles puede tener esta metodología?

------

## Reglas de sintonía

[Ziegler y Nichols](https://controlautomaticoeducacion.com/control-realimentado/ziegler-nichols-sintonia-de-control-pid/) propusieron métodos para sintonizar un controlador PID.

### En lazo abierto
El método está desarrollado para sistemas cuya respuesta ante un escalón unitario tienen forma de $S$., es decir, sin sobreimpulsos ni dinámicas inestables. En este caso, el sistema puede aproximarse mediante un modelo de primer orden y tiempo muerto:

$$\frac{K}{\tau s + 1}e^{-Ls}$$


como se muestra en la figura.



A partir de este modelo, se sugieren los siguientes valores para los compensadores:

| Controlador | $K_C$               | $\tau_I$       | $\tau_D$ |
|-------------|---------------------|----------------|----------|
| P           |$\frac{\tau}{KL}$    | $\infty$       | $0$      |
| PI          |$0.9\frac{\tau}{KL}$ |$\frac{L}{0.3}$ | $0$      |
| PID         |$1.2\frac{\tau}{KL}$ |$2L$            | $0.5L$   |



```python
ym = df.Y
ts = df.tiempo
interact(resultplot1, K=(0.0, 10.0), tau=(0., 10.), 
         theta=(0., 10.), y0=(0., 10.));
```


    interactive(children=(FloatSlider(value=5.0, description='K', max=10.0), FloatSlider(value=5.0, description='t…



```python
p_inicial = [1.8, 2.0, 1.0, 5.0]
[K_2, tau_2, theta_2, y0_2], pcov = scipy.optimize.curve_fit(fopdt, ts, ym, p_inicial)
[K_2, tau_2, theta_2, y0_2]
```

A partir del modelo, se calculan los parámetros del PID.

$$K_C = 1.2 \frac{\tau}{KL}$$
$$\tau_I = 2 L$$
$$\tau_D = 0.5 L$$



```python
KC2 = 1.2*tau_2/(K_2*theta_2)
TI2 = 2*theta_2
TD2 = 0.5*theta_2
[KC2,TI2,TD2]
```


```python
GC_PID2 = KC2*(1 + control.tf(1,[TI2,0]) + control.tf([TD2,0],[0.001*TD2,1]))
GC_PID2
```




$$\frac{1.243 s^2 + 3.319 s + 2.217}{0.0005604 s^2 + 1.497 s}$$




```python

_, yLCPID2 = control.step_response(control.feedback(GC_PID2*GP,1),ts)

plt.plot(ts,yLA,
         ts,yLC1,
         ts,yLCPID1,
         ts,yLCPID2)
plt.legend(('Lazo abierto','Realimentado','PID directo','ZyN abierto'))
plt.grid()

```

### En lazo cerrado
El método está desarrollado para sistemas que lleguen a oscilación con al ser realimentados. 

Se cierra el lazo eliminando las acciones derivativas e integradoras, es decir, se configura 

$$\tau_I = \infty$$ 
$$\tau_D = 0$$ 

Se varía la ganancia proporcional hasta obtener una respuesta en lazo cerrado que sea una oscilación de amplitud constante. 

- La ganancia del compensador que logra esta respuesta se conoce como ganancia crítica $K_u$.
- El periodo de dicha oscilación es el periodo crítico $P_u$.

Con estos valores se encuentran los parámetros del PID

| Controlador   |	$K_p$   |	 $\tau_I$       |  $\tau_d$     |
|---------------|-----------|-------------------|---------------|
| P	            | $0.5K_u$  | 	$\infty$        |	$0$         |
| PI	        | $0.45K_u$ |$\frac{1}{1.2}P_u$ | 	$0 $        |
| PID	        | $0.6K_u$  | $0.5P_u$          |  $0.125P_u$     |




```python
GP
```




$$\frac{s + 2}{2 s^3 + 3 s^2 + 4 s + 1}$$




```python
def resultplotKfeed(K):
    plt.figure(figsize=(21, 7))
    ts = np.linspace(0,5,1000)
    _, yLCPIDk = control.step_response(control.feedback(K*GP,1),ts)
    plt.plot(ts,yLCPIDk)
    plt.grid(True)
    
interact(resultplotKfeed, K=(0.0, 20.0))
```


    interactive(children=(FloatSlider(value=10.0, description='K', max=20.0), Output()), _dom_classes=('widget-int…





    <function __main__.resultplotKfeed(K)>




```python
Ku = 9.4
Pu = 2.4
KC3 = 0.6*Ku
TI3 = 0.5*Pu
TD3 = 0.125*Pu
GC_PID3 = KC3*(1 + control.tf(1,[TI3,0]) + control.tf([TD3,0],1))
GC_PID3
```




$$\frac{2.03 s^2 + 6.768 s + 5.64}{1.2 s}$$




```python
_, yLCPID3 = control.step_response(control.feedback(GC_PID3*GP,1),ts)
plt.figure(figsize=(14, 7))

plt.plot(ts,yLA,
         ts,yLC1,
         ts,yLCPID1,
         ts,yLCPID2,
         ts,yLCPID3)
plt.legend(('Lazo abierto','Realimentado','PID directo','ZyN abierto','ZyN cerrado'))
plt.grid()
```




- ¿Cuál es la función de ransferencia en lazo cerrado desde la Referencia hacia le señal de control?


$$\frac{Y}{Y_{sp}} = \frac{G_C G_P}{1+G_C G_P}$$

$$E = Y_{sp}-Y$$
$$E + Y = Y_{sp}$$
$$Y = Y_{sp}- E$$
$$Y = Y_{sp}- \frac{U}{G_C}$$

$$U = E G_C$$
$$Y = U G_P$$
$$Y_{sp}- \frac{U}{G_C} = U G_P$$
$$Y_{sp} = U G_P +  \frac{U}{G_C}$$
$$Y_{sp} = U (G_P +  \frac{1}{G_C})$$
$$Y_{sp} = U (\frac{G_C G_P}{G_C} +  \frac{1}{G_C})$$
$$Y_{sp} = U (\frac{G_C G_P + 1}{G_C})$$
$$\frac{Y_{sp}}{U} = (\frac{G_C G_P + 1}{G_C})$$



```python
GC_PID2/(1+GP*GC_PID2)
```




$$\frac{0.001394 s^7 + 3.729 s^6 + 15.54 s^5 + 29.01 s^4 + 31.7 s^3 + 18.25 s^2 + 3.319 s}{6.281 \times 10^{-7} s^7 + 0.003357 s^6 + 4.489 s^5 + 8.596 s^4 + 17.67 s^3 + 15.5 s^2 + 6.637 s}$$



- ¿Cómo se comportan los distintos controladores ante una señal escalón?


```python
# código necesario

_, uLCPID1 = control.step_response(control.feedback(sys1=GC_PID2,
                                                    sys2=GP),ts)
plt.figure(figsize=(14, 7))

plt.plot(ts,uLCPID1)
plt.legend(('PID directo'))
plt.grid()
```

## Sintonización basada en criterios integrales

El comportamiento de los lazos de control puede evaluarse mediante la señal de error $e(t) = ref(t) - Y(t)$. A continuación algunos índices basados en el error.

| Índice                                               |        Definición          |
|------------------------------------------------------|------------------------------|
| IAE:Integral del valor absoluto del error            |   $$\int_0^\infty |e(t)|dt$$ |
| ISE:Integral del cuadrado del error                  |   $$\int_0^\infty e(t)^2dt$$ |
| ITAE:Integral de tiempo por valor absoluto del error |   $$\int_0^\infty t|e(t)|dt$$|

Una técnica de diseño consiste en encontrar los parámetros de un compensador PID que minimicen alguno de los índices mencionados con base en un modelo para el sistema. Puede explorar más [aquí](https://www.researchgate.net/publication/260058488_METODOS_DE_SINTONIZACION_DE_CONTROLADORES_PID_QUE_OPERAN_COMO_SERVOMECANISMOS), [aquí](https://www.researchgate.net/publication/260058590_Optimizacion_del_desempeno_de_los_reguladores_y_servomecanismos_PID) y [aquí](https://www.researchgate.net/publication/255641510_SINTONIZACION_DE_CONTROLADORES_PI_Y_PID_UTILIZANDO_LOS_CRITERIOS_INTEGRALES_IAE_E_ITAE)




- Encuentre y simule un controlador PID que minimice el índice ITAE para el sistema


```python
GP
```




$$\frac{s + 2}{2 s^3 + 3 s^2 + 4 s + 1}$$




```python
# código necesario
```

- ¿Cómo se comportan los distintos controladores ante una señal escalón?


```python
# código necesario
```

0e5082411887c2b556d716cc19814dc79aa07d9f
ipynb
Jupyter Notebook
PID.ipynb
pierrediazp/Control
2a185eff5b5dc84045115009e62296174d072220
MIT
PID.ipynb
pierrediazp/Control
2a185eff5b5dc84045115009e62296174d072220
MIT
PID.ipynb
pierrediazp/Control
2a185eff5b5dc84045115009e62296174d072220
MIT
2021-11-18T13:08:36.000Z
2021-11-18T13:08:36.000Z
Qwen/Qwen-72B
1. YES
2. YES