Restricciones en el MPC

En esta entrada comenzaremos a entender una de las potencialidades que tiene un control MPC con relación a un control clásico, y es el poder involucar las restricciones del proceso en el calculo de la ley de control.

Antes de comenzar, te hago la inviración a que veas el CURSO GRATIS DE CONTROL MPC

Constrained MPC – Restricciones

Como habíamos mencionado en vídeos pasados, los procesos de una industria presentan un sin numero de restricciones, los cuales debemos involucrar en nuestro proyecto de control predictivo.

Video aula en Youtube

Función de Costo

Recordando la función de costo vista en el video del  GPC, debíamos minimizarla usando la función de quadprod de Matlab:

Sin Restricción

\nabla J = 0 \to \Delta u = K(w – f)

con

K ={({G^T}{Q_\delta }G + {Q_\lambda })^{ – 1}}{G^T}{Q_\delta }

(Primera fila de la matiz)

Con Restricción

Para usar un problema con restricciones debemos usar un algoritmo de programación cuadrática que considere las restricciones dentro de mi función de optimización, para eso debemos usar la siguiente formulación patrón, y así toda restricción que deseemos involucrar en la optimización, debemos expresarla de las formas que establece dicha función patrón siempre en función del incremento de control. Las restricciones pueden ser de desigualdad o de igualdad.

J =\dfrac{1}{2}\Delta {u^T}H\Delta u + {f_0}^T\Delta u s.a\ A\Delta u \le b\\{A_{eq}}\Delta u = {b_{eq}}

con

H = {G^T}{Q_\delta }G + {Q_\lambda } {f_0} = {(f – w)^T}{Q_\delta }G

Comencemos viendo diferentes tipos de restricciones que podemos involucrar en el MPC.

Entradas relacionadas con MPC

Restricción en el incremento de la acción de control

Cualquier técnica MPC considera las restricciones en la ley de control.

La restricción del incremento de control esta definida como:

Esa restricción debemos expresarla en la forma de desigualdad que se mostró en la ecuación patrón osea:

Tomemos el primer intervalo:

Representandolo de forma matricial tenemos:

Entonces tenemos que para el incremento de control máximo esta dado por:

Por otro lado tomando la restricción del incremento de control mínimo:

en forma matricial tendríamos que

en forma reducida podemos expresarlo como:

Como el sentido de la desigualdad para el caso del incremento de control mínimo no es el mismo que el de la formula patrón, multiplicamos por menos uno en ambos lados de la desigualdad para cambiarle el sentido a la misma:

Asi conseguimos llevarla a la misma representación. Por último vamos a mezclar los dos extremos de cada restricción en una restricción única.

Creo una matriz con un bloque de matrices identidades, multiplicado por menor igual un bloque de vectores de Umaximos sobre un bloque de vectores de Uminimos asi:

Y listo conseguimos representar la restricción de desigualdad del incremento de control.

Restricción en la señal de control

Podemos representar en forma matricial la ley de control

De forma rápida, podemos representar la matriz triangular inferior que acompaña el incremento de control como T, asi la restricción está representada por:

Hasta aca, las restricciones vistas son definidas por el propio controlador. Osea el mismo controlador sabe cuales van a ser sus limitaciones para aplicar sobre el proceso. Por otro lado si pensamos en las restricciones de la salida de la planta, esa restricción será basada en el modelo.

Entradas de Control Clásico

Restricciones en la salida

Como las restricciones son hechas en base al modelo, debemos asegurarnos de tener un muy buen modelo para garantizar que las restricciones de salida sean cumplidas.

La restricción de salida debe representarse en función de ∆u

sabemos que la salida se representa por:

Asi la restricción viene dada por:

Restricción de sobre impulso

Esta restricción típicamente actúa únicamente cuando hay cambios de referencia.

Si la referencia es positiva, siempre la variable debe ser menor que la referencia, si la referencia es negativa, la variable será mayor que la referencia. O sea que la restricción siempre va a depender del sentido de la referencia.

Generalmente se define un intervalo N01 y N02. Pues normalmente la restricción debe ser aplicado cuando la variable se acerca  a la referencia, y no es necesario usar mi algoritmo de búsqueda desde las primeros periodos de muestreo. Si no supiera definir esos horizontes, puedo tomar todo el horizonte completo, la única diferencia es que mi algoritmo va a tener que calcular todo el horizonte completo sin necesidad.

Comportamiento Monotonico

Se fuerza a que la respuesta tenga un comportamiento estrictamente creciente o estrictamente decreciente.

Para el caso estrictamente creciente se tiene que:

El y(k) seria la salida actual y y’ seria el vector y(k)  sin⁡  la ultima fila.

Restricción Terminal

• En esta restricción se especifica que al final del horizonte se tenga un comportamiento especifico, por ejemplo:
• y(k)=w(k)
• Su ventaja, es que puedo garantizar que mi variable va a seguir la referencia
• Desventaja, puede ser que no sea factible. Por ejemplo que sea físicamente imposible alcanzar la referencia. Tambien, por ejemplo si tuviera un horizonte pequeño digamos 4 muestras, no voy a poder garantizar que en esas 4 muestras voy a alcanzar la referencia. Por otro lado si tengo un horizonte de control pequeño, el va a tener una variación de la acción de control muy brusca

Esta es una restricción de IGUALDAD, entonces solo tomamos la ultima fila de  las matrices, ya que eso son los datos del final del horizonte, llamamos estas filas como , , .

Código

A continuación se muestra un código ilustrativo de como programar las restricciones en MATLAB

Código de Implementación:

