Saltar al contenido
Control Automático Educación
Control Predictivo

Control GPC

En esta entrada aprenderemos la teoría sobre el control GPC, Haré una explicación breve sobre que consiste este tipo de control, pero principalmente explicaré todo detalladamente en los 3 videos que veremos a continuación. Antes de empezar, si no eres muy experto en el tema, te recomiendo mirar inicialmente los Videos del DMC (Dando Click aqui) Debido a que la teoria que alli se aborda, es fundamental, para entender el funcionamiento del Control Predictivo GPC.

Contenidos ocultar
1. CONTROL GPC
1.1. ¿Cual es la diferencia del GPC y el DMC?
1.2. Historia
1.3. Objetivo del GPC
1.4. Modelo
1.5. Función de Costo
1.6. Control Predictivo Basado en Modelo – DMC
1.7. Control NMPC – Matlab – Simulink
1.8. MPC advanced control applied to a gas compression system of an offshore platform
2. Control Predictivo GPC
3. Código de Implementación:
3.1. Código Principal
3.2. Función Diofantina

Al final de este post, encontraras unos ejemplos prácticos de como programar un GPC junto con su respectivo video que explica claramente cada linea de código implementada.

CONTROL GPC

¿Cual es la diferencia del GPC y el DMC?

  • Para el estudio del control predictivo GPC (Generalized Predictive Controller ó Control Predictivo Generalizado) al igual que el DMC, vamos a necesitar obligatoriamente tener un modelo matemático que represente el comportamiento dinámico de mi sistema. La diferencia es que en el GPC mi modelo ya NO van a ser los coeficientes de mi respuesta al Escalón (Tal y como vimos que sucedía con el DMC) si no que ahora mi modelo va a ser basado en Función de Transferencia Discreta, con la cual vamos a calcular una predicción de lo que va a suceder en el futuro.

Historia

Para ubicarnos un poco, seria conveniente conocer que el GPC inicialmente fue propuesto en dos artículos realizados por los autores D. W. CLARKE, C. MOHTADI y P. S. TUFFS llamado:

  1. Generalized Predictive Control Algorithm Part I, the Basic
  2. Generalized Predictive Control Part II. Extensions and Interpretations

Objetivo del GPC

  • El principal objetivo del control GPC es encontrar una secuencia de control óptimo.
  • ¿Qué quiere decir con encontrar un control óptimo?
  • Por ejemplo, poder llevar mi variable al lugar que yo deseo, variando mínimamente mi acción de control

Modelo

Para el modelo del GPC, se utiliza una estructura de un modelo CARIMA, donde se tiene típicamente una función de transferencia junto con un termino adicional que representa algún modelo desconocido como por ejemplo una perturbación, que será una dinámica estocástica (Es algo que no sabremos como se comportará en el futuro)

\Delta A(z^{-1})y(k)=B(z^{-1})\Delta u(k-1)+C(z^{-1})\varepsilon (k)

  • A(z^{-1})=1-a_1 z^{-1}+a_2 z^{-2}+\cdots+ a_n z^{-n}
  • B(z^{-1})=b_0+b_1 z^{-1}+b_2 z^{-2}+\cdots+b_n z^{-n}
  • C(z^{-1})=c_0+c_1 z^{-1}+c_2 z^{-2}+\cdots+ c_n z^{-n}
  • y(k)=salida de la planta
  • u(k)=entrada de la planta
  • \varepsilon(k)=ruido blanco de media nula
  • \Delta=1-z^{-1}

Función de Costo

La función costo o función objetivo que pondera el error del sistema y las acciones de control del GPC esta dado por la siguiente ecuación:

J=\sum_{j=N1}^{N2}\delta (j)\left [ \hat{y}(k+j/k)-W(k+j) \right ]^2 + \sum_{j=1}^{Nu}\lambda (j)\left [ \Delta u(k+j-1) \right ]^2

Control Predictivo GPC

A continuación vamos a realizar un ejemplo de un control GPC sin Restricciones, para esto vamos a utilizar 3 funciones de transferencia distintas para observar el comportamiento y funcionamiento del control predictivo GPC.

Para este ejemplo usaremos un horizonte de predicción (N2=3) y un horizonte de control (Nu=3)

Función de transferencia 1 (Ejemplo 1)

G(z^{-1})=\dfrac{y(k)}{u(k)}=\dfrac{0.5z^{-1}}{1-0.5z^{-1}}

Función de transferencia 2 (Ejemplo 2)

G(z^{-1})=\dfrac{y(k)}{u(k)}=\dfrac{0.07056z^{-1} + 0.0621z^{-2}}{1-0.7788z^{-1}}z^{-1}

Función de transferencia 3 (Ejemplo 3)

G(z^{-1})=\dfrac{y(k)}{u(k)}=\dfrac{0.4z^{-1}+0.2z^{-2}}{1-0.6z^{-1}+0.4z^{-2}}z^{-2}

Código de Implementación:

Descargar Códigos del GPC

Código Principal

%% Control GPC
% Sergio Andres Castaño G.
% controlautomaticoeducacion.com
%___________________________________________________________________
clc
clear all 
close all

w=menu('Seleccione el ejemplo','Ejemplo 1','Ejemplo 2','Ejemplo 3');
if w==1
    %Ejemplo 1
    T=1;
    B=[0 0.5];
    A=[1 -0.5];
    %Retardo de la planta
    d=0; %Ejemplo 1
    
    %% Parametros de sintonia del GPC
    %Ventana de prediccion
    N=3;    %Ejemplo 
    N1=d+1;   %Horizonte Inicial
    N2=d+N;   %Horizonte final de salida
    Nu=N;   %Horizonte de entrada
    lambda=1;   %Parametro de ponderacion
    delta=1;
 
end

if w==2
    %Ejemplo 2
    T=2;
    B=[0 0.07056 0.0621];
    A=[1 -0.7788];
    %Retardo de la planta
    d=1; %Ejemplo 2
    %% Parametros de sintonia del GPC
    %Ventana de prediccion
    N=3;    %Ejemplo 
    N1=1;   %Horizonte Inicial
    N2=N;   %Horizonte final de salida
    Nu=N;   %Horizonte de entrada
    lambda=1;   %Parametro de ponderacion
    delta=1;
end

if w==3
    %Ejemplo 3
    T=0.5;
    B=[0 0.4 0.2];
    A=[1 -0.6 0.4];
    %Retardo de la planta
    d=2; %Ejemplo 3
    %% Parametros de sintonia del GPC
    %Ventana de prediccion
    N=3;    %Ejemplo 
    N1=d+1;   %Horizonte Inicial
    N2=d+N;   %Horizonte final de salida
    Nu=N;   %Horizonte de entrada
    lambda=1.5;   %Parametro de ponderacion
    delta=1;
end
  
ftz=filt(B,A,T);
ftz.iodelay=d

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:N1+N-1,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(N1,N)'; %Vector de controles pasados (Pertenece a la respuesta libre)
j=2;
for i=N1:N2
 aux=conv(En(i,1:end),B);
 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
 uG(j-1,1:N1)=aux(j:N1+j-1);
 j=j+1;
end



g=conv(E,B); %Calcula polinomio g

G=zeros(N,Nu);                       % Inicializa la matriz G
 
for k=1:Nu
    G(k:end,k)=g(1:N-k+1);         % Forma a matriz G                  
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=160; %Numero de interacciones
inc_u=0;
u_ant(1:10) = 0;
u(1:20+d) = 0; ym(1:20+d) = 0; r(1:20+d) = 0;

% Referência
r(d+5:40) = 1; r(41:80) = 3; r(81:120) = 2; r(120:nit) = 2;

%Perturbacion
do(1:119)= 0;do(120:nit) = 0.1;

lduf=size(F); %Averiguo la longitud del polinomio F
lduf=lduf(1,2);

luG=size(uG); %Longitud de uG
luG=luG(1,2);
duf=zeros(1,luG); %Vector incrementos de control pasados libre

for k=9+d:nit
     % Salida del proceso
     if w==1
        ym(k)=B(1)*u(k-1-d)-A(2)*ym(k-1)+do(k); %Ejemplo 1
     end
     if w==2
        ym(k)=B(1)*u(k-1-d)+B(2)*u(k-2-d)-A(2)*ym(k-1)+do(k); %Ejemplo 2
     end
     if w==3
        ym(k)=B(1)*u(k-1-d)+B(2)*u(k-2-d)-A(2)*ym(k-1)-A(3)*ym(k-2)+do(k); %Ejemplo 3
     end

     
     %Calculo de la respuesta libre
     free=0;
     for i=1:lduf
         free=free+ym(k-i+1)*F(:,i);
     end
     
     if luG==1 %Si uG es un vector columna
         free=free+duf*uG; %Invierto el producto, para que sea realizable
     else %Si uG es una matriz
         free=free+uG*duf';
     end
     %free=ym(k)*F(:,1) + ym(k-1)*F(:,2) + ym(k-2)*F(:,3) +uG*duf';
      
     
     %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
     aux_u=u_ant(1:length(B)-1);
     u_ant=[u(k) aux_u];

    %actualiza duf
    % duf= [du(k-1) du(k-2) ..... du(k-N)]

    aux_2=duf(1:end-1);
    duf=[inc_u aux_2];

