Tutorial de graficación de plano de fase en SLI

La gráfica del plano de fase o diagrama de fase de un Sistema Lineal e Invariante en el tiempo (SLI) proporciona información útil del comportamiento del sistema.  Con él podemos tener un panorama amplio del comportamiento de la variable y su rapidez, su estabilidad, así como de los puntos de máximo y mínimo del sistema.

Un resumen de las gráficas de los planos de fase en sistemas estables sub-amortiguados, críticamente-amortiguados y sobre-amortiguados con excitación constante (escalón unitario) se pueden ver en las siguientes figuras:

Sistema críticamente-amortiguado

Sistema críticamente-amortiguado

Sistema sub-amortiguado

Sistema sub-amortiguado

Sistema sobre-amortiguado

Sistema sobre-amortiguado

Para más detalles sobre el plano de fase, consulte el libro guía .

Obtener el diagrama de fase de un sistema de manera analítica no siempre una tarea fácil y normalmente conduce a una función implícita en las variables involucradas.  Por esta razón se suele usar programas de software para su graficación.

A continuación vamos a exponer con algunos ejemplos,  cómo obtener el plano de fase, usando los programas Máxima, Python y Matlab.

Ejemplo 1:

Considere el sistema sub-amoritguado: y(t)=e^{-t}\sin(4t) .

Derivando, se obtiene: y'(t)=e^{-t}4\cos(4t)-e^{-t}\sin(4 t)

Para graficar el plano de fase:

Máxima:

La función plot2d() tiene la opción de realizar una gráfica paramétrica de dos funciones definidas.  En este caso, el parámetro es la variable t.  El código se muestra a continuación:

(%i1) y(t):=exp(-t)*sin(4*t);Maxima
      v(t):=diff(y(t),t)$
      plot2d ([parametric, y(t), v(t),
      [t,0, 10*2*%pi/4], [nticks, 20]])$

Cuyo resultado es:

Maxima-planofase1

Como se puede ver en el código, no es necesario tener la expresión para la derivada de y(t), pues el operador diff() realiza esta operación.  El rango de tiempo usado, en este caso, son 10 períodos de la señal y la con la opción nticks se define el número de puntos a representar en el gráfico (opcional).

Python:

El código correspondiente, sin incluir las librerías es:

def y(t): return exp(-t)*sin(4*t)python
def v(t): return exp(-t)*(4*cos(4*t)-sin(4*t))
t=linspace(0,10*2*pi/4,500)
plot(y(t),v(t))

Cuyo resultado es:

Plano de fase ipython

Matlab:

El código es:

>> t=[0:0.01:10*pi/4]; Matlab
>> y=exp(-t).*sin(4*t); v=exp(-t).*(4*cos(4*t)-sin(t));
>> plot(y,v)

Cuyo resultado es:

PlanofaseMatlab

Ejemplo 2:

Considere un SLI (en reposo) regido por el operador L_2(D)=D^2+6D+5, el cual es excitado con una entrada correspondiente al escalón unitario u(t).  Dibuje el diagrama de fase del sistema.

Solución:

El PVI correspondiente al problema planteado es: (D^2+6D+5)y(t)=1 con y(0)=y'(0)=0 para t>0. Cuya solución, mediante Máxima se encuentra con la secuencia de comandos:

eq:'diff(y,t,2)+6*'diff(y,t)+5*y=1$Maxima
ode2(eq,y,t)$
ic2(%,t=0,y=0,'diff(y,t)=0);

Cuyo resultado es:

y=-%e^-t/4+%e^(-5*t)/20+1/5Maxima

Como puede verse, el SLI es de tipo sobre-amortiguado y su respuesta en estado estacionario es 1/5=0.2.

Para graficar el plano de fase, hacemos:

y(t):=-%e^-t/4+%e^(-5*t)/20+1/5;Maxima
plot2d([parametric,y(t),diff(y(t),t)],[t,0,10]);

Cuyo resultado es:

planofase-eje2

Como puede verse en la gráfica, el estado final de y(t) es 0.2=1/5 (eje horizontal) y la rapidez máxima de la variable (y'(t)) es superior a 0.13 (eje vertical) cuando y(t) adquiere un valor aproximado de 0.05.

También podemos usar la función rhs() para graficar el resultado directamente, mediante el siguiente código:

eq:'diff(y,t,2)+6*'diff(y,t)+5*y=1$
ode2(eq,y,t)$
sol:ic2(%,t=0,y=0,'diff(y,t)=0);
y(t):=rhs(sol)$Maxima
plot2d([parametric,y(t),diff(y(t),t)],[t,0,10]);