MIMO MPC

En esta entrada veremos un poco el funcionamiento de un controlador predictivo basado en modelo (MPC) en su resolución GPC (Generalized Predictive Control) aplicado en procesos multivariables o por sus siglas en ingles MIMO (Multi-Input Multi-Output).

Antes de comenzar, si te interesa aprender más sobre la teoría del control y en especial sobre control avanzado, te invito a que mires nuestro curso gratuito de Control Predictivo y que te suscribas al canal de YouTube.

Sistema Multivariable – MIMO

Una representación basica en diagrama de bloques para un sistema MIMO 2×2 puede verse en la siguiente figura:

En este ejemplo especifico vemos un sistema que tiene 2 entradas (u1 y u2) y dos salidas (y1 y y2). Cada una de esas entradas puede afectar en cierta proporción a cada una de las salidas del sistema.

CARACTERÍSTICAS DEL MPC MIMO

• No se preocupa por realizar desacoplo de variables ya que trata el problema de forma general
• Se deben escoger adecuadamente los pares entrada-salida (Una técnica para decidir esto es la matriz RGA)
• Es importante escalar adecuadamente TODAS las variables del sistema, debido a que normalmente se trabaja con distintas variables (Temperatura, Presión, Nivel, Flujo, etc)

FUNCION OBJETIVO

La función objetivo viene siendo igual a como hemos venido trabajando en el caso SISO. Donde se tiene el seguimiento de referencia y las acciones de control, todo eso ponderado.

La diferencia es que ahora las ponderaciones no son escalares, dado a que tengo más de una variable en mi sistema. Y  para cada salida voy a tener horizontes diferentes.

La propiedad aplicada en la función objetivo es una Norma 2.

En otras palabras como para cada salida voy a tener diferentes horizontes, suponiendo a modo de ejemplo que tengo dos salidas (y1 que va desde el horizonte N11 hasta N12) y (y2 que va desde el horizonte N21 hasta N22).

Mi matriz de predicciones “y” va a contener las dos salidas, es decir que la salida estará dividida por sectores y que para cada sector voy a tener que ponderar de forma diferente, quiere decir que mi vector de ponderación también estará dividido por sectores:

Y se hace exactamente lo mismo con las ponderaciones y con el vector del incremento de control, el cual también estará dividido por sectores:

En forma general de la función objetivo podemos ver que los tamaños de cada vector o matriz viene dado por:

• Vector de predicción tamaño = Suma()
• Vector de control tiene tamaño = Suma()
• Matriz de Ponderación y bloque diagonal

Para montar la matriz de coeficientes del sistema, se deben dejar todas las entradas del sistema constantes y aplicar un escalón únicamente en una entrada y ver como se comporta cada una de las salidas del sistema. Para la primera variable manipulada podría calcular los primeros sectores de la matriz G:

Y asi sucesivamente hasta llenar toda la matriz (Este procedimiento sirve para determinar la matriz G tanto para el DMC cuanto para el GPC).

Caso GPC

Analicemos como se estructura nuestro modelo MIMO expresado en su forma CARIMA:

Donde “y” son las salidas, “P” es la función de transferencia MIMO y “u” son las entradas.

Cada función de transferencia que conforma el elemento “P” podria expresarse de la siguiente manera:

Y como es un sistema MIMO, es muy probable que cada una de las funciones de transferencia presenten retardos diferentes. Por eso lo más conveniente es tomar para cada salida del sistema MIMO (Cada fila de la matriz P) el mínimo retardo que se presente en esa salida (esa fila):

Por ejemplo si tenemos la siguiente matriz de transferencia 2×2:

Podemos extraer el retardo mínimo de cada salida y agruparlo en la diagonal principal de otra matriz, es decir llevarlo a la siguiente representación:

Donde “P” es la matriz de transferencia, “D” es la matriz con los retardos mínimos y “G” es la matriz de transferencia que contiene los retardos efectivos del sistema:

