Control Predictivo

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:

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.

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]

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.

Ejemplo

Comparte esto: