Saltar al contenido
Control Automático Educación

GPC con Perturbaciones

En esta entrada vamos a aprender como involucrar el modelo de perturbaciones conocidas en la estructura del GPC y con esto podriamos pensar en realizar con control con acción FeedForward.

 

Para el caso de un control convencional tenemos la siguiente representación:

Control Convencional

Si tengo un sensor que mida la entrada q, podría hacer un controlador para ayudar al control principal para rechazar esa perturbación. En otras palabras podría realizar un control FeedForward.

Control FeedForward

El problema de esta abordaje, es que muchas veces ese bloque FeedForwrd muchas veces no es realizable, por ejemplo si el sistema G fuera de fase no mínima, el sistema se volvería inestable, o si tuviera más polos en G que en Gq, va a quedar con más ceros y por ende queda No realizable. Pero en términos generales, si yo consigo llevar mi sistema a esta forma, voy a poder rechazar mis perturbaciones de una mejor forma

 

En el control GPC, tendremos que modificar el modelo CARIMA para poder incluir el efecto de la perturbación. Esto quiere decir que básicamente el control GPC no cambia, solo que ahora con ayuda de mi modelo, el controlador tendrá que efectuar las predicciones teniendo en consideración la presencia de una perturbación medible dentro del lazo de control.

 

Ay(k)=Bu(k-1)+\dfrac{C}{\Delta}\varepsilon(k)+Dq(k)

 

Veamos un Ejemplo con dos funciones de transferencia, una para la perturbación y otra para la planta:

 

G_q=\dfrac{1-0.9z^{-1}}{1-0.8z^{-1}}

 

 

G=\dfrac{0.5z^{-1}}{1-0.5z^{-1}}

Vemos que como en este caso, los denominadores NO son iguales, debemos obtener un polinomio común para que sea nuestro polinomio A, asi los polinomios A,B y C quedan como:

A=(1-0.5z^{-1})(1-0.8z^{-1})=1-1.3z^{-1}+0.4z^{-2}

 

 

B=0.5(1-0.8z^{-1})=0.5-0.4z^{-1}

 

 

D=(1-0.5z^{-1})(1-0.9z^{-1})=1-1.4z^{-1}+0.45z^{-2}

 

 

Multiplicando el modelo CARIMA por E_j las predicciones pasan a ser>

 

 

E_j\tilde{A}y(k+1)=E_jB\Delta u(k-1+j)+E_jD\Delta q(k+j)+E_jC\varepsilon(k+j)

 

 

[1-z^{-j}F_j]y(k+1)=E_jB\Delta u(k-1+j)+E_jD\Delta q(k+j)+E_jC\varepsilon(k+j)

 

 

\hat{y}(k+1|k)=F_jy(k)+E_jB\Delta u(k-1+j)+E_jD\Delta q(k+j)

 

 

El término E_jD\Delta q(k+j) puede contener información de perturbaciones pasadas y de perturbaciones futuras. Veamos que la perturbación es algo conocido, si bien puede que yo no conozca las perturbaciones futuras, eventualmente voy a terminar conociéndolas cuando entren al sistema, así entonces yo colocaría las perturbaciones futuras igual a cero (En el caso de no conocerlas) y voy a tener valores de perturbación pasada. Osea que las perturbaciones van a pertenecer únicamente a la respuesta FORZADA del sistema.

 

Sabemos que por definición la respuesta forzada es la respuesta provocada por una variación de la señal de control. Osea que las perturbaciones NO van a provocar nada en la variación de control, por eso es respuesta libre.

 

 

E_jD\Delta q(k+j)=L\Delta q(k+j)+F_q

 

 

Donde

 

 

L\Delta q(k+j) son los terminos del futuro (Si no se conocen entonces se coloca \Delta q(k+j)=0)

 

F_q son los terminos del pasado.

 

L es igual que el G solo que con los coeficientes de la perturbación con una excitación del tipo escalón

 

Sabemos que las predicciones viene dado por:

Y=G\Delta u +\tilde{f}
\tilde{f}=f+L\Delta q(k+j)+F_q

 

f es la respuesta libre normal del MPC

 

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




 

Ejemplo

Proyectar un control predictivo GPC teniendo en consideración la perturbación medible del sistema. La función de transferencia del proceso, cuanto de la perturbación se describen a contrinuación:

 

G_p=\dfrac{0.2713z^{-3}}{1- 0.8351z^{-1}}

 

G_q=\dfrac{0.05z^{-3}}{1- 0.91z^{-1}}

 

Código de implementación:

Recuerda que para poder visualizar el código de implementación, debes compartir el contenido de este post con cualquiera de los siguientes 3 botones, así ayudas a que más personas se beneficien de esta información.

[sociallocker id=948]