end

nm=nit;
t = 0:T:(nm-1)*T;
figure
subplot(2,1,1),plot(t,r,'--k',t,ym,'-r','Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
legend('y_r','y','Location','SouthEast')
grid on;
hold
subplot(2,1,2),plot(t,u,'b','Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('u')
grid on;

Función Diofantina

% Example extracted from the book:
% Control of Dead-time Processes, Springer-Verlag, 2007
% by Julio Normey-Rico & Eduardo F. Camacho
%
%---------------------------------------------------------- %
%                                                           %
% Solution of the Diophantine equation                      %
%                                                           %
%---------------------------------------------------------- %
%[E,F] = diophantine(A,N,d)
% A = Denominator
% N = Prediction Horizon
% d = Transport delay

function[E,F] = diophantine(A,N,d)

%clear all

% Computes polynomials E(z^-1) e F(z^-1) 

% delta = 1-z^(-1)

delta = [1 -1];

% A = 1 + a1 z^(-1) + ... + a2 z(-na)

%A = [1 -0.8]

%A = [1 -0.905]

%A = [1 -1.8 0.81]

% Ã = Adelta

AD = conv(A,delta);

% note that nAD = n~a + 1

nAD = size(AD);
nAD = nAD(2);

% compute horizons

 N1 = d +1;
 N2 = d + N;

% Compute F(z^-1)

% inilialization vector f

f(1,:)= [1 zeros(1,nAD-2)];

% i = 0 ... nã-1

for j = 1: N2;

% Note that for i = 1 corresponds to f(j,0)

for i = 1:nAD-2
   f(j+1,i) = f(j,i+1)-f(j,1)*AD(i+1);
end
   f(j+1,nAD-1) = -f(j,1)*AD(nAD);
end

F = f(1+N1:1+N2,:);

% Computes E(z^-1)

E = zeros(N2);
e(1) = 1;               % for the special case  1/~A

E(1,1) = e(1);

for i = 2: N2
    e(i) = f(i,1);
    E(i,1:i)=e;
end

E = E(N1:N2,:);
Respuesta control GPC

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 Predictivo Basado en Modelo MPC

Entradas relacionadas

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.

Comentarios (23)

Buenos días Sergio , presento una duda con el programa:

la formula estándar del predictivo(GPC) es 𝑦=𝐺𝑢+𝐹 𝑧−1 𝑦 𝑡 +𝐺′ 𝑧−1 Δ𝑢 𝑡−1.

No encuentro en el programa en que variable se guardan los valores de 𝐺′ completos, en el programa encuentro que el 𝐺′ es la primera columna de este calculo pero no logro percibir donde están los otros valores, me podrías indicar donde se encuentran guardados.

Gracias.

Responder

Yesid, en el código ese 𝐺′ se llama uG.
Se calcula en la linea 82 del código.
En el último video que está en este post, se explica más o menos ese cálculo, lo puedes ver a partir del minuto 27.

Responder

Sergio saludos donde te encuentres, he implementado el GPC en un plc siemens S71200 para el control de nivel. Hasta el momento funciona pero me gustaría mejorarlo y requiero de tu ayuda.

Responder

Hola José, cual es tu duda? Saludos!

Responder

Bueno, tengo una duda con el u(k)=u(k-1)+Δu, si eso debe estar en un ciclo (FOR) o un ciclo infinito (WHILE), talvez algún momento podamos conectarnos por algún medio para compartirte mi pantalla y lo veas.

Gracias por la contestación.

Responder

el calculo de la ley de control, va en un ciclo while, pues está debe ser calculada en cada periodo de muestreo, para ser mandada a la planta. Por lo tanto siempre se debe calcular en cuanto el plc esté energizado.

Responder

Hola Sergio, tengo la función de transferencia de un inversor monofásico de segundo orden con un cero, sin retardo. Cuando simule sin retardo d=0, es decir que tome a B directamente del numerador de la función de transferencia, con esto sigue la referencia perfectamente. Sin embargo en teoría no estaría agregando el retardo del Modelo CARIMA Bz^1 y debería agregar retardo d=1 para cumplir con el modelo, cuando agrego el retardo ya no controla para nada. En el primer caso el vector de incrementos de control pasado es un vector esto es posible por el cero de la función de transferencia. Estoy poco confundido con ese tema. De igual manera muchas gracias por tus valiosos conocimientos.
Saludos.

Responder

Hola Alex, deberia ser parecido al ejemplo 2 presentado aqui. Dale un vistazo nuevamente pues corregí el código, ya que vi que no estaba usando la muestra atrasada del mocelo CARIMA, pero ya lo corregí. Saludos

Responder

Sergio, podrias decirme, que se estudia en un doctorado de controles industriales en la universidad de Brasil.?

Responder

Hola Hely, depende de lo que te quieras especializar en tu doctorado, la idea es siempre encontrar asi sea una pequeña inovación que permita avanzar en la ciencia. En mi caso yo hago un doctorado en ingeniería química, en la linea de especialización de control, modelado y simulación de procesos. En mi caso hice varias disciplinas de control clasico, control avanzado, optimización, analisis de sistemas y algunas obligatorias de la ingeniería química (las cuales me costaron algún trabajo, pues yo no soy ing. químico) Luego escoges un tema y te dedicas a investigar dicho tema por 3 años y defiendes tu tesis.

Responder

En Brasil se puede hacer un doctorado en ingeniería a distancia?

Responder

Bueno eso no lo sé, tendrías que investigar. Pero en un doctorado no necesitas estar presente para hacerlo, solo en el primer periodo cuando haces las disciplinas. Luego puedes continuarlo en tu país, si no necesitas de laboratorios o cosas asi. Eso si, tendrías que viajar algunas veces para hacer defensas de cualificacion y tesis. Y Claro depende mucho del orientador que tengas.

Responder

hola sergio, saludos desde Colombia , una pregunta necesito hacer un control predictivo a un helicóptero de 3 DOF, es un sistema no lineal, cual estrategia de control predictivo me recomienda?

Responder

Hola Yosman, todos los sistemas son No lineales, va depender mucho de la calidad del modelo que tengas de tu sistema. Para el caso del helicóptero, creo que talvez un NMPC sería lo más adecuado para controlarlo. Claro que no se cual sea el grado de no linealidad que tenga tu sistema. Si no es muy grande, puedes usar un MPC lineal, o puedes linealizar en cada punto de equilibrio. Las posibilidades son muchas. Investiga en la literatura que se ha hecho hasta el momento en ese tipo de sistemas. Éxitos.

Responder

Hola serio buenos días,estuve investigando y han manejado control NMPC, una pregunta de casualidad tiene códigos de matlab implementados sobre nmpc que me pueda facilitar, agradezco su colaboración seria de gran ayuda.
Muchas gracias

Responder

Los estoy preparando para montarlos al curso. Pronto los publicaré

Responder

Hola cual es la diferencia entre MPC y GPC

Responder

Hola. MPC es como decir la categoria. Es decir MPC (Model Predictive Control) es la categoria del controlador que utiliza algoritmos de optimización para el calculo de la ley de control y dentro de esa categoria existen varios tipos de controles predictivos como el DMC, GPC, DTC-GPC, GMV, ETC.

Responder

Hola Sergio,como podría aplicar tu explicación sobre el sistema GPC a un sistema multivariable?
Saludos y enhorabuna por tu trabajo

Responder

para un sistema MIMO se hace un poco mas largo el proceso. Basicamente las salida de predicción y=G.dU+f debes hacerlo como una matriz que contenga todas las matrices del proceso. Si por ejemplo es un proceso 2×2, vas a tener 4 matrices que representan la respuesta al escalon. asi

| G11 G12 | |dU1| |f1 |
Y= | G21 G22 | |dU2| + |f2 |

Donde G11 seria la matriz con los coeficientes de la respuesta al escalón. Y así va. Es basicamente lo mismo, solo que debes ampliar las matrices.

Responder

Hola Sergio, me gustaría saber como poder implementar el ejemplo anterior (GPC) como un sistema MIMO de 3×3.
Gracias

Responder

Estimado Sergio, con cual de los tipos de control GPC o DMC es optimo hacer un control multivariable.
Saludos

Responder

Hola Gabriel, puedes emplear cualquiera de los dos. Con el GPC seria fácil, pues puedes trabajar directo con la matriz de funciones de transferencia, Si trabajas en cambio con DMC, debes manejar una arreglo de matrices que almacenen toda la dinámica del proceso MIMO. Igual, ambos se pueden optimizar, existe una función en matlab llamada QUADPROD, que sirve justamente para optimizar, llevando en cuenta las restricciones de tu proceso. Todo esto lo pienso explicar en la pagina, solo que preparar una clase de control predictivo me lleva muchísimo tiempo. Apenas me desocupe un poco de mi tesis, continuo con los videos. Saludos.

Responder