Ahora se debe aplicar el mismo procedimiento que se aplica en el caso de las perturbaciones. Si analizamos bien el problema, veremos que el sistema MIMO es un sistema que presenta interacciones entre variables y que cada interacción, nosotros podriamos considerarla como si fuera una perturbación medible.

O sea que podriamos pensar en un desacoplo por FeedForward. Entonces aplicando el mismo concepto, debemos obtener el mínimo común multiplo de cada salida como se muestra a continuación:

Con el minimo común multiplo, pasaría a calcular los nuevos polinomios “B” los cuales incluiran los retardos efectivos.

Haciendo todas estas operaciones puedo recudir mi problema a la solución de varias ecuaciones MISO (Multi-Input Single-Output) o sea multiples entradas y una sola salida, lo cual es mucho más sencillo de calcular.

Por ultimo, podemos ver que las predicciones del sistema MIMO van a partir desde la formulación del modelo multivariable:

Para poder implementar un GPC MIMO tenemos dos posibilidades equivalentes:

1. Realizar la implementación de Manera Recursiva (Como lo vimos en el primer video de GPC SISO)
2. Utilizando las Ecuaciones Diofantinas – Se repite el método aprendido, pero con cada uno de los modelos del sistema MIMO.

CALCULO DE LAS PREDICCIONES POR ECUACIONES DIOFANTINAS (Método 1)

Básicamente es similar al procedimiento hecho en el caso SISO. Para un modelo CARIMA de n salidas por m entradas puedo usar el siguiente polinomio identidad que me permite obtener la predicciones optimas de cada salida y(k+1):

La matriz sera una matriz diagonal , el polinomio tendrá un orden de  (numero de predicciones en el horizonte de predicción menos uno) y el polinomio  tendrá un odern de  es decir su tamaño depende del polinomio A relacionado a cada salida del sistema MIMO.

Retomando nuestro ejemplo con un horizonte de predicción de 4 y un horizonte de control igual a 2 Miremos las predicciones para la salida y1 comenzando desde el retardo mínimo d1

Donde las ecuaciones diofantinas E y F son:

El producto entre el polinomio E y los polinomios B para la primera salida son:

Para entender bien las predicciones voy a separar cada entrada en matrices diferentes, para ordenar los términos en las predicciones de la parte forzada y la parte libre.

CALCULO DE LAS PREDICCIONES DE FORMA RECURSIVA (Método 2)

De esa forma llegamos a que nuestra ecuación de salida para la variable 1 (y1) viene dado por la siguiente ecuación:

Y que podriamos hacer todas las predicciones hasta el horizonte de predicción N=4, note que de forma recursiva siempre en una predicción debo usar el valor que calculé en la predicción anterior, por ejemplo en la predicción de (y1(t+2)) tengo que usar el valor numérico que obtuve en la predicción anterior o sea en (y1(t+1)).

Podemos hacer el mismo procedimiento para la salida y2 y calcular todas las predicciones que queramos

Vemos que como en esta segunda salida tenemos un atraso minimo de 1, nuestra predicción comienza en d+1=2 hasta nuestro horizonte de predicción.

Ley de Control

Sabemos que el Calculo de la ley de control para el caso irrestricto viene dado por:

donde W es la referencia y f es la respuesta libre.

Sabíamos que el primer termino de esta ley de control basta con calcularlo una única vez, o sea que podríamos decir que:

llegando a la ley de control:

Y que de la matriz K solo tomabamos la primera fila, por causa del horizonte deslizante. Para un proceso Multivariable, vamos a tomar también la primera fila de cada bloque asociado a cada horizontes de control, asi:

Asi observamos que debemos tomar toda la primera fila completa de K de cada horizonte de control, para poder formar nuestra ley de control . En la imagen, habían 3 incrementos de control, por lo tanto se tomaran 3 K diferentes para poder formar las 3 leyes de control.

Recordando que cada k es un vector fila.

(Agradecimientos a Deinis Muñoz Muñoz por contribuir con las imagenes y formulas de la ley de control)