Para ver el código de implementación debes compartir el contenido de este post, de esa forma ayudaras con que este sitio web pueda seguir continuando en pie brindando mucho más contenido gratuito.

[sociallocker id=948]

%% Control GPC
% Sergio Andres Castaño G.
% controlautomaticoeducacion.com
%___________________________________________________________________
clc
clear all 
close all
rest=menu('Parametros de Simulacion','Restriccion de aumento de control',...
    'Restriccion en señal de control','Restriccion en la Salida',...
    'Restriccion de Sobre Impulso','Restriccion Monotonica','Sin Restricciones');
gps=tf(50,[1 0 25]);
T=0.1;
d=0;
ftz=c2d(gps,T);
[B,A]=tfdata(ftz,'v');
ftz=filt(B,A,T)


%ftz.iodelay=d
%% Parametros de sintonia del GPC
%Ventana de prediccion
N=11;    %Ejemplo 
N1=d+1;   %Horizonte Inicial
N2=d+N;   %Horizonte final de salida
Nu=N;   %Horizonte de entrada
lambda=50;   %Parametro de ponderacion
delta=1;

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
G0=zeros(N,Nu);
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=300; %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+15:80) = 1; r(81:140) = 4; r(141:200) = 2; r(201:nit) = 3;
do(1:nit)=0;
%Perturbacion
%do(1:259)= 0;do(260: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
r1=1;
for k=9+d:nit
     % Salida del proceso
       ym(k)=B(1)*u(k-1-d)+B(2)*u(k-2-d)...
            -A(2)*ym(k-1)-A(3)*ym(k-2);
 
     
     %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';
     
     switch rest
         case 1
             %----------- Restriccion del incremento de Control -----------------%
             I=eye(Nu); %Matriz identidad 
             InUmax=0.05*ones(Nu,1); %Incremento maximo
             InUmin=-0.05*ones(Nu,1); %Incremento maximo

             a=[I; -I];
             b=[InUmax; -InUmin];
             H=(G'*Qd*G+Ql);
             Fo=(free-r(k))'*Qd*G;
             %Calculo do Controle
             %options = optimset('LargeScale','off');
             [x,fval,exitflag] = quadprog(H,Fo,a,b,[],[],[],[],[]);
             inc_u=x(1);
      
         case 2 
             %----------- Restriccion en la Señal de Control -----------------%
             Triang=tril(ones(Nu));  
             T=ones(Nu,1);
             ub=1.5; %Maxima señal de control
             u_max=ub*ones(Nu,1);  
             u_min=-0*ones(Nu,1);
             a=[Triang; -Triang];
             b=[u_max-T*u(k-1); T*u(k-1)-u_min];
             
             H=(G'*Qd*G+Ql);
             Fo=(free-r(k))'*Qd*G;
             %Calculo do Controle
             options = optimset('LargeScale','off');
             [x,fval,exitflag] = quadprog(H,Fo,a,b,[],[],[],[],[],options);
             inc_u=x(1);
         case 3 
             %----------- Restriccion en la Salida -----------------%
             yb=3.2; %Maxima señal de salida
             y_max=yb*ones(N,1);  
             y_min=-0*ones(N,1);
             a=[G; -G];
             b=[y_max-free; free-y_min];
             
             H=(G'*Qd*G+Ql);
             Fo=(free-r(k))'*Qd*G;
             %Calculo do Controle
             options = optimset('LargeScale','off');
             [x,fval,exitflag] = quadprog(H,Fo,a,b,[],[],[],[],[],options);
             inc_u=x(1);
         case 4 
             %----------- Restriccion Sobre Impulso -----------------%
             if r(k)>r(k-1)
                 r1=1;
             end
             if r(k)<r(k-1)
                 r1=-1;
             end
             if r(k)==r(k-1)
                 r1=r1;
             end
             refP=(r(k))*ones(N,1);  %Referencia 
             %refN=-(r(k))*ones(Nu,1);  %Referencia Negativa

             a=r1*[G];% -G];
             b=[r1*refP-r1*free];%; free-refN];
             
             H=(G'*Qd*G+Ql);
             Fo=(free-r(k))'*Qd*G;
             %Calculo do Controle
             options = optimset('LargeScale','off');
             [x,fval,exitflag] = quadprog(H,Fo,a,b,[],[],[],[],[],options);
             inc_u=x(1);
             
          case 5   
              %----------- Restriccion Monotonica -----------------% 
             for i=N:-1:2
                G0(i,1:end)=G(i-1,1:end);                     
             end
             if r(k)>r(k-1)
                 r1=-1;
             end
             if r(k)<r(k-1)
                 r1=1;
             end
             if r(k)==r(k-1)
                 r1=r1;
             end
             
             a=r1*(G-G0);
             b=r1*[ym(k)-free(1);free(1:end-1)-free(2:end)];
             
             
             H=(G'*Qd*G+Ql);
             Fo=(free-r(k))'*Qd*G;
             %Calculo do Controle
             options = optimset('LargeScale','off');
             [x,fval,exitflag] = quadprog(H,Fo,a,b,[],[],[],[],[],options);
             inc_u=x(1);
         case 6
          inc_u=K1*(r(k)-free);
     end
     
     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)
stairs(t,r,'--k','Linewidth',2),hold on
stairs(t,ym,'-r','Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
legend('y_r','y','Location','SouthEast')
grid on;
hold
subplot(2,1,2)
stairs(t,u,'b','Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('u')
grid on;

[/sociallocker]

Si te gusto la información no olvides en compartir en redes sociales, también te puedes sucribir a mi canal de youtube o darle me gusta a mi pagina en Facebook.

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.