Saltar al contenido
Control Predictivo

Control Predictivo GPC

En este post, profundizaremos en la teoría del Control Predictivo Generalizado (GPC), un tema esencial para quienes están interesados en el control avanzado de procesos. Comenzaremos con una breve explicación sobre qué es el GPC y por qué es relevante, seguido de una serie de videos que ilustrarán detalladamente su funcionamiento.

Si eres nuevo en el tema, te recomiendo comenzar viendo nuestros videos sobre el Control DMC (Dynamic Matrix Control), ya que proporcionan fundamentos cruciales para entender el GPC. Haz clic aquí para ver los videos del DMC.

Contenidos mostrar

¿Qué es el Control GPC?

El Control Predictivo Generalizado, conocido por sus siglas GPC, es una técnica de control avanzada que utiliza modelos matemáticos para predecir y optimizar la respuesta del sistema en futuras etapas de control. A diferencia del DMC, donde el modelo se basa en los coeficientes de la respuesta al escalón, el GPC utiliza una Función de Transferencia Discreta para anticipar las dinámicas futuras.

¿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

El GPC fue inicialmente propuesto por D. W. Clarke, C. Mohtadi y P. S. Tuffs en dos artículos fundamentales:

Estos documentos sentaron las bases teóricas del GPC y sus extensiones.

Objetivo del GPC

El principal objetivo del GPC es determinar una secuencia de control óptima que permita alcanzar un objetivo específico mientras minimiza la variación en las acciones 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
  • N1 y N2 son los límites del horizonte de predicción.
  • Nu es el horizonte de control.
  • δ(j) y λ(j) son pesos que equilibran la importancia entre el error de predicción y el esfuerzo de control.
  • W(k+j) representa el valor deseado o referencia en el tiempo futuro (k+j).

Importancia de los Pesos en la Función de Costo

Los pesos 𝛿(𝑗)δ(j) y 𝜆(𝑗)λ(j) en la función de costo son cruciales porque permiten al diseñador del control ajustar la respuesta del controlador. Un mayor peso en 𝛿(𝑗)δ(j) enfatiza la importancia de seguir la trayectoria deseada, mientras que un mayor 𝜆(𝑗)λ(j) pone énfasis en minimizar el cambio en la señal de control, lo que es particularmente útil para evitar acciones de control excesivas que podrían desestabilizar el sistema o acelerar el desgaste de los actuadores.

Ecuaciones Diofantinas en el Control GPC

Las ecuaciones diofantinas juegan un papel clave en el cálculo de los parámetros de control en el GPC. Estas ecuaciones permiten resolver el problema de control de manera que se satisfaga la ecuación del modelo del sistema junto con la minimización de la función de costo.

Aplicación de las Ecuaciones Diofantinas

En el contexto del GPC, las ecuaciones diofantinas ayudan a determinar los valores futuros de la señal de control 𝑢(𝑘) que minimizarán la función de costo dada. La solución típica involucra:

  1. Establecer la relación deseada entre la salida futura y la señal de control usando una función de transferencia deseada.
  2. Resolver la ecuación diofantina para encontrar los polinomios que describen cómo las entradas afectarán las salidas dentro del horizonte de predicción.

Ecuación Diofantina Típica en GPC

Una forma de ecuación diofantina que se puede utilizar en el GPC es:

A(z ^{−1} )Y(z)=B(z ^{−1} )U(z)+C(z ^{−1} )E(z)

Donde:

  • Y(z), U(z), E(z) son las transformadas Z de la salida, entrada, y error de predicción, respectivamente.
  • A, B, C son polinomios que representan el comportamiento del sistema y la estructura del controlador.

Esta ecuación se resuelve para 𝑈(𝑧), proporcionando la ley de control que minimiza la función de costo sobre el horizonte de predicción y control establecidos.

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=5;   %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=2;   %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=3;   %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,N,d); %Calculo de la funcion Diofantina
E=En(end,:); %Polinomio E seria la ultima fila arrojada por la funcion


%% 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-N1+1,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
     

     free=free+uG*duf';
     %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

    %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,:);
GPC Control MPC

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 (25)

Buenas noches desde España, Sergio

A la hora de ejecutar las simulaciones 2 y 3 los resultados son incorrectos porque a la hora de calcular la diofantina, se está llamando a la función con el parametro delay=0, por lo que los coeficientes obtenidos son incorrectos. Me he dado cuenta reimplementando el código en Python.

Muchisimas gracias por el contenido que nos ofreces, me está siendo de inestimable ayuda! un saludo

Responder

Lo he revizado y he montado algunas correcciones que hice hace tiempo solo que olvidé montarlas aquí. Éxitos Guzmán.

Responder

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