Predictor de Smith Filtrado

En esta entrada veremos una modificación de la estructura del Predictor de Smith que vimos en la entrada pasada (Click aqui para ver la clase del Predictor de Smith), dicha modificación se conoce como el Predictor de Smith Filtrado (PSF) que tiene como principal caracteristica corregir las principales limitaciones del Predictor de Smith Original. Las cuales son acelerar la dinámica de rechazo de perturbación y controlar procesos integradores y inestables.

Preditor de Smith Filtrado

La estructura del predictor de Smith filtrado es representada asi:

Sabemos que las relaciones de entrada y salida vienen dadas por las siguientes ecuaciones, las cuales son prácticamente las mismas vistas en la entrada del Predictor de Smith:

Vemos que la dinámica del Filtro Fr(s) solo aparece en las respuestas de las perturbaciones q(t) y n(t).

Si tomamos la respuesta frente a la perturbación de carga (q(t)), notaremos la mejora que el filtro de predicción Fr(s) introduce en la estructura del PSF:

Predictor de Smith Filtrado en un Invernadero

El siguiente trabajo muestra la aplicación del predictor de Smith filtrado multivariable en un invernadero.

S. A. C. Giraldo, R. C. C. Flesch, and J. E. Normey-Rico, “Multivariable greenhouse control using the filtered  Smith predictor,” Journal of Control, Automation and Electrical Systems, pp. 1–10, 2016. (click aca)

Ejemplo PSF de primer orden

Si tenemos el siguiente sistema:

y Tenemos un controlador primario PI

Si diseñamos un control rápido PI por cancelación de polos osea que obtendremos la siguiente respuesta en lazo cerrado:

donde es la constante de tiempo en lazo cerrado dado por: .

dividimos el filtro en y resuelvo

El objetivo es conseguir que el polinomio cancele los polos lentos de la planta , para eso aplico la propiedad de polinomios que dice que dado un polinomio y otro polinomio notamos que para cancelar el polinomio el polinomio B debe tener la misma raiz que el polinomio A , osea toda raiz de A es raiz de B.

Para cumplir la condición 1 debo imponer que:

Para cumplir la condición 2 debo imponer que:

Como voy a matar una dinámica de primer orden, diseño un filtro de primer orden:

este numerador es usado para cumplir la condición 1.

el denominador del filtro es usado para determinar el tiempo de respuesta de rechazo de perturbación.

Reemplazando y en la ecuación de la condición 1 se tiene que:

Así, de esta ecuación la única incógnita es , entonces solo bastaría con despejarla y hallar su valor.

Determinando el valor de note que la ecuación restante daría algo parecido con:

Daría un polinomio cualquiera, que dentro de su dinámica va a tener el polo lento de la planta y una dinámica , de esa forma note como quedaría la dinámica de rechazo de perturbación:

note como se cancelan los polos lentos del modelo:

Puede notarse que ahora la dinámica de rechazo de perturbación está dado por el dominador del filtro y por los polos de lazo cerrado.

Segunda opción en Calculo del Filtro

Una segunda opción seria diseñar un filtro de segundo orden que además de cancelar los polos lentos, también cancele la dinámica de lazo cerrado. De esa forma se escogería diseñar un filtro con la siguiente forma:

Diseñando este filtro mi respuesta de rechazo de perturbación va a tener una velocidad diferente a la velocidad del seguimiento de referencia.

Aqui se aplica el mismo concepto, debo hacer con que el numerador cancele los polos lentos de la planta , para eso evalúo el denominador con la raíz del polinomio (polo) y encuentro el valor de .

así, se tiene que:

Sistemas Integradores o Inestables

Para un caso integrador o mismo para un caso inestable, para efectos de analisis se puede utilizar la misma estructura del PSF, pero para efectos de implementación práctica, no se puede utilizar la misma estructura del PSF. Debido a que esa estructura es internamente inestable. Ya que si observamos la dinámica de rechazo de perturbaciones, veremos que si la planta () tiene polos inestables, la estructura será internamente inestable:

Para solucionar este problema, basta con tener una estructura equivalente al PSF, la cual la conoceremos como estructura estable:

Donde

Para que la nueva estructura sea internamente estable, hay que asegurarse de que tanto el filtro y el filtro sean estables. En el caso del filtro de es fácil de asegurar su estabilidad, dado que sus polos son elegidos libremente por el diseñador. Por otra parte, el filtro , se puede observar que si tiene polos integradores o polos inestable, el único grado de libertad que se tiene para cancelar esos polos es por medio del filtro , evitando así la inestabilidad interna. Por esta razón, el termino siempre debe ser calculado con el fin de cancelar las dinámicas inestables o integradoras del modelo rápido de la planta, , donde se debe adicionar una raíz en (z-1) para que de esta forma el filtro  tenga ganancia unitaria, y el filtro tiene que ser implementado con las cancelaciones ya realizados. Por lo tanto, la salida predita será una predicción estable de la salida real . (Todo esto es explicado en detalle en el video del ejemplo 2, a continuación.)