Retomando nuestro ejemplo con N=4 y Nu=2, la matriz K seria dado por:

Donde la matriz anterior podemos obtener las dos ganancias para el calculo de los incrementos de control. La ganancia 1 () seria la primera fila de la matriz. La ganancia 2 seria tomar la primera fila donde comienza el bloque para la segunda salida, y para poder descubrir cual fila corresponde debemos observar el horizonte de control anterior (de la ganancia 1 o primera salida). Sabemos que el horizonte de control de ambas salidas es  y , por lo tanto la ganancia esta en la fila:  . Entonces debo tomar toda la fila 3 para la ganancia 2 ( )

Ejemplo

Considere la siguiente columna de destilación la cual se muestra en la siguiente figura. (Este ejemplo es derivado de el modelamiento de Wood and Berry)

A continuación se muestra la función de transferencia que representa el modelo, se pide realizar con control MIMO GPC. (Discretice el modelo con Ts=0.5). Este ejemplo es el que hemos venido realizando en toda esta entrada, con todos los cálculos que ya vimos en la parte superior.

[magicactionbox id=”1184″]

A continuación se muestra el código de implementación en MATLAB, recuerda que para ver el código debes compartir el contenido de este post por redes sociales, para que más personas puedan beneficiarse de esta información.

Ejemplo Forma Diofantina

Ejemplo Forma Recursiva

[sociallocker id=948]

Para que los códigos funciones, debes bajar las siguientes librerias y colocarlas en la carpeta de MATLAB de “mis documentos” o simplemente colocarlas en la carpeta donde estas trabajando tu código MIMO GPC.

>> DESCARGAR CODIGOS (CLICK ACÁ) <<

GPC MIMO POR DIOFANTINAS:

%% MIMO GPC for Wood an Berry column
% Sergio Andres Castaño Giraldo
% https://controlautomaticoeducacion.com/mimo-mpc/
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - 2017
% Universidade Federal de Santa Catarina
% Florianópolis - 2016
% Por Ecuaciones Diofantinas
%_____________________________________________________________________

clc
close all
clear all

%% Define o processo do reator:

Ps=[tf(12.8,[16.7 1]) tf(-18.9,[21 1]);tf(6.6,[10.9 1]) tf(-19.4,[14.4 1])];
%Ps.iodelay=[1 3;7 3];
Ps.iodelay=[0 1;1 0.5];

Ts=0.5;                     %Tiempo de Muestreo
Pz=c2d(Ps,Ts);              %Processo discreto

[Bp,Ap,dp]=descompMPC(Pz);     %Descompone num, den e atraso de Pz em celdas

%% Parametros de sintonia del GPC
%Encuentro cual es el retardo minimo en mi función de transferencia
dmin=[100 100]';
for i=1:2
    dmin(i)=min(dp(i,:));
end

N=[4;4];                              %Ventana de Predicción
N1=[dmin(1)+1 dmin(2)+1];               %Horizonte Inicial
N2=[dmin(1)+N(1) dmin(2)+N(2)];         %Horizonte final
Nu=[2;2];                               %Horizonte de control
lambda=[10 100];                        %Ponderación de la acción de control
delta=[1 1];                            %Ponderación del seguimiento de referencia

%Matrices en bloque de los parametros de poderación
Ql=blkdiag(lambda(1)*eye(Nu(1)),lambda(2)*eye(Nu(2)));  %(lambda -> Acción de Control)
Qd=blkdiag(delta(1)*eye(N(1)),delta(2)*eye(N(2)));      %(delta  -> Seguimiento de Referencia)

%% Calculo de la ecuacion Diofantina
%Se obtiene la matriz B y A (Numerador y denominador sistema MIMO)
% la matriz A es el minimo común denominador del sistema
[B,A] = BA_MIMO(Bp,Ap);
%Función que Calcula Los Polinomios de la ecuación Diofantina
[E,En,F] = diophantineMIMO(A,N,dmin);

