Control por asignación de polos (RST Incremental)

Hola Controleros y Controleras, hoy vamos a aprender como funciona un Control RST Incremental (acción integral dentro del lazo de control). Recordando que el control RST es un método elegante de como usar una colocación de polos.

Antes de comenzar, podrías darle un vistazo al:

CURSO GRATUITO DE CONTROL DE PROCESOS AVANZADOS

Tabla de Contenido

1 Control RST Incremental
2 Proyecto de Control RST Incremental
- 2.1 Ley de Control
3 Ejemplo Control RST Incremental
- 3.1 Ejemplo 1
- 3.2 Ejemplo 2
  - 3.2.1 Código en Matlab:

Control RST Incremental

El control RST Incremental es basicamente el mismo control RST posicional visto anteriormente, solo que presenta una pequeña modificación en el lazo de control que altera la ley de control y altera la identidad polinomial de la siguiente manera.

Por eso será importante que veas esa entrada anterior, donde explicamos en detalle que es un RST y mostramos algunos aspectos teoricos del controlador.

A continuación vamos a ver la modificación (adición de una acción integradora) en el modelo de la planta del proceso y en la acción de control:

Planta:

A(z^{-1})\Delta y(t)=z^{-d}B(z^{-1})\Delta u(t)

Donde:

\Delta=1-z^{-1}\rightarrow \Delta u(t)-u(t-1)

Ley de Control

R(z^{-1})\Delta u(t)=T(z^{-1})y_i(t)-S(z^{-1})y(t)

Con esto, el diagrama de bloques del sistema en lazo cerrado queda de la siguiente manera:

Notese que en el Lazo de control se adiciono un Integrador que acompaña al polinomio R y es denotado como Δ, este termino garantiza que el sistema siga la referencia y también que pueda sobrellevar perturbaciones del sistema como el do(t) que esta ingresando a la salida de la planta.

Identidad:

La identidad polinomial con el ingreso del termino integrador en el controlador RST (Reference Signal Tracking) o simplemente Seguidor de Referencia viene dado por:

A(z^{-1})\Delta R(z^{-1})+z^{-d}B(z^{-1})S(z^{-1})=P_{MF}

Tamaño del Polinomio R

N_R=N_B+d-1

Tamaño del Polinomio S

N_S=N_A

$N_B$ Tamaño del Polinomio B Numerador

$N_A$ Tamaño del Polinomio A Denominador

Control PID con microcontrolador PIC

Controlador PI por Asignación de Polos

Enmascarar control RST en control PID

Proyecto de Control RST Incremental

Acomodando la ley de control con el modelo estimado de la planta se obtiene:

y(t)=\dfrac{z^{-d}BT}{A\Delta R+z^{-d}BS}Y_i(t)=\dfrac{z^{-d}BT}{P_{MF}(z^{-1})}

Para garantizar el seguimiento del set point el polinomio T debe ser igual a la suma de todos los coeficientes del polinomio S

T(1)=S(1)

Para un sistema de Segudo Orden por ejemplo, la identidad polinomial quedaria de la siguiente manera:

Gm(z)=\dfrac{b_0+b_1z^{-1}}{1+a_1z^{-1}+a_2z^{-2}}z^{-1}

n_a=2;n_b=1;d=1

n_r=n_b+d-1=1+1-1=1\rightarrow (R_0+R_1z^{-1})

n_s=n_a=2\rightarrow (S_0+S_1z^{-1}+S_2z^{-2})

(1+a_1z^{-1}+a_2z^{-2})(1-z^{-1})(R_0+R_1z^{-1})+z^{-1}(b_0+b_1z^{-1})(S_0+S_1z^{-1}+S_2z^{-2})=(1+p_1z^{-1}+p_2z^{-2})

P1 y P2 son los polos asignados.

Si multiplicamos los terminos del denominador (A) con el integrador (Δ) tenemos:

\bar{A}=A(z^{-1})\Delta=1+\bar{a_1}z^{-1}+\bar{a_2}z^{-2}