Ejemplo 1 (Caso Estable Continuo)

Comencemos con un proceso estable con la siguiente característica (Mismo ejemplo visto en la entrada del predictor de Smith):

El controlador PI es hecho colocando el parámetro Ti igual al tao de la planta, osea Ti=6. Determinamos entonces que deseamos que la respuesta en lazo cerrado sea un 75% más rápido que la respuesta en lazo abierto asi . De esa forma obtenemos el siguiente controlador PI:

La respuesta en lazo cerrado viene dado por:

Sabemos que la dinámica de rechazo de perturbación viene dada por:

reemplazamos los valores que conocemos. Ya tenemos el lazo cerrado y tenemos la dinamica de la planta en lazo abierto:

Ahora diseñamos el filtro cumpliendo las dos condiciones necesarias. Para este caso diseñamos un filtro de segundo orden, para cancelar la dinámica de lazo cerrado.

Ya tenemos el valor de el cual cancela la dinámica de lazo cerrado. Los polos del filtro los definimos como , ese parámetro es de libre elección, en este caso yo lo escogí más rápido que la dinámica en lazo cerrado. Con esos valores procedemos a encontrar

De esa forma, el filtro del predictor de Smith, viene dado por:

Podemos notar como el rechazo de perturbaciones es acelerado gracias a la inclusión del filtro, el cual consigue cancelar la dinámica de la planta en lazo abierto:

Respuesta del Sistema:

En la grafica de la salida se compara el comportamiento del Predictor de Smith y el Predictor de Smith Filtrado, puede verse como el PSF acelera la dinámica de rechazo de perturbación. El código de implementación está abajo al final del post, para que lo puedas copiar e implementar en tu propio computador.

 Ejemplo caso Integrador (igual para caso Inestable)

Para entender mejor el concepto de sintonia del filtro, se presentará a continuación un ejemplo. Se asume un proceso integrador, propuesto por el siguiente modelo nominal :

Se discretiza el proceso con un retenedor de orden cero (ZOH) y se obtiene el siguiente sistema: , obtenido con s, y se desea diseñar un filtro para rechazar perturbaciones del tipo escalón.

Dividimos el modelo del proceso en terminos estables y terminos inestables como:

Parte estable=
Parte inestable=
Numerador=

desmembrabdo el bloque como numeradores y denominadores de cada uno de sus terminos, dicho bloque puede ser reescrito como:

El ingeniero tiene la libertad de escoger arbitrariamente el denominador del filtro con raíces dentro del circulo unitario con un orden igual o mayor al numero de polos NO deseados del sistema. En este ejemplo dado que tenemos tan solo «1» polo indeseado en el origen, vamos a escoger el  filtro de orden uno. Para eso seleccionamos el polo del filtro igual a o,9. Así, todos los polos de lazo abierto del modelo del proceso que esten fuera del circulo de radio 0,9 serán cancelados por el filtro predictor garantizando con esto la eliminación del polo en el origen.

El polo del filtro entonces es seleccionado como de manera de no ser mas lento que la dinamica deseada para el rechazo de perturbaciones.

La idea es que la parte del numerador, , del bloque se cancele con la parte del denominador inestable, , del mismo bloque. Reemplazando valores tenemos que:

Reescriviendo la ecuación en terminos de potências positivas, se tiene que

El polinomio es un residuo el cual debe ser seleccionado de forma que complete el orden del polinomio de la izquierda. Como el denominador del filtro fue escogido de primer orden, el numerador del filtro también debe ser de primer orden, de forma que el filtro sea propio, osea, . Eso resulta en la igualdad:

Igualando termino a termino y solucionando el sistema de ecuaciones, se llega a que ; ; ; ; ; e . Así, el filtro del preditor es dado por

Así el bloque es dado por:

Para este proceso se sintoniza un controlador PI con los siguientes parametros:

y  Obteniendo un controlador como:

y en su version discreta por tustin

Esquema en Simulink:

Resultado:

Espero esta información te haya sido de utilidad, no olvides compartir para de esa manera permitir que más personas se beneficien de esta información.

Códigos en MATLAB:

A continuación están todos los codigos del predictor de Smith Filtrado implementados en MATLAB y SIMULINK para que tu mismo los puedas reproducir en tu casa y aprendas una forma de hacerlo. Para tener acceso a los códigos solo debes compartir el contenido de este post, de esa forma ayudas a que otras personas conozcan sobre esta información y además ayudas a que este sitio web continue aportando contenido gratuito y de calidad. Gracias!