%% Matriz de la Respuesta Forzada G
[G] = MatrizG(E,B,N,Nu,dp);

%Sistema sin Restricciones
M=inv(G'*Qd*G+Ql)*G'*Qd;    %Ganancia M
K1=M(1,:);                  %Tomo la primera columna unicamente (para control u1)
K2=M(Nu(1)+1,:);            %Tomo valores de M desde [N(1)+1] una muestra despues del primer 
                            %horizonte de control (para control u2)

%Determino el Vector con los controles pasados utilizando la funcion:
uG = deltaUFree(B,En,N,dp);

%% Loop de Controle
% inicializa parametros de Simulacion
nit=400;                        %Numero de interacciones
deltaU1(1:nit) = 0;             % Inicializa controle incremental
deltaU2(1:nit) = 0;

% Señales del processo MIMO 2x2 (Inicialización)
%Salidas de cada bloque MIMO
yp11(1:nit) = 0; yp12(1:nit) = 0; yp21(1:nit) = 0; yp22(1:nit) = 0; 
%Salidas MIMO
y1(1:nit) = 0; y2(1:nit) = 0; 
%Referencias
r1(40:nit) = 0.8; r2(150:nit) = 0.4; 
%Acciones de Control
u1(1:nit) = 0; u2(1:nit) = 0; 

% Determino el tamaño de los polinomios F
lduf1=size(F{1}); %Averiguo la longitud del polinomio F
lduf2=size(F{2}); %Averiguo la longitud del polinomio F
lduf1=lduf1(1,2);
lduf2=lduf2(1,2);
F1=F{1};
F2=F{2};

%_________ Vectores de Control Pasado _____________________________
%Determino el numero de columnas de los vectores pasados
luG11=size(uG{1,1}); 
luG12=size(uG{1,2}); 
luG21=size(uG{2,1}); 
luG22=size(uG{2,2}); 

luG11=luG11(1,2);
luG12=luG12(1,2);
luG21=luG21(1,2);
luG22=luG22(1,2);

duf11=zeros(1,luG11); %Vector incrementos de control pasados libre
duf12=zeros(1,luG12); %Vector incrementos de control pasados libre
duf21=zeros(1,luG21); %Vector incrementos de control pasados libre
duf22=zeros(1,luG22); %Vector incrementos de control pasados libre

%Comienzo el lazo de Control
for k=40:nit
     % Salida del proceso Real MIMO
     yp11(k) = transferencia2(Bp{1,1},Ap{1,1},dp(1,1),u1,yp11,k);
     yp12(k) = transferencia2(Bp{1,2},Ap{1,2},dp(1,2),u2,yp12,k);
     yp21(k) = transferencia2(Bp{2,1},Ap{2,1},dp(2,1),u1,yp21,k);
     yp22(k) = transferencia2(Bp{2,2},Ap{2,2},dp(2,2),u2,yp22,k);

     y1(k)=yp11(k)+yp12(k);
     y2(k)=yp22(k)+yp21(k);

      %Puedes simular tambien el proceso de la siguiente forma:
      %t = 0:Ts:(k-1)*Ts;
      %yest=lsim(Ps,[u1(1:k);u2(1:k)],t,'zoh');
      %y1(k)=yest(end,1);
      %y2(k)=yest(end,2);
        
     %Calculo de la respuesta libre
     free1=0;
     free2=0;
     for i=1:lduf1
         free1=free1+y1(k-i+1)*F1(:,i);
     end
     for i=1:lduf2
         free2=free2+y2(k-i+1)*F2(:,i);
     end
     
     %Adiciono vectores de controles pasados
     free1=free1+uG{1,1}*duf11'+uG{1,2}*duf12';
     free2=free2+uG{2,1}*duf21'+uG{2,2}*duf22';
        
     %Calculo del incremento de Control
     %Proyeto sin referencias futuras
     deltaU1(k)=K1*[(r1(k)*ones(1,length(free1(1,1:end)))'-free1);(r2(k)*ones(1,length(free2(1,1:end)))'-free2)];
     deltaU2(k)=K2*[(r1(k)*ones(1,length(free1(1,1:end)))'-free1);(r2(k)*ones(1,length(free2(1,1:end)))'-free2)];
     
     %Calculo de la ley de control
     if k==1
        u1(k)=deltaU1(k);
        u2(k)=deltaU2(k);
     else
        u1(k)=u1(k-1)+ deltaU1(k);
        u2(k)=u2(k-1)+ deltaU2(k);
     end
     
    %Actualizo vectores de Controles Pasados 
    aux_2=duf11(1:end-1);
    duf11=[deltaU1(k) aux_2];
    
    aux_2=duf12(1:end-1);
    duf12=[deltaU2(k) aux_2];
    
    aux_2=duf21(1:end-1);
    duf21=[deltaU1(k) aux_2];
    
    aux_2=duf22(1:end-1);
    duf22=[deltaU2(k) aux_2];

end

%Grafico Resultados
nm=nit;
t = 0:Ts:(nm-1)*Ts;
figure
subplot(2,1,1),plot(t,r1,t,r2,t,y1,t,y2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
grid on;
hold
subplot(2,1,2),plot(t,u1,t,u2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('u')
grid on;

MIMO RECURSIVO

%% MIMO GPC for Wood an Berry column
% Sergio Andres Castaño Giraldo
% https://controlautomaticoeducacion.com/mimo-mpc/
% Universidade Federal de Rio de Janeiro
% Rio de Janeiro - 2017
% Universidade Federal de Santa Catarina
% Florianópolis - 2016
% Por forma Recursiva
%_____________________________________________________________________
clc
close all
clear all
%% Define o processo do reator:
Ps=[tf(12.8,[16.7 1]) tf(-18.9,[21 1]);tf(6.6,[10.9 1]) tf(-19.4,[14.4 1])];
%Ps.iodelay=[1 3;7 3];
Ps.iodelay=[0 1;1 0.5];
Ts=0.5;                     %Tiempo de Muestreo
Pz=c2d(Ps,Ts);              %Processo discreto
[Bp,Ap,dp]=descompMPC(Pz);     %Descompone num, den e atraso de Pz em celdas
%% Parametros de sintonia del GPC
%Encuentro cual es el retardo minimo en mi función de transferencia
dmin=[100 100]';
for i=1:2
dmin(i)=min(dp(i,:));
end
N=[10;10];                              %Ventana de Predicción
N1=[dmin(1)+1 dmin(2)+1];               %Horizonte Inicial
N2=[dmin(1)+N(1) dmin(2)+N(2)];         %Horizonte final
Nu=[2;2];                               %Horizonte de control
lambda=[10 100];                        %Ponderación de la acción de control
delta=[1 1];                            %Ponderación del seguimiento de referencia
%Matrices en bloque de los parametros de poderación
Ql=blkdiag(lambda(1)*eye(Nu(1)),lambda(2)*eye(Nu(2)));  %(lambda -> Acción de Control)
Qd=blkdiag(delta(1)*eye(N(1)),delta(2)*eye(N(2)));      %(delta  -> Seguimiento de Referencia)
%% Calculo de la ecuacion Diofantina
%Se obtiene la matriz B y A (Numerador y denominador sistema MIMO)
% la matriz A es el minimo común denominador del sistema
[Bt,A] = BA_MIMO(Bp,Ap);
%Calculo de las Diofantinas
[E,En,F] = diophantineMIMO(A,N,dmin);
%% Matriz de la Respuesta Forzada G
[G] = MatrizG(E,Bt,N,Nu,dp);
%Sistema sin Restricciones
M=inv(G'*Qd*G+Ql)*G'*Qd;    %Ganancia M
K1=M(1,:);                  %Tomo la primera columna unicamente (para control u1)
K2=M(Nu(1)+1,:);            %Tomo valores de M desde [N(1)+1] una muestra despues del primer 
%horizonte de control (para control u2)
%% Loop de Controle
%% inicializa parametros de Simulacion
nit=400; %Numero de interacciones
%Inicializa vector de las respuestas libres
y_free1=[];
y_free2=[];
%Inicializa vectores de los incrementos de control
deltaU1(1:nit) = 0;            
deltaU2(1:nit) = 0;
% Señales del processo MIMO 2x2 (Inicialización)
%Salidas de cada bloque MIMO
yp11(1:nit) = 0; yp12(1:nit) = 0; yp21(1:nit) = 0; yp22(1:nit) = 0; 
%Salidas MIMO
y1(1:nit) = 0; y2(1:nit) = 0; 
%Referencias
r1(40:nit) = 0.8; r2(150:nit) = 0.4; 
%Acciones de Control
u1(1:nit) = 0; u2(1:nit) = 0; 
%Atil (Producto entre polinomio A y el integrador)
At{1,1}=conv(A{1,1},[1 -1]);
At{2,2}=conv(A{2,2},[1 -1]);
%Determino el tamaño de cada polinomio Atil
na1=length(At{1,1})-1;                       
na2=length(At{2,2})-1;   
%Determino el tamaño de cada polinomio B proveniente del calculo del minimo
%comun denominador de la funcion "BA_MIMO"
nb1_11=length(Bt{1,1});                       
nb1_21=length(Bt{2,1});                       
nb1_12=length(Bt{1,2});                       
nb1_22=length(Bt{2,2});                       
%Comienzo el Lazo de Control
for k=40:nit
% Salida del proceso Real
yp11(k) = transferencia2(Bp{1,1},Ap{1,1},dp(1,1),u1,yp11,k);
yp12(k) = transferencia2(Bp{1,2},Ap{1,2},dp(1,2),u2,yp12,k);
yp21(k) = transferencia2(Bp{2,1},Ap{2,1},dp(2,1),u1,yp21,k);
yp22(k) = transferencia2(Bp{2,2},Ap{2,2},dp(2,2),u2,yp22,k);
y1(k)=yp11(k)+yp12(k);
y2(k)=yp22(k)+yp21(k);
%Puedes simular tambien el proceso de la siguiente forma: 
%t = 0:Ts:(k-1)*Ts;
%yest=lsim(Ps,[u1(1:k);u2(1:k)],t,'zoh');
%y1(k)=yest(end,1);
%y2(k)=yest(end,2);
%% Creo los vectores que conforman la RESPUESTA LIBRE
%Creo el vector de controles pasados basandome en la ecuacion de
%salida del modelo de la planta actual y(k). O sea coloco en un vector
%todos los incrementos de control de cada una de las entradas, algo como:
%[deltaU(k-1) deltaU(k-2) ... deltaU(k-nb)] donde "nb" es
%la longitud del vector B (Numerador del modelo)
%Controles pasados de y1 en el instante actual
deltaU11_pas=deltaU1(k-dp(1,1)-1:-1:k-dp(1,1)-nb1_11);           
deltaU12_pas=deltaU2(k-dp(1,2)-1:-1:k-dp(1,2)-nb1_12); 
%Controles pasados de y2 en el instante actual
deltaU21_pas=deltaU1(k-dp(2,1)-1:-1:k-dp(2,1)-nb1_21);           
deltaU22_pas=deltaU2(k-dp(2,2)-1:-1:k-dp(2,2)-nb1_22); 
%Creo un vector de salidas pasadas, similar al vector de controles
%pasados solo que en este caso con las salidas hasta la longitud del
%polinomio A (denominador del modelo)
y_pas1=y1(k:-1:k-na1+1);                         % Vetor de saída
y_pas2=y2(k:-1:k-na2+1);                         % Vetor de saída
%Ahora realizo mis predicciones hasta el horizonte de predicción N2
n=1; %Indice para la salida libre
for j=1:N2(1)
if j<=dp(1,1) 
%Comienzo a desplazar los deltaU pasados a medida que aumento
%la ventana de predicción, hasta que la predicción en este
%vector llegue al retardo
deltaU11_pas=[deltaU1(k-dp(1,1)-1+j) deltaU11_pas(1:end-1)];
elseif j>dp(1,1) 
%Una vez la ventana de predicción es mayor que el retardo, en
%la ecuación de salida dejan de existir los deltaU pasados y
%comenzaran a aparecer los deltaU futuros, por eso comienzo a
%cerar los deltaU pasados.
deltaU11_pas=[0 deltaU11_pas(1:end-1)];
end
if j<=dp(1,2) 
deltaU12_pas=[deltaU2(k-dp(1,2)-1+j) deltaU12_pas(1:end-1)];
elseif j>dp(1,2) 
deltaU12_pas=[0 deltaU12_pas(1:end-1)];
end
%Hago el calculo de la respuesta libre con el polinomio Atil y el
%polinomio B, junto con las salidas pasadas y los deltaU pasados
y_free1(n)=-y_pas1*At{1,1}(2:end)'+deltaU11_pas*Bt{1,1}(1:end)'...
+deltaU12_pas*Bt{1,2}(1:end)';% Armazena yfr
%Actualizo el vector de salidas pasadas
y_pas1=[y_free1(n) y_pas1(1:end-1)];                    % Atualiza yfr
n=n+1;
end
%Ahora realizo mis predicciones hasta el horizonte de predicción
n=1; %Indice para la salida libre
for j=1:N2(2)
if j<=dp(2,1)
deltaU21_pas=[deltaU1(k-dp(2,1)-1+j) deltaU21_pas(1:end-1)];
elseif j>dp(2,1)
deltaU21_pas=[0 deltaU21_pas(1:end-1)];
end
if j<=dp(2,2)
deltaU22_pas=[deltaU2(k-dp(2,2)-1+j) deltaU22_pas(1:end-1)];
elseif j>dp(2,2)
deltaU22_pas=[0 deltaU22_pas(1:end-1)];
end
y_free2(n)=-y_pas2*At{2,2}(2:end)'+deltaU21_pas*Bt{2,1}(1:end)'...
+deltaU22_pas*Bt{2,2}(1:end)';% Armazena yfr
%Actualizo el vector de salidas pasadas
y_pas2=[y_free2(n) y_pas2(1:end-1)];                    % Atualiza yfr
n=n+1;
end
%% Calculo del Incremento de Control   
%Calculo los incrementos de control utilizando unicamente los datos de
%las respuestas libres que se encuentran dentro de los horizontes de
%predicción N1 y N2 y hago el producto con la ganancia calculada de la
%respuesta forzada
deltaU1(k)=K1*[(r1(k)*ones(1,length(y_free1(N1(1):N2(1))))'-y_free1(N1(1):N2(1))');(r2(k)*ones(1,length(y_free2(N1(2):N2(2))))'-y_free2(N1(2):N2(2))')];
deltaU2(k)=K2*[(r1(k)*ones(1,length(y_free1(N1(1):N2(1))))'-y_free1(N1(1):N2(1))');(r2(k)*ones(1,length(y_free2(N1(2):N2(2))))'-y_free2(N1(2):N2(2))')];
%% Calculo la acción de Control
if k==1
u1(k)=deltaU1(k);
u2(k)=deltaU2(k);
else
u1(k)=u1(k-1)+ deltaU1(k);
u2(k)=u2(k-1)+ deltaU2(k);
end
end
%% Grafico los Resultados
nm=nit;
t = 0:Ts:(nm-1)*Ts;
figure
subplot(2,1,1),plot(t,r1,t,r2,t,y1,t,y2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Saida');
grid on;
hold
subplot(2,1,2),plot(t,u1,t,u2,'Linewidth',2)
xlabel('Tempo (s)');
ylabel('Controle');
legend('u')
grid on;

[/sociallocker]