\bar{A}=A(z^{-1})\Delta=(1+a_1z^{-1}+a_2z^{-2})(1-z^{-1})

\bar{A}=A(z^{-1})\Delta=1+(a_1-1)z^{-1}+(a_2-a_1)z^{-2}-a_2z^{-3}

\begin{bmatrix} 1 &b_0 &0 &0 \\ \bar{a_1} &b_1 &b_0 &0 \\ \bar{a_2} &0 &b_1 &b_0 \\ \bar{a_3} &0 &0 &b_1 \end{bmatrix}\begin{bmatrix} R_1\\ S_0\\ S_1\\ S_2 \end{bmatrix}=\begin{bmatrix} p_1-\bar{a_1}\\ p_2-\bar{a_2}\\ -\bar{a_3}\\ 0 \end{bmatrix}

T=S_0+S_1+S_2

Recordar que $R_0=1$

Ley de Control

R(z^{-1})\Delta u(t)=T(z^{-1})y_i(t)-S(z^{-1})y(t)

(1+R_1)(1-z^{-1})u(t)=TY_i(t)-(S_0+S_1z^{-1}+S_2z^{-2})y(t)

Despejando U y aplicando transformada inversa Z

u(k)=TY_i(k)-S_0y(k)-S_1y(k-1)-S_2y(k-2)-(R_1-1)u(k-1)-R_1u(k-2)

Ejemplo Control RST Incremental

Ejemplo 1

Usando el ejemplo numero 2 visto en la entrada anterior, se hace el mismo tratamiento para un RST incremental usando las ecuaciones anteriormente descritas que son para un sistema de segundo orden. Se presenta entonces el codigo en matlab.

clc
clear all
close all
%% Simulacion
% Condições iniciais
clear all;  clc; close all;
nit = 400; umin = 0; umax = 4.9; ts = 0.1;
u(1:10) = 0; erro(1:10) = 0; y(1:10) = 0;
% Referência
r(1:100) = 1; r(101:200) = 4; r(201:300) = 3; r(301:nit) = 3;
%% Planta
T=0.1;                          %Tempo da Amostragem
num=1;                          %Numerador Continuo
den=[0.1 1.1 1];                %Denominador Continuo
gp=tf(num,den);                 %Função de transferencia
ftz=c2d(gp,T);                  %Planta Discreta
[B,A]=tfdata(ftz,'v');          %Divide Numerador em B e denominador em A
nb=length(B);                   %Ordem do numerador
na=length(A)-1;                 %Ordem do Denominador
nr=nb-1;                        %Ordem do polinomio R
ns=na;                        %Ordem do Polinomio S

gpz=filt(B,A,T,'name','Planta');      %Funcao de transferenca discreta

%% PMF
%Alocaçao de polos por condições de disenho.
Ts=1;                     %Tempo de estabelecimento desejado
ep=0.6;                     %Epsilon 
wn=4/(ep*Ts);               %

%Polos conjugados
za=exp(-ep*wn*T);
th=57.3*wn*T*sqrt(1-ep^2);
polos=za*[cosd(th)+i*sind(th),cosd(th)-i*sind(th)];

% Polinomio desejado, com polo insignificante igual a zero
Am=poly([polos 0 0]);

%Delata A. Multiplica por integrador
DA=conv(A,[1 -1]);

%Resolução de Equações metodo matricial
M=[1 B(2) 0 0; DA(2) B(3) B(2) 0; DA(3) 0 B(3) B(2);DA(4) 0 0 B(3)];
q=[Am(2)-DA(2); Am(3)-DA(3);Am(4)-DA(4);Am(5)];
X=inv(M)*q;

R=[1 X(1)];             %Polinomio R
S=X(2:end)';            %Polinomio S

%Polinomio T
Tz=sum(S);

%Ingresa la perturbacion en el instante 300
d0(1:300)=0;
d0(301:nit)=0.05;