Código para el Caso Estable:

% Predictor de Smith
% Sergio Andres Castaño Giraldo
% http:\\controlautomaticoeducacion.com
%_________________________________________________________________________

clc
clear all
close all

Ts=0.4;             %Tiempo de Muestreo
FS=18;
%% Proceso Real P(s)
nP=0.12;                %Numerador de la Planta
dP=[6 1];               %Denominador de la Planta
P=tf(nP,dP);            %Processo real P(s)
L=3.0;                  %Retardo de P(s)
P.iodelay=L;            %Aplico el retardo en P(s)
Pd=c2d(P,Ts);           %Discretizo el Proceso con retendedor de Orden Cero P(z)
[B,A]=tfdata(Pd,'v');   %Divido el proceso discreto en numerador y denominador
d=Pd.iodelay;          %Obtengo el retardo en discreto

%% Modelo Nominal Pn(s)
nPn=0.12;                %Numerador del Proceso Nominal
dPn=[6 1];               %Denominador del Proceso Nominal
Gn=tf(nP,dP);            %Modelo Nominal sin retardo Gn(s)
Ln=3.0;                  %Retardo de Pn(s)
Pn=Gn;
Pn.iodelay=Ln;           %Aplico el retardo en P(s)
Pdn=c2d(Pn,Ts);           %Discretizo el Proceso nominal con retendedor de Orden Cero Pn(z)
Gnd=Pdn;
Gnd.iodelay=0;          %Modelo Nominal sin retardo Gn(z)
[Bgn,Agn]=tfdata(Gnd,'v');  %Divido el modelo discreto en numerador y denominador
[Bn,An]=tfdata(Pd,'v');  %Divido el proceso discreto en numerador y denominador
dn=Pdn.iodelay;           %Obtengo el retardo en discreto

%% Controlador Primario por cancelamiento de polos
Ti=dPn(1);      %Selecciono el parametro Ti=Tn (Parametro Integral)
Kn=nPn;         %Ganancia de la planta
Tr=dPn(1)*0.75; %Constante de tiempo que deseo en lazo cerrado (75% del tao de la planta)
Kc=Ti/(Kn*Tr);  %Parametro proporcional del controlador PI
C=tf(Kc*[Ti 1],[Ti 0]); %Controlador PI [C(s)]
Cd=c2d(C,Ts,'tustin');  %Discretizo el controlador
[Bc,Ac]=tfdata(Cd,'v');   %Divido el controlador discreto en numerador y denominador

%% Filtro de Robustez
Tr=4.5;    %Dinamica en lazo Cerrado
Tf=3;       %Dinamica del Filtro
beta=(1-(1-Tf/Ti)^2*exp(-Ln/Ti))*Ti;
Frs=tf(conv([Tr 1],[beta 1]),conv([Tf 1],[Tf 1]));
Fr=c2d(Frs,Ts,'tustin');
%Fr=filtro_siso(Pdn,0.5,0.91,1); %Calculo del filtro de robustez del PSF
[Bf,Af]=tfdata(Fr,'v');   %Divido el filtro discreto en numerador y denominador

%% Simulacion
nit=600;        %Numero de Interacciones
ref(1:nit)=0;   %Vector de SetPoint
ref(27:end)=1;  %Escalon unitario en el Set Point
rq(1:nit)=0;
rq(220:nit)=10;   %Estimulo de perturbacion de carga
rn(1:nit)=0;
rn(450:nit)=0.3; %Estimulo da perturbacion de salida

%Variables de simulacion
y1(1:nit) = 0; y(1:nit) = 0; u(1:nit) = 0;e(1:nit) = 0;yp(1:nit) = 0;
yt(1:nit) = 0; yp(1:nit) = 0; ep(1:nit) = 0;yn(1:nit) = 0; up(1:nit) = 0;
yf(1:nit) = 0; %Salida Filtrada
%Loop de Control
for i=1:2
for k=d+5:nit
    
    %Salida del proceso real P(z) [y(t)]
    y1(k)= B(1)*up(k-1-d)+B(2)*up(k-2-d)-A(2)*y1(k-1);
    y(k)=y1(k)+rn(k);
    
    %Salida del proceso nominal Pn(z)
    yn(k)= Bn(1)*u(k-1-dn)+Bn(2)*u(k-2-dn)-An(2)*yn(k-1);
    
    %Diferencia entre proceso real y proceso nominal (error de prediccion)
    ep(k)=y(k)-yn(k);
    
    if i==2 %Si es igual a 2 utilice el filtro del PSF
       yf(k)= Bf(1)*ep(k-1)+Bf(2)*ep(k-2)+Bf(3)*ep(k-3)-Af(2)*yf(k-1)-Af(3)*yf(k-2);
    end
    
    %Salida del modelo rapido Gn(z)  [y(k+d)]
    yt(k)= Bgn(1)*u(k-1)+Bgn(2)*u(k-2)-Agn(2)*yt(k-1);
    
    %Suma entre error de prediccion y salida predicha
    if i==1
        yp(k)=yt(k)+ep(k);
    else
        yp(k)=yt(k)+yf(k);
    end
    
    %Error de la ley de control
    e(k)=ref(k)-yp(k);
    
    %ley de control
    u(k)=u(k-1)+ Bc(1)*e(k)+Bc(2)*e(k-1);  
    
    %Ingreso de perturbacion en la entrada
    up(k)=u(k)+rq(k);
