9 Métodos empíricos y semiempíricos
9.1 Método del ensayo y error
Para realizar la sintonía experimentalmente se sigue un procedimiento similar al siguiente:
Se pone el controlador en manual y se eliminan las acciones integral y derivativa (\(\tau_I \to \infty\) y \(\tau_D = 0\))
Fijar la banda proporcional a un valor elevado, por ejemplo, 200 (ganancia proporcional pequeña)
Poner el controlador en automático
Realizar un pequeño cambio en la consigna o en la carga y seguir la respuesta de la variable controlada hasta que se alcance una respuesta estacionaria. Al ser la ganancia tan baja, la respuesta debe ser lenta
Reducir la banda proporcional a la mitad (doblar la ganancia) y realizar un nuevo cambio en la consigna o en la carga
Repetir el paso 5. hasta que la respuesta sea muy subamortiguada y oscilatoria. Esta es la ganancia última
Fijar la banda proporcional al doble de su último valor (la mitad de la ganancia proporcional)
Iniciar una acción integral reduciendo el valor de \(\tau_I\) a su mitad. Realizar pequeñas perturbaciones o cambios en la consigna y observar el efecto
Encontrar el valor de \(\tau_I\) que hace que la dinámica del sistema sea muy subamortiguada y fijar \(\tau_I\) al doble de ese valor
Aumentar el valor de \(\tau_D\) y realizar cambios en la consigna o en la carga. Encontrar el valor de \(\tau_D\) que proporciona la mayor acción de control sin amplificar el ruido en el proceso de medida de la señal
Reducir el valor de la banda proporcional en pasos de un 10% hasta que se cumplan las especificaciones deseadas
Es importante decir que este método funciona para casi todos los lazos de control pero no funciona con aquellos sistemas que son inestables para ganancias elevadas y para ganancias bajas, es decir, aquellos que solo sean estables para valores de ganancia proporcional intermedios.
9.2 Métodos de criterio único
En este caso se trata de seleccionar uno de los múltiples parámetros que describen el comportamiento de un lazo de control retroalimentado y utilizarlo para el diseño del controlador. Se elegirán los parámetros del controlador en función del rendimiento deseado para ese parámetro.
Es posible tomar como parámetro el error permanente, de manera que el resultado será que se desean eliminar errores que permanezcan durante periodos largos de tiempo. Este es un criterio de estado estacionario.
También se puede seleccionar un criterio dinámico, como puede ser la frecuencia de oscilación de la parte transitoria de la respuesta. También se pueden tomar parámetros como puede ser el sobrepaso, el tiempo necesario para alcanzar el ±5 % de un cierto valor, el tiempo de levantamiento, etc. Es frecuente que la selección de un cierto criterio haga imposible que se cumpla otro. Por ejemplo, si se trata de reducir el sobrepaso reduciendo la ganancia proporcional, se aumenta el tiempo de levantamiento.
Uno de los criterios más frecuentes es seleccionar los parámetros del controlador para que la razón de disminución tome el valor de \(\frac{1}{4}\) (Figura 9.1).
9.3 Método del criterio integral con el tiempo
En este método se toma la respuesta dinámica del sistema de control. En lugar de caracterizarla con un solo parámetro, como en el apartado anterior, se utilizará toda la curva desde \(t = 0\) hasta \(t \to \infty\). La sintonía del controlador se realizará minimizando una de las siguientes integrales con el tiempo:
Integral del error al cuadrado (Integral of the square error):
\[\mathrm{ISE} = \int_0^{\infty} \varepsilon^2 (t) \mathrm{d}t\]
Se minimiza esta integral cuando se desea eliminar errores grandes, ya que estos son los que más contribuyen al valor de la integral.
EjemploEn el problema 9.4 se minimiza la ISE utilizando Julia para sintonizar un lazo de control.
Integral del valor absoluto del error (Integral of the absolute valor of the error):
\[\mathrm{IAE} = \int_0^{\infty} | \varepsilon (t) | \mathrm{d}t\]
En este caso se trata de eliminar errores pequeños. Estos errores podrían ser muy difíciles de eliminar minimizando la ISE, ya que al elevarlos al cuadrado se hace todavía más pequeños respecto a los errores grandes.
Integral del valor absoluto del error ponderado con el tiempo (Integral of the time-weighted absolute error):
\[\mathrm{ITAE} = \int_0^{\infty} t | \varepsilon (t) | \mathrm{d}t\]
Se utiliza cuando se desea eliminar errores muy persistentes en el tiempo, ya que la integral amplifica los errores que permanecen durante tiempos largos, incluso cuando se trata de errores pequeños.
9.4 Método de la curva de respuesta. Método de Cohen y Coon
Una vez seleccionado el controlador que se va a utilizar para el lazo de retroalimentación es necesario dar valores a los parámetros del mismo para completar el diseño, es lo que se conoce como realizar la sintonía del controlador. Para la selección se puede recurrir a métodos semiempíricos como es el método de Cohen y Coon o el método de Ziegler-Nichols, que se trata en el siguiente apartado.
El método de Cohen y Coon también es conocido como el método de la curva de reacción del proceso y se basa en el método de Ziegler-Nichols de lazo abierto. Para ponerlo en práctica hay que abrir el bucle de retroalimentación desconectando el elemento final de control (Figura 9.2).
Una vez abierto se produce un cambio en escalón de altura \(A\) en la variable \(c (t)\), que actúa sobre el elemento final de control. Se registra el valor medido de la variable controlada \(y_m (t)\) respecto al tiempo. De esta manera se obtiene la curva de respuesta del proceso (CRP). La función de transferencia que relaciona la variación de la entrada, \(c(s)\), y la respuesta, \(y_m (s)\), es:
\[G_{\mathrm{CRP}} = \frac{y_m (s)}{c (s)} \approx G_f G_p G_m\]
Esta función de transferencia no solo depende de la dinámica del proceso, también incluye a aquellos elementos que pueda haber entre la entrada y la salida. Típicamente, incluirá al elemento final de control y al medidor, aunque puede haber otros elementos.
Cohen y Coon observaron que la mayoría de procesos presentaban un CRP de aspecto parecido a una sigmoide, que se podía aproximar a un proceso de primer orden con un retraso:
\[G_{\mathrm{CRP}} = \frac{y_m (s)}{c (s)} \approx \frac{K}{\tau s + 1} \mathrm{e}^{- t_d s}\]
donde \(K\), \(\tau\) y \(t_d\) se determinan a partir del análisis de la curva de respuesta del proceso (Figura 9.3):
\[\begin{aligned} K &= \frac{B}{A} = \frac{\text{Salida en estado estacionario}}{\text{Entrada en estado estacionario}} \\ \tau &= \left(1-\frac{1}{\mathrm{e}}\right) B\\ t_d &= \text{Tiempo transcurrido hasta que el sistema responde} \end{aligned}\]
Para realizar la sintonía del controlador se recurre a las siguientes ecuaciones:
Controladores P:
\[K_c = \frac{1}{K} \frac{\tau}{t_d} \left( 1 + \frac{t_d}{3 \tau} \right)\]
Controladores PI:
\[\begin{aligned} K_c &= \frac{1}{K} \frac{\tau}{t_d} \left( 0.9 + \frac{t_d}{12 \tau}\right)\\ \tau_I &= t_d \frac{30 + 3 \frac{t_d}{\tau}}{9 + 20 \frac{t_d}{\tau}} \end{aligned}\]
Controladores PID:
\[\begin{aligned} K_c &= \frac{1}{K} \frac{\tau}{t_d} \left( \frac{4}{3} + \frac{t_d}{4\tau} \right)\\ \tau_I &= t_d \frac{32 + 6 \frac{t_d}{\tau}}{13 + 8 \frac{t_d}{\tau}}\\ \tau_D &= t_d \frac{4}{11 + 2 \frac{t_d}{\tau}} \end{aligned}\]
En Hägglund y Åström (2002) y en Åström y Hägglund (2004) se proponen nuevos valores para la sintonía con el objetivo de que presenten una respuesta menos oscilatoria:
Controladores PI:
\[\begin{aligned} K_c &= \frac{0.14}{K_p} + \frac{0.28 \tau}{t_d}{K_p}\\ \tau_I &= 0.33 t_d + \frac{6.8 \tau t_d}{10 t_d + \tau} \end{aligned}\]
Controladores PID
\[\begin{aligned} K_c &= \frac{1}{K_p}\left(0.2+0.45\frac{\tau}{t_d}\right)\\ \tau_I &= \frac{0.4 t_d + 0.8 \tau}{t_d + 0.1 \tau} t_d\\ \tau_D &= \frac{0.5 \tau t_d}{0.3 t_d + \tau} \end{aligned}\]
9.5 Método de Ziegler y Nichols
Al contrario que el método anterior, esta es una técnica de lazo cerrado. El procedimiento es el siguiente:
Llevar el sistema al nivel de operación deseado, que normalmente corresponderá con las condiciones de diseño
Utilizando solo el control proporcional y con el ciclo cerrado, introducir un cambio en la consigna y variar la ganancia proporcional hasta que el sistema varíe continuamente. La frecuencia de oscilación sostenida es la frecuencia de cruce \(\omega_{co}\). La razón de amplitudes de la respuesta en la frecuencia de cruce es \(RA_{co}\)
Se calculan las siguientes magnitudes:
\[\begin{aligned} \text{Ganancia última: } & K_u = \frac{1}{RA{co}} & \\ \text{Periodo último de oscilaciones sostenidas: } & P_u = \frac{2 \pi}{\omega_{co}} & \end{aligned}\]
Ziegler y Nichols recomendaron los siguientes parámetros para controladores por retroalimentación:
P:
\[K_c = \frac{K_u}{2}\]
PI:
\[\begin{aligned} K_c &= \frac{K_u}{2.2}\\ \tau_I &= \frac{P_u}{1.2} \end{aligned}\]
PID:
\[\begin{aligned} K_c &= \frac{K_u}{1.7}\\ \tau_I &= \frac{P_u}{2}\\ \tau_D &= \frac{P_u}{8} \end{aligned}\]
En Hägglund y J.Åström (2004) se propone otra sintonía para el controlador PI con objeto de obtener una respuesta menos oscilatoria:
\[\begin{aligned} K_c &= K_u 0.15\\ \tau_I &= P_u 0.17 \end{aligned}\]
En O’Dwyer (2009) se puede encontrar una extensísima lista con ecuaciones similares, con reglas para realizar la sintonía de controladores PI y PID en las más diversas situaciones.
9.6 Problemas
Problema 9.1 ★
Sea el sistema de control representado en la figura:
donde \(G_1=\frac{1}{s+1}\) y \(G_2 = \exp(-1.02 s)\).
- Si \(G_c = K_c\), determinar el error permanente de la respuesta para una entrada en escalón unidad.
- Para eliminar el error permanente se recomienda que el controlador sea PID. ¿Qué valores de diseño recomendaría para los parámetros del controlador PID? Se sugiere usar el método de Ziegler-Nichols.
Solución
- En primer lugar, encontraremos la función de transferencia del lazo de control, \(G=\frac{C}{R}\), y, a continuación calcularemos el error permanente.
Empezamos cargando el paquete ClaseControl.jl::
using ClaseControlBuscamos la función de transferencia \(G\) y calculamos el error permanente:
# Definición de la variable Kc
@syms s, t::real
@syms Kc::positive=>"K_c"
# Definición de las funciones de transferencia del lazo
G1 = 1/(s+1)
G2 = exp(-1.02s)
Gc = Kc
# Función de transferencia del lazo de control
G = Gc*G1*G2/(1+Gc*G1*G2)
# Cambio en la consigna en forma de escalón unidad
R = 1/s
# Respuesta del lazo de control
C = G*R
# Cálculo del offset
Offset = sympy.limit(s*(R-G*R), s, 0)\(\frac{1}{K_{c} + 1}\)
- Para el método de Ziegler-Nichols, en primer lugar, hay que encontrar la frecuencia de cruce. Es decir, aquella frecuencia que provoca un retraso de \(-\pi\) rad. A partir de este valor, se puede encontrar la ganancia última \(K_u\), que marca el límite de estabilidad, y el periodo último \(P_u\) en esas condiciones.
Podemos encontrar \(\omega_{co}\) a partir del diagrama de Bode:
salida = ClaseControl.bode(lambdify(G1*G2); co=true)
salida.figsalida.wco1.995576456564642Es importante destacar que en el primero de los gráficos realmente se está representando \(\frac{RA}{K_c}\).
Con el valor de la frecuencia de cruce, se puede calcular el periodo último de oscilación:
\[P_u =\frac{2 \pi}{\omega_{co}}\]
Pu = 2*pi/salida.wco3.1485565419006822La ganancia última es el inverso de la razón de amplitudes para la frecuencia de cruce:
\[K_u = \frac{1}{RA_{co}}\]
Ku = 1/salida.RAco2.2321123166173993Sustituyendo los valores \(K_u\) y \(P_u\) en las fórmulas de Ziegler y Nichols se encuentra la sintonía del controlador PID:
Kc = Ku/1.71.3130072450690584τ_I = Pu/21.5742782709503411τD = Pu/80.3935695677375853Problema 9.2 ★
Determinar los parámetros efectivos de un sistema a partir de la curva de reacción que se indica y calcular la frecuencia crítica y la ganancia máxima. Determinar los ajustes de los parámetros de un controlador PID según el método de Ziegler-Nichols y compararlos con los que se obtienen directamente de la curva de reacción.
| Tiempo (min) | Resp. (ua) |
|---|---|
| 0 | 0 |
| 1 | 0 |
| 2 | 0 |
| 3 | 4 |
| 4 | 10 |
| 5 | 19 |
| 6 | 27 |
| 7 | 35 |
| 8 | 41 |
| 9 | 45 |
| 30 | 50 |
Solución
A partir de los datos de la curva de reacción se pueden obtener los parámetros de la función de transferencia de lazo abierto de varias maneras. En cualquier caso se supone que la curva de reacción del proceso sigue un modelo de primer orden con retraso.
Determinación de los parámetros de la curva de respuesta del proceso
Se representan los puntos disponibles y se encuentran en el gráfico los puntos \(B\), \(\tau\) y \(t_d\), tal como se muestra en la figura a continuación:
t = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 30]
Resp = [0, 0, 0, 4, 10, 19, 27, 35, 41, 45, 50]
scatter(t, Resp, legend=false, xlabel="t (min)",
ylabel="Resp (ua)", dpi=350)A partir del gráfico se observa que \(B\) toma un valor de 50 ua. Suponiendo que se haya realizado una entrada en forma de escalón unidad (\(A=1\)), se puede obtener directamente el valor de la ganancia del proceso (\(K_p=\frac{B}{A}\)):
Kp = 5050Para encontrar el retraso, determinamos cuando empieza a responder el proceso, tal como muestra la Figura 9.3. Para ellos, hemos considerado la recta que pasa entre los puntos (4 min, 10 ua) y (6 min, 27 ua):
# Cálculo de la pendiente (a) y de la ordenada en el origen (b)
a = (27-10)/(6-4)
b = 27 - a*6
# Representación de la recta
plot!(t->a*t+b, 2.0, 7.0)Como consecuencia el valor del retraso será el tiempo que hace que esa recta tenga un valor \(Resp = 0\):
td = -b/a2.823529411764706Para determinar la constante de tiempo del proceso (\(\tau_p\)), deberemos encontrar el tiempo para el que la \(Resp\) es \(B (1-1/\mathrm{e})\)). El tiempo que encontraremos será \(\tau_p + t_d\), por lo que será necesario restar el valor del retraso que hemos encontrado más arriba:
# Utilizaremos un algoritmo de interpolación lineal para encontrar el
# tiempo que hace que Rep = B*(1-exp(-1))
using Interpolations, Roots
itp = LinearInterpolation(t, Resp);
# Resolvemos la ecuación numéricamente
Tp = find_zero(t->itp(t)-50*(1-exp(-1)), 8)-td3.7522240809137792Comprobamos la bondad del ajuste representando la respuesta del proceso de primer orden con un retraso utilizando los valores obtenidos a partir de la curva de respuesta del proceso:
R_t(t) = Kp*(1-exp(-(t-td)/Tp))*(t>td)
scatter(t, Resp, legend=false, xlabel="t (min)", ylabel="Resp (ua)")
plot!(R_t, lw=2)Se comprueba que se obtiene un ajuste razonable considerando la simplicidad del método utilizado. La razón más probable por la que el ajuste no es mejor es porque no se trate realmente de un proceso de primer orden. Probablemente, un proceso de segundo orden sobreamortiguado con un retraso proporcionaría mejor ajuste.
Sintonía con el método de Cohen y Coon
Para encontrar la sintonía por el método de lazo abierto o de la curva de respuesta del proceso, solo hay que sustituir en las fórmulas propias de un controlador PID:
Kc_cc = 1/Kp*Tp/td*(4/3+td/4/Tp)0.040437671875296795Ti_cc = td*(32+6*td/Tp)/(13+8*td/Tp)5.420678851042273Td_cc = td*4/(11+2*td/Tp)0.9031688839998492Sintonía por el método de Ziegler y Nichols
Para aplicar este método deberemos encontrar el periodo último de oscilación y la ganancia última considerando un controlador proporcional.
Obtendremos estos valores considerando la siguiente función de transferencia de lazo abierto:
\[G_{OL} = K_c G_{CRP}\]
Al dibujar el diagrama de Bode, representaremos \(\frac{RA}{K_c}\), lo que permitirá encontrar el parámetro \(M\), a partir del que calcularemos la ganancia última \(K_u\). El valor de \(P_u\) lo calcularemos a partir de la frecuencia de cruce:
salida = ClaseControl.bode(s->Kp/(Tp*s+1)*exp(-td*s); wmax=1, co=true)
salida.figPor lo tanto:
wco = salida.wco0.6873341673074302M = salida.RAco18.075878202106487Lo que supone:
Ku = 1/M0.055322346655525934Pu = 2*pi/wco9.141383632641746La sintonía del controlador será:
Kc_zn = Ku/1.70.03254255685619173Ti_zn = Pu/24.570691816320873Td_zn = Pu/81.1426729540802183Comparativa de las dos sintonías
Estas son las sintonías obtenidas:
| Cohen-Coon | Ziegler-Nichols | |
|---|---|---|
| \(K_c\) | 0.040 | 0.032 |
| \(\tau_I\) | 4.4 | 4.6 |
| \(\tau_d\) | 0.90 | 1.1 |
Como se puede observar la sintonía obtenida por el método de Cohen y Coon es más agresiva, tiene una mayor ganancia proporcional y una menor constante de tiempo integral. Por otro lado, la acción de control derivativa de la sintonía propuesta por el método de Ziegler y Nichols es más elevada.
Para comparar el rendimiento de ambas sintonías se puede realizar una cambio en la consigna en forma de escalón unidad:
Gcrp = Kp/(Tp*s+1)*exp(-td*s)
# Sintonía de Cohen-Coon
Gc_cc = Kc_cc*(1+1/(Ti_cc*s)+Td_cc*s)
G_cc = Gc_cc*Gcrp/(1+Gc_cc*Gcrp)
# Sintonía de Ziegler-Nichols
Gc_zn = Kc_zn*(1+1/(Ti_zn*s)+Td_zn*s)
G_zn = Gc_zn*Gcrp/(1+Gc_zn*Gcrp)\(\frac{50 \cdot \left(0.0371854995761881 s + 0.0325425568561917 + \frac{0.00711983178126118}{s}\right) e^{- 2.82352941176471 s}}{\left(1 + \frac{50 \cdot \left(0.0371854995761881 s + 0.0325425568561917 + \frac{0.00711983178126118}{s}\right) e^{- 2.82352941176471 s}}{3.75222408091378 s + 1}\right) \left(3.75222408091378 s + 1\right)}\)
En este caso la transformada inversa de Laplace no parece tener solución analítica, por lo que la realizaremos numéricamente utilizando el método de Talbot:
using InverseLaplace
# Constructores de las transformadas inversas de Laplace
y_cc = ILT(lambdify(G_cc*1/s))
y_zn = ILT(lambdify(G_zn*1/s))
# Valores de tiempo para los que haremos la simulación
tt = range(0, 25, step=0.1)
# Cálculo de las respuesta de los lazos de control
y_cc_t = map(y_cc, tt)
y_zn_t = map(y_zn, tt)
# Representación de las simulaciones realizadas
plot(tt, y_cc_t, label="Cohen-Coon", xlabel="t", ylabel="R(t)", lw=2, dpi=350)
plot!(tt, y_zn_t, label="Ziegler-Nichols", lw=2)Se puede comprobar que la sintonía de Cohen-Coon tiene una respuesta ligeramente más rápida (tiempo necesario en alcanzar el valor estacionario por primera vez). Pero tiene como inconvenientes un mayor sobrepaso y una respuesta más subamortiguada.
Problema 9.3 ★
Determinar la ganancia de un controlador proporcional para que la razón de disminución de lazo cerrado sea 1/4. La función de transferencia que describe el proceso es:
\[G_p(s)=\frac{1}{s^2+3s+1}\]
Las funciones del medidor y del elemento final de control son iguales a la unidad.
Solución
El lazo de control que propone el problema es:
Por tanto, la función de transferencia que representará la dinámica del bucle de control es:
\[G = \frac{G_c G_p}{1 + G_c G_p}\]
Sustituyendo y operando se encuentra:
\[G = \frac{\frac{\mathrm{Kc}}{K_c + 1}}{\frac{1}{K_c + 1} s^2 + \frac{3}{K_c + 1} s + 1}\]
Al tratarse de un sistema de segundo orden, la constante de tiempo es:
\[\tau = \sqrt{\frac{1}{K_c + 1}}\]
Para un sistema de segundo orden la razón de disminución es:
\[RD = \exp \left( - \frac{2 \pi \zeta}{\sqrt{1 - \zeta^2}} \right)\]
La razón de disminución debe ser de \(\frac{1}{4}\), lo que supone que el coeficiente de amortiguamiento es, tras resolver la ecuación:
\[\zeta = \frac{\sqrt{5727}}{\sqrt{11920 \pi^2 + 5727}}\]
Conocido \(\zeta\) calcular la ganancia del controlador resulta trivial:
\[K_c = \frac{11920 \pi^2}{5727}\]
Problema 9.4 ★
Seleccionar la ganancia y tiempo integral de un controlador PI, empleando el criterio de hacer mínima la ISE. Considerar un cambio en escalón unidad para la consigna. El proceso a controlar es de primer orden con ganancia 10 y constante de tiempo 1.0. Asumir que las funciones de transferencia del medidor y del elemento final de control son iguales a la unidad. Los parámetros seleccionados deben cumplir las siguientes restricciones:
\[\begin{cases} 100 \ge K_c \ge 1\\ 10 > \tau_I > 0.1 \end{cases}\]
Solución
Para minimizar la integral del error al cuadrado simularemos el proceso y evaluaremos numéricamente el valor de la integral. Utilizaremos un algoritmo de optimización para encontrar la sintonía óptima.
En primer lugar, cargaremos las bibliotecas necesarias. Son las siguientes:
SymPy.jl: Para cálculo simbólico, la utilizaremos para realizar los cálculos previos con las funciones de transferencia.Plots.jl: Para la representación de gráficos.InverseLaplace.jl: Para evaluar numéricamente las transformadas inversas de Laplace, la utilizaremos para realizar las simulaciones del lazo de control.Trapz.jl: Para calcular numéricamente la integral utilizando el método de los trapecios.Optim.jl: Para realizar el proceso de optimización.
using SymPy, InverseLaplace, Plots, Optim, TrapzEmpezaremos definiendo las variables necesarias para poder realizas los cálculos con Sympy:
@syms Kc::positive=>"K_c" Ti::positive=>"tau_I" s t::real(K_c, tau_I, s, t)A continuación definimos las funciones de transferencia del proceso y del controlador PI:
Gc = Kc*(1+1/Ti/s)
Gp = 10/(s+1)\(\frac{10}{s + 1}\)
En este caso nos interesa conocer como cambia el error \(\varepsilon(t)\) en función de un cambio en la consigna \(y_{sp}\). Podemos redibujar el lazo de control de la manera siguiente para poner en evidencia que la consigna será la entrada y el error la salida del bloque equivalente:
Se trata del mismo lazo de control, lo único que hemos hecho ha sido redibujarlo para indicar con mayor claridad la entrada y la salida que formarán el bloque equivalente:
La función de transferencia del bloque equivalente será:
G = 1/(1+Gc*Gp)\(\frac{1}{\frac{10 K_{c} \left(1 + \frac{1}{s \tau_{I}}\right)}{s + 1} + 1}\)
La transformada de Laplace del error para un cambio en la consigna en forma de escalón unidad es:
e_s = G*1/s\(\frac{1}{s \left(\frac{10 K_{c} \left(1 + \frac{1}{s \tau_{I}}\right)}{s + 1} + 1\right)}\)
El siguiente paso consiste en definir la función que calculará la ISE:
function ise(x)
# Los parámetros del controlador se expresan como el vector x, que
# tiene las siguientes compoenentes (Kc, Ti)
# Tiempo de la simulación
tsim = range(0, 0.5; length=100)
# Cálculo de error(t) al cuadrado
# 1. Definición de la transformada inversa de Laplace numérica
error = ILT(lambdify(e_s(Kc=>x[1], Ti=>x[2])))
# 2. Obtención del error y se eleva al cuadrado
e2sim = map(error, tsim).^2
# Convertimos los valores NaN en 0.0, para evitar problemas
# con trapz()
for i in range(1, length(e2sim), step=1)
if isnan(e2sim[i])
e2sim[i] = 0.0
end
end
# Cálculo de la integral
ise = trapz(tsim, e2sim)
endise (generic function with 1 method)Finalmente, solo queda plantear la optimización, para la que se minimizará el valor de ISE considerando los límites de los parámetros del controlador que propone el enunciado del problema. Como valores iniciales se toma \(K_{c_0}=2\) y \(\tau_{I_0}=9\):
result = optimize(ise, [1.0, 0.1], [100.0, 10.0], [2.0, 9.0]) * Status: success
* Candidate solution
Final objective value: 2.040391e-07
* Found with
Algorithm: Fminbox with L-BFGS
* Convergence measures
|x - x'| = 4.07e-07 ≰ 0.0e+00
|x - x'|/|x'| = 4.07e-09 ≰ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 4.47e-10 ≤ 1.0e-08
* Work counters
Seconds run: 51 (vs limit Inf)
Iterations: 5
f(x) calls: 162
∇f(x) calls: 162Creamos dos nuevas variables con los valores de los parámetros del controlador que minimizan la ISE:
Kc_optim, Ti_optim = Optim.minimizer(result)2-element Vector{Float64}:
99.99999999955283
0.9046308215714349Finalmente, comprobamos la sintonía del controlador simulando la respuesta del lazo de control para un cambio de la consigna en forma de escalón unidad:
y_s = Gc*Gp/(1+Gc*Gp)*1/s
y_t = ILT(lambdify(y_s(Kc=>Kc_optim, Ti=>Ti_optim)))
tsim = range(0, .025; length=100)
ysim = map(y_t, tsim)
plot(tsim, ysim, legend=false, xlabel="t", ylabel="y(t)")Problema 9.5 ★
Seleccionar la ganancia de un controlador proporcional utilizando el criterio de la razón de disminución 1/4. El proceso a controlar es:
\[G_p(s) = \frac{10}{(s+2)(2s+1)}\]
Asumir que \(G_m(s) = G_f(s) = 1\). Realizar también la sintonía utilizando el criterio del ISE mínimo con un cambio en la consigna en escalón unidad. En ambos casos se debe satisfacer la condición de que \(100 \ge K_c \ge 0.1\). Comparar las sintonías y explicar las diferencias entre ellas.
Solución
Razón de disminución 1/4
El bucle de control propuesto por el enunciado del problema es:
La función de transferencia que representa la dinámica de este lazo es:
\[G = \frac{K_c G_p}{1 + K_c G_p} = \frac{\frac{10 \mathrm{Kc}}{10 \mathrm{Kc} + 2}}{\frac{2}{10 \mathrm{Kc} + 2} s^2 + \frac{5}{10 \mathrm{Kc} + 2} s + 1}\]
Se trata de un sistema de segundo orden, por tanto:
\[\begin{aligned} & \tau^2 = \frac{2}{10 \mathrm{Kc} + 2} & \\ & 2 \tau \zeta = \frac{5}{10 \mathrm{Kc} + 2} & \end{aligned}\]
Como criterio de sintonía se toma que la razón de disminución de este lazo de control debe ser de 1/4:
\[RD = \exp \left( \frac{- 2 \pi \zeta}{\sqrt{1 - \zeta^2}} \right) = \frac{1}{4}\]
Resolviendo la ecuación anterior se obtiene:
\[\zeta = 0.2154\]
Conocido el coeficiente de amortiguamiento, ya solo queda resolver el sistema de ecuaciones anterior para obtener el valor de la ganancia proporcional:
\[K_c = 6.53\]
Minimizar ISE
Para obtener la sintonía óptima, aquella que minimiza la integral del error al cuadrado, del controlador se ha utilizado el programa siguiente:
using InverseLaplace, Optim, Trapz
# Definición de la función de transferencia del proceso
G_p(s) = 10/((s+2)*(2s+1))
# Definición de la función coste: ISE
function ISE(Kc)
# Función de transferencia del lazo de control
G(s) = Kc*G_p(s)/(1+Kc*G_p(s))
# Respuesta del lazo de control para un cambio en la consigna
# en forma de escalón unidad
y(s) = G(s)*1/s
# Tiempo de la simulación
tsim = range(0.0, 2.0; length=200)
# Transformada inversa de Laplace numérica para obtener y(t)
y_t = ILT(y, 80)
# Cálculo de los valores simulados de la variable controlada
ysim = map(y_t, tsim)
# Eliminación de los NaN
for i in 1:length(ysim)
if isnan(ysim[i])
ysim[i] = 0.0
end
end
# Cálculo de la integral del error al cuadrado
ISE = trapz(tsim, (1.0 .-ysim).^2)
end
# Optimización univariante
result= optimize(ISE, 0.1, 100.0)Results of Optimization Algorithm
* Status: success
* Algorithm: Brent's Method
* Search Interval: [0.100000, 100.000000]
* Minimizer: 9.415487e+01
* Minimum: 2.003754e-01
* Iterations: 32
* Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
* Objective Function Calls: 33La ganancia proporcional obtenida es:
Kc_optim = Optim.minimizer(result)94.15486825948487El valor obtenido de la integral del error al cuadrado es:
Optim.minimum(result)0.20037541519520863Podemos compararlo con el valor obtenido con la ganancia proporcional obtenida para una razón de disminución de 1/4:
ISE(6.53)0.22714470820459098Se observa una ligera reducción en el ISE.
Finalmente podemos comparar las respuestas para las dos sintonías obtenidas:
using Plots
function ysimul(tsim, Kc)
G(s) = Kc*G_p(s)/(1+Kc*G_p(s))
y(s) = G(s)*1/s
y_t = ILT(y, 80)
ysim = map(y_t, tsim)
end
tsim = range(0.0, 2.0; length=200)
ysim_optim = ysimul(tsim, Kc_optim)
ysim_rd = ysimul(tsim, 6.53)
plot(tsim, ysim_optim, lw=2, label="min(ISE)", xlabel="t", ylabel="y(t)")
plot!(tsim, ysim_rd, lw=2, label="RD=1/4")
hline!([1.0], label="Consigna")Problema 9.6
Repetir el problema anterior utilizando la técnica de Ziegler-Nichols en lugar de minimizar el ISE. Comparar los resultados obtenidos con los del problema anterior.
Problema 9.7
Considerar un lazo de control con las siguientes funciones de transferencia:
\[\begin{aligned} G_f &= 5\\ G_p &= \frac{10}{s+4}\\ G_m &= \frac{1}{10s+1} \end{aligned}\]
Realizar la sintonía de un controlador PI utilizando la técnica de Cohen-Coon.
Dibujar la curva real de reacción del proceso junto con una aproximación de primer orden con retraso.
Problema 9.8
Repetir el problema anterior con un controlador PID y las siguientes funciones de transferencia:
\[\begin{aligned} G_f &= 1\\ G_p &= \frac{\mathrm{e}^{-s}}{s+10}\\ G_m &= \frac{5 \mathrm{e}^{-0.1 s}}{0.01s+1} \end{aligned}\]
Problema 9.9 ★
Se obtiene experimentalmente la curva de reacción de un proceso y se obtienen los siguientes datos:
| Tiempo (min) | Variable manipulable | Salida medida |
|---|---|---|
| -2 | 100 | 200 |
| -1 | 100 | 200 |
| 0 | 150 | 200.1 |
| 0.2 | 150 | 201.1 |
| 0.4 | 150 | 204.0 |
| 0.6 | 150 | 227.0 |
| 0.8 | 150 | 251.0 |
| 1.0 | 150 | 280.0 |
| 1.2 | 150 | 302.5 |
| 1.4 | 150 | 318.0 |
| 1.6 | 150 | 329.5 |
| 1.8 | 150 | 336.0 |
| 2.0 | 150 | 339.0 |
| 2.2 | 150 | 340.5 |
| 2.4 | 150 | 341.0 |
Usando esos valores:
Aproximar la respuesta de lazo abierto a un sistema de primer orden con retraso.
Seleccionar los parámetros de un controlador PI utilizando la técnica de Cohen-Coon.
Solución
Una de las primeras tareas cuando se dispone de datos experimentales es representarlos gráficamente:
using Plots, LaTeXStrings
t = [-2.0, -1.0, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4]
c = [100, 100, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150]
ym = [200, 200, 200.1, 201.1, 204.0, 227.0, 251.0, 280.0, 302.5, 318.0, 329.5, 336.0, 339.0, 340.5, 341.0]
scatter(t, c, label=L"c(t)", legend=:topleft, xlabel=L"t")
scatter!(t, ym, label=L"y_m(t)")Observamos que el cambio de la variable manipulable (\(c(t)\)) tiene forma de escalón y que la respuesta del proceso (\(y_m(t)\)) tiene un aspecto que parece ser compatible con la de un sistema de primer orden.
Para poder aplicar el método de la curva de reacción del proceso para poder identificar los parámetros del proceso equivalente de primer orden (\(K\), \(\tau\) y \(t_d\)) deberemos, en primer lugar, definir las variables de desviación. A partir del análisis del gráfico o de la tabla de datos que los valores en estado estacionario son:
\[\begin{aligned} c_e = 100\\ y_{m,e} = 200 \end{aligned}\]
Por lo tanto las variables de desviación serán:
\[\begin{aligned} C(t) = c(t) - c_e = c(t) - 100\\ Y_m = y_m(t) - y_{m,e} = y_m(t) - 200 \end{aligned}\]
Además, solo consideraremos los puntos del 3 hasta el último, es decir, a partir del momento que se introduce el escalón:
tdata = t[3:end]
C = c[3:end] .- c[1]
Ym = ym[3:end] .- ym[1]13-element Vector{Float64}:
0.09999999999999432
1.0999999999999943
4.0
27.0
51.0
80.0
102.5
118.0
129.5
136.0
139.0
140.5
141.0La representación gráfica de las variables de desviación queda como:
scatter(tdata, C, label=L"C(t)", legend=:topleft, xlabel=L"t")
scatter!(tdata, Ym, label=L"Y_m(t)")En el gráfico se observa que la altura del escalón es:
A = C[1]50y como la respuesta en estado estacionario de la curva de reacción del proceso es:
B = Ym[end]141.0se puede encontrar fácilmente la ganancia de la curva de reacción del proceso:
K = B/A2.82Estimamos la posición del retraso en:
td = 0.450.45Encontramos la el valor del retraso mediante interpolación lineal para el punto en el que \(Y_m = B (1-1/\mathrm{e})\):
using Interpolations, Roots
itp = LinearInterpolation(tdata, Ym)
f(x)=itp(x)-B*(1-1/exp(1))
T = find_zero(f, 1)-td0.6311466559540144De manera que se obtienen los siguientes parámetros:
# Ganancia de Gcrp
K2.82# Constante de tiempo de Gcrp
round(T, digits=3)0.631# Retraso de Gcrp
td0.45Se puede comprobar la bondad del ajuste representando la curva de respuesta del proceso junto con los datos experimentales:
scatter(tdata, Ym, legend=:bottomright, label="Exp")
ysim_2(t, K, T, td) = K*A*(1-exp(-(t-td)/T))*(t-td>0)
tsim = tdata[1]:0.01:tdata[end]
plot!(tsim, ysim_2.(tsim, K, T, td), label="CRP", lw=2)Se observa que el ajuste empeora a medida de que pasa el tiempo. Con toda probabilidad es debido a que el experimento ha finalizado demasiado pronto.
También podríamos haber obtenido los parámetros de \(G_{CRP}\) mediante regresión no lineal utilizando la biblioteca LsqFit:
# using Pkg; Pkg.add("LsqFit")
using LsqFitDefinimos el modelo para el ajuste, que se corresponde con la respuesta de un sistema de primer orden con retraso para una entrada en escalón de magnitud \(A\):
# p[1] es la ganancia
# p[2] es la constante de tiempo
# p[3] es el retraso
@. model(t, p) = p[1]*A*(1-exp(-(t-p[3])/p[2]))model (generic function with 1 method)Realizamos el ajuste no lineal:
fit = curve_fit(model, tdata, Ym, [1.0,1.0,1.0]; autodiff=:forwarddiff)LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}, Vector{LsqFit.LMState{LsqFit.LevenbergMarquardt}}}([4.802447943261538, 2.146641188060415, 0.1704415269233574], [-19.942848286106212, 2.183738765096502, 20.352977620725124, 16.547890600002994, 10.035218301146834, -3.0331319717963225, -11.018766634602585, -13.295603638258598, -12.748777231056735, -8.273640153386225, -1.27485459757807, 6.334435963762047, 14.13336127501833], [-4.131819547143312 9.615498215011632 -121.103259778628; 0.6837635314098534 -1.519202519332268 -110.32987707273732; … ; 30.57491464739143 -41.08736156945702 -43.45764057737317; 32.30297612964045 -41.12092291916432 -39.591635696159315], true, Iter Function value Gradient norm
------ -------------- --------------
, Float64[])Para facilitar el análisis, vamos a crear variables con los resultados del ajuste. Usaremos el subíndice lsq para indicar que son los valores obtenidos por regresión no lineal:
K_lsq = fit.param[1]
round(K_lsq, digits=3)4.802T_lsq = fit.param[2]
round(T_lsq, digits=3)2.147td_lsq = fit.param[3]
round(td_lsq, digits=3)0.17Al representar los dos modelos, podemos comprobar que el ajuste mediante regresión no lineal se ve todavía más afectado por la falta de datos a tiempos largos.
scatter(tdata, Ym, legend=:bottomright, label="Exp")
ysim_3(t, K, T, td) = K*A*(1-exp(-(t-td)/T))*(t-td>0)
plot!(tsim, ysim_3.(tsim, K, T, td), label="CRP", lw=2)
plot!(tsim, ysim_3.(tsim, K_lsq, T_lsq, td_lsq), label="LSQ", lw=2)Quizás de manera sorprendente hemos encontrado que el ajuste utilizado por una técnica que podríamos considerar más rudimentaria nos proporciona un modelo experimental que tiene un comportamiento mejor que el proporcionado por un método más sofisticado y complejo, como es la regresión no lineal. La mayor sofisticación de la técnica no implica mejores resultados.
Calcular los parámetros del controlador PI resulta una tarea trivial:
Kc = 1/K*T/td*(0.9+td/12/T)
round(Kc, digits=3)0.477Ti = td*(30+3*td/T)/(9+20*td/T)
round(Ti, digits=3)0.622Problema 9.10
Utilizando los datos del problema anterior contestar a las siguientes preguntas:
Sintonizar un controlador PI utilizando la técnica de Ziegler-Nichols.
Comparar los resultados de la sintonía obtenidos con la técnica de Cohen-Coon y a los obtenidos al minimizar la ISE para un cambio en la consigna en escalón unidad.
Comparar la tolerancia de las diferentes sintonías a errores en la ganancia, constante de tiempo o tiempo muerto. ¿Cuál de ellas posee una tolerancia mayor?
Problema 9.11 ★
La curva de reacción de un proceso de un sistema de control de temperatura proporciona los siguientes valores: \(K = 10\), \(\tau = 2\) min y \(t_d = 0.1\) min. Responder:
Realizar la sintonía mediante la técnica de Ziegler-Nichols.
Comparar la sintonía anterior con la obtenida con la técnica de Cohen-Coon.
Asumir que los valores obtenidos con el método de la curva de reacción del proceso no son muy fiables. Calcular qué porcentajes de error de los valores \(K\), \(\tau\) y \(t_d\) puede tolerar la sintonía de Ziegler-Nichols sin volverse inestable.
Solución
- Para obtener la curva de reacción de un proceso, hay que abrir el lazo de control (desconectando el controlador). A continuación hay que introducir un cambio en escalón en el elemento final de control. Para ello se puede utilizar un generador de funciones que simule al controlador, de manera que la variable \(c\) tenga una forma de escalón. Si se trata de una válvula puede ser algo tan sencillo como abrir o cerrar la válvula con rapidez. Evidentemente, el cambio generado debe ser conocido, por ejemplo, en el caso de la válvula se debe conocer el recorrido del émbolo. Como consecuencia en el registrador a la salida del proceso se obtiene la curva de respuesta del proceso, tal como muestra la figura siguiente:
A partir del análisis de la curva obtenida, siguiendo el método propuesto por Cohen y Coon, se obtiene la función de transferencia de la curva de respuesta del proceso, que engloba la dinámica del elemento final de control, proceso y medidor:
Los parámetros indicados en el enunciado del problema se han obtenido de esta manera.
Por tanto, el lazo de control se puede representar según el siguiente bucle:
Para realizar la sintonía del controlador utilizando la técnica de Zielger-Nichols hay que sustituir el controlador existente por un controlador proporcional. A continuación hay que buscar la frecuencia de cruce y la ganancia última del controlador. En este caso la función de transferencia de lazo abierto del bucle de control es:
\[G_{OL} (s) = K_c \frac{10}{2 s + 1} \mathrm{e}^{- 0.1 s}\]
El diagrama de Bode de este lazo de control es:
sol = ClaseControl.bode(s->10/(2s+1)*exp(-0.1s); wmin=0.1, wmax=100, co=true,
RAlabel="RA/Kc");
sol.figround(sol.wco, sigdigits=4)16.02round(sol.RAco, sigdigits=4)0.312La frecuencia de cruce tiene un valor aproximado de 16 rad/min y que \(M\) = 0.31, lo que implica que la ganancia última tiene un valor aproximado de 3.2.
Mediante las ecuaciones de la razón de amplitudes y el desfase de la función de transferencia de lazo abierto también se puede resolver el problema:
\[\begin{aligned} RA_{OL} &= K_c \frac{10}{\sqrt{1 + 4 \omega^2}} \\ \varphi_{OL} &= \mathrm{atan} (- 2 \omega) - 0.1 \omega \end{aligned}\]
En primer lugar hay que encontrar la frecuencia de cruce resolviendo la siguiente ecuación:
\[- \pi = \mathrm{atan} (- 2 \omega_{co}) - 0.1 \omega_{co}\]
Se obtiene una frecuencia de cruce:
\[\omega_{co} = 16.02 \text{ rad/min}\]
Utilizando la frecuencia de cruce se calcula el periodo último (de oscilaciones sostenidas):
\[P_u = \frac{2 \pi}{\omega_{co}} = 0.392 \text{ min}\]
La ganancia última, aquella que marca el límite de estabilidad (\(RA = 1\)), se obtendrá resolviendo la ecuación:
\[1 = K_u \frac{10}{\sqrt{1 + 4 \omega_{co}^2}}\]
La ganancia última es:
\[K_u = 3.206\]
Al no indicar el problema el tipo de controlador a sintonizar, se supondrá que se trata de un PID. Según los valores recomendados por Ziegler y Nichols la sintonía del controlador será:
\[\begin{aligned} K_c &= \frac{K_u}{1.7} = 1.89\\ \tau_I &= \frac{P_u}{2} = 0.196 \text{ min}\\ \tau_D &= \frac{P_u}{8} = 0.049 \text{ min} \end{aligned}\]
- El cálculo de la sintonía utilizando el método de Cohen y Coon es directo ya que la función de transferencia de la curva de reacción del proceso es conocida. La sintonía del controlador PID será:
\[\begin{aligned} K_c &= \frac{1}{K} \frac{\tau}{t_d} \left( \frac{4}{3} + \frac{t_d}{4 \tau} \right) = 2.69\\ \tau_I &= t_d \frac{32 + 6 \frac{t_d}{\tau}}{13 + 8 \frac{t_d}{\tau}} = 0.241 \text{ min}\\ \tau_D &= t_d \frac{4}{11 + 2 \frac{t_d}{\tau}} = 0.036 \text{ min} \end{aligned}\]
Se puede observar que la sintonía propuesta por el método de Ziegler-Nichols tiene una ganancia proporcional menor. Como compensación la acción integral es más intensa, lo que supondrá una respuesta menos amortiguada (más oscilaciones).
c) En el caso de que los parámetros de la curva de respuesta del proceso no sean muy fiables es necesario realizar un estudio de sensibilidad.
La función de transferencia de lazo abierto en este caso es:
\[G_{OL} = 1.89 \left( 1 + \frac{1}{0.196 s} + 0.049 s \right) \frac{K}{\tau s + 1} \mathrm{e}^{- t_d s}\]
Por tanto, la razón de amplitudes y el desfase serán:
\[\begin{aligned} RA_{OL} &= 1.89 \sqrt{\left( 0.049 \omega - \frac{1}{0.196 \omega} \right)^2 + 1} \frac{K}{\sqrt{1 + \tau^2 \omega^2}}\\ \varphi_{OL} &= \mathrm{atan} \left( 0.049 \omega - \frac{1}{0.196 \omega} \right) + \mathrm{atan} (- \tau \omega) - t_d \omega \end{aligned}\]
El lazo de control será estable si \(RA_{OL} \le 1\) para la frecuencia de cruce (\(\varphi_{OL} (\omega_{co}) = -\pi\)). Como simplificación se va a considerar la influencia de cada uno de los parámetros del proceso por separado.
Ganancia del proceso
En este caso:
\[\begin{aligned} RA_{OL} &= 1.89 \sqrt{\left( 0.049 \omega - \frac{1}{0.196 \omega} \right)^2 + 1} \frac{K}{\sqrt{1 + 4 \omega^2}}\\ \varphi_{OL} &= \mathrm{atan} \left( 0.049 \omega - \frac{1}{0.196 \omega} \right) + \mathrm{atan} (- 2 \omega) - 0.1 \omega \end{aligned}\]
El diagrama de Bode de este sistema:
Kc = 1.89
Ti = 0.196
Td = 0.049
sol_K = ClaseControl.bode(s->Kc*(1+1/(Ti*s)+Td*s)*1/(2s+1)*exp(-0.1s);
wmin=10, wmax=100, co=true);
sol_K.figLa frecuencia de cruce, \(\omega_{co}\), tiene un valor de 23.4 rad. La ganancia del proceso que será el límite de estabilidad (\(RA_{OL} (K) = 1\)) se obtendrá resolviendo la siguiente ecuación:
\[1 = 1.89 \sqrt{\left( 0.049 \omega_{co} - \frac{1}{0.196 \omega_{co}} \right)^2 + 1} \frac{K}{\sqrt{1 + 4 \omega_{co}^2}}\]
Se obtiene un valor de \(K\) de 18.1 (\(K = \frac{1}{R{co}}\)), lo que supone que la ganancia del proceso puede aumentar un 81.5%.
Constante de tiempo del proceso
En este caso la razón de amplitudes y el desfase del lazo abierto serán:
\[\begin{aligned} RA_{OL} &= 1.89 \sqrt{\left( 0.049 \omega - \frac{1}{0.196 \omega} \right)^2 + 1} \frac{10}{\sqrt{1 + \tau^2 \omega^2}}\\ \varphi_{OL} &= \mathrm{atan} \left( 0.049 \omega - \frac{1}{0.196 \omega} \right) + \mathrm{atan} (- \tau \omega) - 0.1 \omega \end{aligned}\]
El límete de estabilidad vendrá determinado por las siguientes ecuaciones:
\[\begin{aligned} 1 &= 1.89 \sqrt{\left( 0.049 \omega_{co} - \frac{1}{0.196 \omega_{co}} \right)^2 + 1} \frac{10}{\sqrt{1 + \tau^2_u \omega_{co}^2}}\\ -\pi &= \mathrm{atan} \left( 0.049 \omega_{co} - \frac{1}{0.196 \omega_{co}} \right) + \mathrm{atan} (- \tau_u \omega_{co}) - 0.1 \omega_{co} \end{aligned}\]
donde \(\omega_{co}\) es la frecuencia de cruce y \(\tau_u\) es la constante de tiempo límite de estabilidad, ya que se calcula para una razón de amplitudes igual a la unidad para la frecuencia de cruce.
Se pueden resolver estas ecuaciones numéricamente utilizando la biblioteca NLsolve.jl:
# import Pkg; Pkg.add("NLsolve")
using NLsolve
function f!(F, x)
wco = x[1]
Tu = x[2]
F[1] = Kc*sqrt((Td*wco-1/Ti/wco)^2+1)*10/sqrt(1+Tu^2*wco^2)-1
F[2] = atan(Td*wco-1/Ti/wco)+atan(-Tu*wco)-0.1*wco+pi
end
nlsolve(f!, [20.0, 1.0], autodiff=:forward)Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [20.0, 1.0]
* Zero: [23.6561256979143, 1.0975994533724522]
* Inf-norm of residuals: 0.000000
* Iterations: 4
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 5
* Jacobian Calls (df/dx): 5Se obtiene un valor de la frecuencia de cruce de 23.66 rad y una ganancia del proceso última de 1.10.
Una alternativa consiste en buscar la constante de tiempo última del proceso por tanteo, ya sea utilizando los diagramas de Bode o directamente las ecuaciones. Por ejemplo, para una constante de tiempo de 1.1 min se obtiene el siguiente diagrama de Bode:
sol_T = ClaseControl.bode(s->Kc*(1+1/(Ti*s)+Td*s)*10/(1.1s+1)*exp(-0.1s);
wmin=10, wmax=100, co=true);
sol_T.figSe puede apreciar que el lazo de control está al límite de estabilidad. Por tanto se puede disminuir la constante de tiempo del proceso hasta un 45 % (aproximadamente) manteniendo el sistema estable.
Retraso
La situación es análoga a las anteriores:
\[\begin{aligned} RA_{OL} &= 1.89 \sqrt{\left( 0.049 \omega - \frac{1}{0.196 \omega} \right)^2 + 1} \frac{10}{\sqrt{1 + 4 \omega^2}}\\ \varphi_{OL} &= \mathrm{atan} \left( 0.049 \omega - \frac{1}{0.196 \omega} \right) + \mathrm{atan} (- 2 \omega) - t_d \omega \end{aligned}\]
El valor límite del retraso \(t_{du}\) se puede calcular utilizando las siguientes ecuaciones:
\[\begin{aligned} 1 &= 1.89 \sqrt{\left( 0.049 \omega_{co} - \frac{1}{0.196 \omega_{co}} \right)^2 + 1} \frac{10}{\sqrt{1 + 4 \omega_{co}^2}}\\ -\pi &= \mathrm{atan} \left( 0.049 \omega_{co} - \frac{1}{0.196 \omega_{co}} \right) + \mathrm{atan} (- 2 \omega_{co}) - t_{du} \omega_{co} \end{aligned}\]
Resolviendo la primera de las ecuaciones se encuentra la frecuencia de cruce:
\[\left\{\begin{array}{l} \omega_{co} = \pm \frac{500 \sqrt{20 \sqrt{727254571} + 248647}}{\sqrt{2199637153}}\\ \omega_{co} = \pm \frac{500 \sqrt{20 \sqrt{727254571} - 248647} i}{\sqrt{2199637153}} \end{array}\right.\]
La frecuencia de cruce es un número real y positivo, lo que elimina tres de las soluciones:
\[\omega_{co} = \frac{500 \sqrt{20 \sqrt{727254571} + 248647}}{\sqrt{2199637153}} = 9.464 \text{ rad/min}\]
Sustituyendo en la segunda de las ecuaciones, se obtiene:
\[t_{du} = 0.164 \text{ min}\]
El retraso puede aumentar hasta un 64 % y el sistema continuará siendo estable.