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.
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:
- Generalized Predictive Control Algorithm Part I, the Basic
- 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)
=salida de la planta
=entrada de la planta
=ruido blanco de media nula
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 Basado en Modelo – DMC

Control NMPC – Matlab – Simulink

MPC advanced control applied to a gas compression system of an offshore platform
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)
Función de transferencia 2 (Ejemplo 2)
Función de transferencia 3 (Ejemplo 3)
Código de Implementación:
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,:);
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.