%% Control GPC
% Sergio Andres Castaño G.
% controlautomaticoeducacion.com
%___________________________________________________________________
clc
clear all 
close all
p= menu('Perturbacion','Conocida','No Conocida');
T=2;
Ap=[1 -0.835];
Bp=[0 0.271];
Dq=[0 0.05];
Aq=[1 -0.9];
dm=2;
dq=2;
fz=filt(Bp,Ap,T);
fz.iodelay=dm
fzq=filt(Dq,Aq,T);
fzq.iodelay=dq
%Minimo Comun Multiplo de los polinomios A,B,D
D=conv(Ap,Dq);
A=conv(Ap,Aq);
B=conv(Bp,Aq);
%% Parametros de sintonia del GPC
N=10;        %Horizonte de Prediccion
N1=dm+1;    %Horizonte Inicial
N2=dm+N;     %Horizonte final
Nu=5;       %Horizonte de Entrada
lambda=1;   %Parametro de ponderacion del Control
delta=1;    %Parametro de ponderacion del seguimiento a referencia
%Matrices de Ponderacion
Ql=eye(Nu)*lambda;
Qd=eye(N)*delta;
%% Calculo de la ecuacion Diofantina
[En,F] = diophantine(A,N2,0); %Calculo de la funcion Diofantina
E=En(end,:); %Polinomio E seria la ultima fila arrojada por la funcion
F=F(N1:N2,1:end);
%% Determina los coeficientes de control pasados
%Elimino el cero de la primera posicion de B
if B(1)==0
B=B(2:end); 
end
uG=zeros(N,N1); %Vector de controles pasados (Pertenece a la respuesta libre)
Fq=zeros(N,N1); %Vector de perturbaciones pasadas (Pertenece a la respuesta libre)
j=2;
for i=N1:N2
aux=conv(En(i,:),B);
aux1=conv(En(i,:),D);
if length(aux) < N1+j-1  %Si la longitud del auxiliar es menor que j
aux=[aux zeros(1,(N1+j-1)-length(aux))]; %Completar con ceros
end
if length(aux1) < N1+j-1  %Si la longitud del auxiliar es menor que j
aux1=[aux1 zeros(1,(N1+j-1)-length(aux1))]; %Completar con ceros
end
uG(j-1,1:N1)=aux(j:N1+j-1);  %Controles pasados (B siempre comienza en un instante atras)
Fq(j-1,1:N1+1)=aux1(j:N1+j); %Perturbaciones pasada (D puede comenzar desde el instante actual)
j=j+1;
end
g=conv(E,B); %Calcula polinomio g
Lj=conv(D,E); %Calcula polinomio Lj
G=zeros(N,Nu);                       % Inicializa la matriz G
%Crio a Matriz L (Coeficientes Da perturbação)
L=zeros(N,Nu);
for k=1:Nu
G(k:end,k)=g(1:N-k+1);         % Forma a matriz G   
L(k:end,k)=Lj(1:N-k+1);         % Forma a matriz L 
end
Mn=inv(G'*Qd*G+Ql)*G'*Qd; %Calculo de la Funcion de Costo sin Restriccion
%Calculo do controldor K1 (Primeira fila de Mn)
K1=Mn(1,:);
%% Loop de Controle
%% inicializa parametros de Simulacion
nit=120; %Numero de interacciones
inc_u=0;
u_ant(1:10) = 0;
u(1:20+dm) = 0; ym(1:20+dm) = 0; r(1:20+dm) = 0;
% Referência
r(20:nit) = 1;
%Perturbacion
q(1:40)= 0;q(41:80) = 2;q(81:nit) = 0;
lduf=size(F,2);     %Averiguo la longitud del polinomio F
luG=size(uG,2);     %Longitud de uG
duf=zeros(1,luG);   %Vector incrementos de control pasados libre
lFq=size(Fq,2);     %Longitud de Fq
dqp=zeros(1,lFq);   %Delta de Q pasados (Entrada de perturbacion pasadas)
% Determino las magnitudes de los polinomios
na=size(A,2);
nd=size(D,2);
nb=size(B,2);
%Creo los vectores que van a almacenar las acciones pasadas
y_ant=zeros(1,na);
u_ant=zeros(1,nb);
q_ant=zeros(1,nd);
dqf=zeros(1,Nu); %Los DeltaQ futuros
for k=5+dm:nit
% Salida del proceso
ym(k)=B*u_ant'-A(2:na)*y_ant(1:na-1)'+D*q_ant';
%Actualizo vector de salidas pasadas
%y_ant=[y(k-1) y(k-2) y(k-3) ......]      
if na==1
y_ant=ym(k);
else
aux=y_ant(1:na-1);
y_ant=[ym(k) aux];
end
% Calculo de la Respuesta Libre
f1=uG*duf';        %Controles pasados
f2=(y_ant*F')';    %Salidas Pasadas
f3=Fq*dqp';        %Perturbaciones pasadas
f4=L*dqf';         %Perturbaciones Futuras
free=f1 + f2 + f3 + f4;
%Calculo do Controle
%Projeto onde nao tenho as referencias futuras
inc_u=K1*(r(k)-free);
if k==1
u(k)=inc_u;
else
u(k)=u(k-1)+ inc_u;
end
% actualiza vector de control
%u_ant=[u(k-d) u(k-1-d) u(k-2-d) ..... ]
aux=u_ant(1:nb-1);
u_ant=[u(k-dm) aux];
%actualiza duf (incrementos de control)
% duf= [du(k-1) du(k-2) ..... du(k-N)]
aux=duf(1:end-1);
duf=[inc_u aux];
%_______________________________________________________________________%    
%Actualizo vector de perturbaciones pasadas
%y_ant=[y(k-1) y(k-2) y(k-3) ......]   
if na==1
q_ant=q(k-dq);
else
aux=q_ant(1:nd-1);
q_ant=[q(k-dq) aux];
end
%Calcula el incremento de la Perturbacion
inc_q=q(k)-q(k-1);
%Calculo los incrementos de la perturbacion futura
if p==1 && k<=nit-length(dqf)
for i=1:length(dqf)
dqf(i)=q(k+i)-q(k+i-1);
end   
end
%Actualiza los incrementos de la perturbacion pasados
aux=dqp(1:end-1);
dqp=[inc_q aux];
%_______________________________________________________________________% 
end
nm=nit;
t = 0:T:(nm-1)*T;
figure(1)
subplot(2,1,1),plot(t,r,'--k',t,ym,'-r','Linewidth',2)
hold on
subplot(2,1,2),plot(t,u,'b','Linewidth',2)
hold on
grid on;
%% inicializa parametros de Simulacion
inc_u=0;
u_ant(1:10) = 0;
u(1:20+dm) = 0; ym(1:20+dm) = 0; r(1:20+dm) = 0; y1(1:20+dm) = 0;y2(1:20+dm) = 0;
% Referência
r(20:nit) = 1;
%Perturbacion
q(1:40)= 0;q(41:80) = 2;q(81:nit) = 0;
lduf=size(F,2);     %Averiguo la longitud del polinomio F
luG=size(uG,2);     %Longitud de uG
duf=zeros(1,luG);   %Vector incrementos de control pasados libre
lFq=size(Fq,2);     %Longitud de Fq
dqp=zeros(1,lFq);   %Delta de Q pasados (Entrada de perturbacion pasadas)
% Determino las magnitudes de los polinomios
na=size(A,2);
nd=size(D,2);
nb=size(B,2);
%Creo los vectores que van a almacenar las acciones pasadas
y_ant=zeros(1,na);
u_ant=zeros(1,nb);
q_ant=zeros(1,nd);
dqf=zeros(1,Nu); %Los DeltaQ futuros
for k=5+dm:nit
% Salida del proceso
ym(k)=B*u_ant'-A(2:na)*y_ant(1:na-1)'+D*q_ant';
%Actualizo vector de salidas pasadas
%y_ant=[y(k-1) y(k-2) y(k-3) ......]      
if na==1
y_ant=ym(k);
else
aux=y_ant(1:na-1);
y_ant=[ym(k) aux];
end
% Calculo de la Respuesta Libre
f1=uG*duf';        %Controles pasados
f2=(y_ant*F')';    %Salidas Pasadas
free=f1 + f2;
%Calculo do Controle
%Projeto onde nao tenho as referencias futuras
inc_u=K1*(r(k)-free);
if k==1
u(k)=inc_u;
else
u(k)=u(k-1)+ inc_u;
end
% actualiza vector de control
%u_ant=[u(k-d) u(k-1-d) u(k-2-d) ..... ]
aux=u_ant(1:nb-1);
u_ant=[u(k-dm) aux];
%actualiza duf (incrementos de control)
% duf= [du(k-1) du(k-2) ..... du(k-N)]
aux=duf(1:end-1);
duf=[inc_u aux];
%_______________________________________________________________________%    
%Actualizo vector de perturbaciones pasadas
%y_ant=[y(k-1) y(k-2) y(k-3) ......]   
if na==1
q_ant=q(k-dq);
else
aux=q_ant(1:nd-1);
q_ant=[q(k-dq) aux];
end
%_______________________________________________________________________% 
end
figure(1)
hold on
nm=nit;
t = 0:T:(nm-1)*T;
subplot(2,1,1),plot(t,ym,'-g',t,q,'--b','Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
legend('Ref','Y_{GPC-Pert}','Y_{GPC}','q(k)','Location','SouthEast')
grid on;
hold on
subplot(2,1,2),plot(t,u,'r','Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('uq','u')
grid on;

[/sociallocker]

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.