for k = 5:nit
    %Salida de la Planta
    y(k)=B(2)*u(k-1)+B(3)*u(k-2)-A(2)*y(k-1)-A(3)*y(k-2)+d0(k);
    %Ley de Control
    u(k)=Tz*r(k)-S(1)*y(k)-S(2)*y(k-1)-S(3)*y(k-2)-(R(2)-1)*u(k-1)+R(2)*u(k-2);
    %     Saturação
    if u(k) >= umax 
     u(k) = umax;
     elseif u(k) <= umin
     u(k) = umin;
    end
end

t = 0:ts:(nit-1)*ts;

figure
subplot(2,1,1),plot(t,r,'--k',t,y,'-r','Linewidth',2.5)
xlabel('Tempo (s)');
ylabel('Velocidade');
legend('y_r','y','Location','SouthEast')
grid on;
hold
subplot(2,1,2),plot(t,u,'--b','Linewidth',3)
xlabel('Tempo (s)');
ylabel('Controle (volts)');
legend('u')
grid on;

Notese que a los 30 segundos ingresa una perturbación al sistema, y el controlador RST incremental es capaz de reaccionar para volver a llevar la variable al punto de referencia, cosa que no sucedía con el control RST Posicional.

Ejemplo 2

Considere que el sistema térmico del laboratorio de temperatura vista en otra entrada viene representada por la siguiente función de transferencia:

G(s)=\dfrac{1.04}{160s+1}e^{-10s}

Cuya representación discreta con $T_s=8$ viene dada por:

G(z)=\dfrac{b_0 + b_1z^{-1}}{1 – a_1z^{-1}}z^{-d}=\dfrac{0.03828 + 0.01244z^{-1}}{1 – 0.9512z^{-1}}z^{-2}

El polinomio R viene dado por:

n_r=n_b+d-1=1+2-1=2

R_0+R_1z^{-1}+R_2z^{-2}

El polinomio S viene dado por:

n_s=n_a=1

S_0+S_1z^{-1}

El polinomio $\Delta A$

(1+a_1z^{-1})(1-z^{-1})=1+(a_1-1)z^{-1}-a_1z^{-2}=1+\bar{a_1}z^{-1}+\bar{a_2}z^{-2}

Resolviendo la parte izquierda de la identidad polinomial del control RST:

A(z^{-1})\Delta R(z^{-1})+z^{-d}B(z^{-1})S(z^{-1})

(1+\bar{a_1}z^{-1}+\bar{a_2}z^{-2})(R_0+R_1z^{-1}+R_2z^{-2})+z^{-2}(b_0 + b_1z^{-1})(S_0+S_1z^{-1})

Ecuación E2.1

R_0+(R_1+\bar{a_1}R_0)z^{-1}+(\bar{a_1}R_1+R_2+b_0S_0+\bar{a_2}R_0)z^{-2}+(\bar{a_1}R_2+\bar{a_2}R_1+b_0S_1+b_1S_0)z^{-3}+(b_1S_1+\bar{a_2}R_2)z^{-4}

Se establece por medio del método de asignación de polos que el sistema tenga el comportamiento dominante de un sistema subamortiguado. En este caso se debe establecer un polinomio deseado de cuarto orden, se establece un tiempo de establecimiento deseado de $T_{ss}=220$ y un máximo pico de 5% para los dos polos dominantes y se completa con 2 polos insignificantes (rápidos).

M_p=100e^{\dfrac{-\pi \zeta}{\sqrt{1-\zeta^2}}}

5=100e^{\dfrac{-\pi \zeta}{\sqrt{1-\zeta^2}}}

\zeta=0.6901

\omega_n=\dfrac{4}{\zeta T_{ss}}=0.0263

Ecuación de los polos:

s_{1,2}=-\zeta\omega_n\pm j\omega_n\sqrt{1-\zeta^2}

Los polos deseados son (los dos primeros son los polos dominantes dados por la ecuación anterior, y los dos últimos son los polos rápidos):