end
if i==1
    ySP=y;
    uSP=u;
else
    ySPF=y;
    uSPF=u;
end
end

%Grafico
t = 0:Ts:(nit-1)*Ts;
figure(1)
subplot(3,1,1);
hold on;
plot(t,ref,t,ySP,t,ySPF,'Linewidth',3)
legend('Referencia','PS','PSF');
xlabel('Tiempo (s)','FontSize',FS);
ylabel('Salida','FontSize',FS);
set(gca,'fontsize',FS);
grid on;
hold on
subplot(3,1,2);
hold on;
plot(t,uSP,t,uSPF,'Linewidth',3)
xlabel('Tiempo (s)','FontSize',FS);
ylabel('Control','FontSize',FS);
set(gca,'fontsize',FS);
grid on;
subplot(3,1,3);
hold on;
plot(t,rq,t,rn,'Linewidth',3)
legend('Perturbacion de carga','perturbacion de salida');
xlabel('Tiempo (s)','FontSize',FS);
ylabel('Salida','FontSize',FS);
set(gca,'fontsize',FS);
grid on;
hold on

Código para el Caso Integrador y Inestable

>> Descarga el código del Simulink y de Matlab <<

% Predictor de Smith Filtrado
% Sergio Andres Castaño Giraldo
% http:\\controlautomaticoeducacion.com
%_________________________________________________________________________

h=0.2;                        % Tiempo de Muestreo
tsim=140;                     % tiempo de simulación
Fs=18;                        % Tamaño de la letra del plot

%% Descripción del Sistema

% Modelo nominal del proceso
Pn=tf(3.5,[1 0]);
Gn=Pn;

% Modelo nominal del proceso con retardo
Pn.iodelay=1;


%Discretización del Sistema
Pnd=c2d(Pn,h);                    % Modelo discreto
Gnd=c2d(Gn,h);   

%Proceso Real (Caso Ideal)
Pd=Pnd;

%% Parametros del controlador

k11f=0.14; t11f=10;              % Parametros del PI
C=k11f+k11f/t11f*tf([h 0],[1 -1],h);

%--------------- Calculo del filtro  ------------------------------ 
Fr=tf([1.5 -1.4],[1 -0.9],h);

% Calculo del Filtro estable S
S=Gnd-Fr*Pnd;

%Simula utilizando SIMULINK
sim('PSF_Int');

%Grafico
figure(1)
subplot(2,1,1);
plot(t,Q(:,1),t,Q(:,2),'Linewidth',3)
legend('Referencia','PSF','PSF');
xlabel('Tiempo (s)','FontSize',FS);
ylabel('Salida','FontSize',FS);
set(gca,'fontsize',FS);
grid on;
subplot(2,1,2);
hold on;
plot(t,Q(:,3),'Linewidth',3)
xlabel('Tiempo (s)','FontSize',FS);
ylabel('Control','FontSize',FS);
set(gca,'fontsize',FS);
grid on;

BIBLIOGRAFIA

1. S. A. C. Giraldo, R. C. C. Flesch, and J. E. Normey-Rico, “Multivariable greenhouse control using the filtered  Smith predictor,” Journal of Control, Automation and Electrical Systems, pp. 1–10, 2016. (click aca)
2. S. A. C. Giraldo, R. C. C. Flesch, J. E. Normey-Rico, and M. Z. P. Sejas, “Decoupling filtered smith predictor design for multivariable systems with multiple time delays,” in 2016 12th IEEE International Conference on Industry Applications (INDUSCON), Nov 2016, pp. 1–8. (Click aca) 
3. S. A. C. Giraldo, “ESTUDO DE TECNICAS DE SINTONIA DO PREDITOR DE SMITH FILTRADO PARA SISTEMAS MULTIVARIAVEIS COM ATRASO ”  (Portugués) Mestrado, Florianópolis, 2016. (Click aca)
4. J. E. Normey-Rico and E. F. Camacho, Control of dead-time processes. London: Springer, 2007. (click aca)