S_d=[-0.0182 \pm j0.0191, -10.0000, -12.0000]

En tiempo discreto viene dado por:

Z_d=e^{S_dTs}

Z_d=[0.8546 \pm j0.1314,0,0]

Cuyo polinomio deseado es (Ecuación E2.2):

P_{MF}=1 -1.7092z^{-1}+ 0.7476z^{-2}+0z^{-3}+0z^{-4}=1 +\alpha_1z^{-1}+ \alpha_2z^{-2}+\alpha_3z^{-3}+\alpha_4z^{-4}

Igualando la Ecuación E2.1 con la Ecuación E2.2 y expresandolo en forma matricial (matriz sylvester), tenemos que:

\begin{bmatrix} 1 & 0&0&0 &0 \\ \bar{a_1} &1&0&0 &0 \\ \bar{a_2} &\bar{a_1}&1&b_1 &0 \\ 0 &\bar{a_2}&\bar{a_1}&b_2 &b_1\\ 0 &0&\bar{a_2}&0 &b_2 \end{bmatrix}\begin{bmatrix} R_0\\ R_1\\ S_0\\ S_1\\ S_2 \end{bmatrix}=\begin{bmatrix} 1\\ \alpha_1\\ \alpha_2\\ \alpha_3\\ \alpha_4 \end{bmatrix}

El polinomio T (precompensador estático)

T=S_0+S_1

La ley de control viene dado por:

u(k)=TY_i(k)-S_0y(k)-S_1y(k-1)-(R_1-R_0)u(k-1)-(R_2-R_1)u(k-2)+R_2u(k-3)

Código en Matlab:

%% Control RST (Reference Signal Tracking)
% By: Sergio Andres Castaño Giraldo
% https://controlautomaticoeducacion.com/

clc
clear 
close all

Ts=8; %Periodo de Muestreo
%% Modelo del Proceso
k=1.04;tau=160;theta1=10;
theta = theta1 + Ts/2;

%Función de Transferencia
num = k;
den = [tau 1];
G=tf(num,den);
G.iodelay = theta1;

%Discretización
ftz=c2d(G,Ts);                  %Planta Discreta
[B,A]=tfdata(ftz,'v');          %Divide Numerador em B e denominador em A
d=ftz.iodelay;                  %Retardo Discreto
nb=length(B);                   %Orden del numerador
na=length(A)-1;                 %Orden del Denominador
nr=nb+d-1;                        %Orden del polinomio R
ns=na;                          %Ordem del Polinomio S

%% Especificaciones de Diseño
Mp=5; %Maximo Pico
ep=sqrt(((log(Mp/100))^2)/(pi^2+((log(Mp/100))^2))); %Fator de amortiguamiento
Tlc=220; 
if ep<1
    Wn=4/(ep*Tlc)
else if ep==1
        Wn=5.8335/Tlc;
    end
end

Sd=[-ep*Wn+1i*Wn*sqrt(1-ep^2), -ep*Wn-1i*Wn*sqrt(1-ep^2), -10, -12]; %Alocação de Polos
Zd = exp(Sd*Ts);
Am = poly(Zd);

%Delata A. Multiplica por integrador
DA=conv(A,[1 -1])';

%Resolução de Equações metodo matricial
M=sylvester_RST(B,A,d,1);
P=Am';
% P(2:length(DA)) = P(2:length(DA))- DA(2:end);

X=inv(M)*P;

R=X(1:nr);             %Polinomio R
S=X(nr+1:end);            %Polinomio S
%Polinomio T
Tz=sum(S);

%% inicializa parametros de Simulacion
T0 = 25;  %Temperatura Ambiente
nit=120;    %Numero de interacciones
u(1:nit) = 0; 
ym(1:nit) = 0;  
r(1:nit) = T0;

% Referência
r(10:80) = 45; r(81:nit) = 30; 

for k=4:nit
        %Planta
        %ym(k)=B(1)*u(k-d)+B(2)*u(k-1-d)-A(2)*ym(k-1);   
        t = 0:Ts:(k-1)*Ts;
        ym=lsim(G,u(:,1:k),t,'zoh')+T0;
        %Ley de Control
        u(k)=Tz*r(k)-S(1)*ym(k)-S(2)*ym(k-1)-(R(2)-R(1))*u(k-1)-(R(3)-R(2))*u(k-2)+R(3)*u(k-3);

        if u(k) > 100
            u(k) = 100;
        else if u(k)<0
                u(k) = 0;
            end
        end

end


t = 0:Ts:(nit-1)*Ts;
figure
subplot(2,1,1)
stairs(t,r,'--k','Linewidth',2),hold on
stairs(t,ym,'-r','Linewidth',2)
xlabel('Tiempo (s)');
ylabel('Temperatura (C)');
legend('r','y','Location','SouthEast')
grid on;
hold
subplot(2,1,2)
stairs(t,u,'b','Linewidth',2)
xlabel('Tiempo (s)');
ylabel('Heater (%)');
legend('u')
grid on;

Función para el cálculo de la matriz de Sylvester (Guardar en la misma carpeta del código principal)

function M = sylvester_RST(B,A,d,i)
% Sylvester matrix for RST Control:
%    M = sylvester_RST(B,A,d,i)
% M = Sylvester matrix
% B = transfer function numerator (discrete)
% A = transfer function denominator (discrete)
% d = discrete time delay
% i = 1 for incremental form, 0 for positional form
%
% By: Sergio Andres Castaño Giraldo
% https://controlautomaticoeducacion.com/
%


    nb=length(B)-1;                 %Ordem do numerador
    na=length(A)-1;                 %Ordem do Denominador
    nr=nb+d;                        %Ordem do polinomio R
    
    %Delata A. Multiplica por integrador
    if i==1
        DA=conv(A,[1 -1]);
        
        ns=na+1;                        %Ordem do Polinomio S
    else
        DA=A;
        ns=na;                        %Ordem do Polinomio S
    end
    nda=length(DA);
    
    LM=nr+ns; %Cumprimento da matriz
    M=zeros(LM,LM);
    M(1:nda,1)=DA';
    M(1+d:length(B)+d,nr+1)=B';
    
    %Polinomio R
    for k=2:nr
        M(k:end,k)=M(1:end-k+1,1);         % Forma o polinomio R             
    end
    
    %Polinomio S
    for k=1:ns-1
        M(k+1:end,nr+1+k)=M(1:end-k,nr+1);         % Forma o polinomio S                
    end
end

La respuesta del sistema viene dado por:

Control RST — Respuesta del control RST en el Laboratorio de Temperatura TCLAB

Eso es todo por la entrada del dia de hoy, espero les haya gustado y hayan aprendido algo nuevo. Si te ha servido el contenido de esta entrada, de los videos y los códigos de implementación y deseas apoyar mi trabajo invitandome a un café super barato, puedes hacerlo en el siguiente link:

👉 Invitar a Sergio a un Café ☕️

Que esten muy bien, nos vemos en la siguiente entrada.

Volver al Curso de Control de Procesos Avanzados

Mi nombre es Sergio Andres Castaño Giraldo, y en este sitio web voy a compartir una de las cosas que mas me gusta en la vida y es sobre la Ingeniería de Control y Automatización. El sitio web estará en constante crecimiento, voy a ir publicando material sobre el asunto desde temas básicos hasta temas un poco más complejos. Suscríbete al sitio web, dale me gusta a la página en Facebook y únete al canal de youtube. Espero de corazón que la información que comparto en este sitio, te pueda ser de utilidad. Y nuevamente te doy las gracias y la bienvenida a control automático educación.

Control RST Incremental

Planta:

Ley de Control

Identidad:

Proyecto de Control RST Incremental

Ley de Control

Ejemplo Control RST Incremental

Ejemplo 1

Ejemplo 2

Código en Matlab:

Comparte esto: