MPC por Variables de Estado

En esta entrada haremos una introducción al CONTROL MPC calculado por VARIABLES de ESTADO, así que estante atento y mira el código de ejemplo que vemos al final de este post.

Antes de Comenzar, podría interesarte ver el CURSO GRATUITO DE MPC

Contenidos mostrar

MPC por Variables de Estado

Antes que nada, será importante de que entiendas que son y como se representan las variables de estado, y esa información ya esta disponible en este sitio web.

Se comienza describiendo una planta SISO descrita por:

x_m(k+1)=A_mx_m(k)+B_mu(k)

y(k)=C_mx_m(k)

Donde $u(k)$ es la entrada o variable manipulada, $y(k)$ es la salida del proceso y $x_m(k)$ es el vector de estados de dimensión $n_1$ .

Debemos adaptar nuestro modelo en espacio de estados para que cumpla el proposito del diseño en el que se incluye un integrador.

Tenga en cuenta que una formulación general de un modelo de espacio de estado tiene un término directo de la señal de entrada $u(k)$ a la salida $y(k)$ como:

y(k)=C_mx_m(k)+D_mu(k)

Como siempre asumimos que la entrada $u(k)$ no afecta la salida $y(k)$ en el mismo tiempo por el principio de horizonte deslizante. Con esto se dice que $D_m=0$ .

Aplicando una operación de diferencia en ambos lados de la ecuacion, tenemos que:

x_m(k+1)-x(k)=A_m(x_m(k)-x_m(k-1))+B_m(u(k)-u(k-1)) (1)

Podemos expresar esas diferencias en terminos de incrementos:

\Delta x_m(k+1)=A_m\Delta x_m(k)+B_m\Delta u(k)

Note que como ahora la entrada en espacio de estados es $\Delta u(k)$ debemos relacionar $\Delta x_m(k)$ con la salida $y(k)$ . Esto es fácil creando un nuevo vector de estados:

x(k)=\begin{bmatrix} \Delta x_m(k)^T\\ y(k) \end{bmatrix}

De esta forma la salida del sistema es descrita por:

y(k+1)-y(k)=C_m(x_m(k+1)-x_m(k))=C_m\Delta x_m(k+1)

y(k+1)=C_m\Delta x_m(k+1)+y(k)

Reemplazando $\Delta x_m(k)$

y(k+1)=C_mA_m\Delta x_m(k)+C_mB_m\Delta u(k) + y(k)\ \ \ (2)

Colocando juntos la ecuación (1) y (2) tenemos la siguiente representación en espacio de estados:

\overset{x(k+1)}{\overbrace{\begin{bmatrix} \Delta x_m(k+1)\\ y(k+1) \end{bmatrix}}}=\overset{A}{\overbrace{\begin{bmatrix} Am & 0_m^T\\ C_mA_m & 1 \end{bmatrix}}}\overset{x(k)}{\overbrace{\begin{bmatrix} \Delta x_m(k)\\ y(k) \end{bmatrix}}}+\overset{B}{\overbrace{\begin{bmatrix} B_m\\ C_mB_m \end{bmatrix}}}\Delta u(k)

y(k)=\overset{C}{\overbrace{\begin{bmatrix} 0_m & 1 \end{bmatrix}}}\begin{bmatrix} \Delta x_m(k)\\ y(k) \end{bmatrix}\ \ \ \ \ (3)

Donde $0_m=\begin{bmatrix}0 & 0 & 0 & \cdots & 0\end{bmatrix}$ de dimensión $n_1$ y las matrices A,B,C son conocidas como las matrices aumentadas, las cuales son empleadas para el diseño del control predictivo.

Veamos entonces como sucede la predicción de los estados:

x(k+1|k)=Ax(k)+B\Delta u(k)

x(k+2|k)=Ax(k+1|k)+B\Delta u(k+1|k)=A^2x(k)+AB\Delta u(k)+B\Delta u(k+1|k)

\vdots

x(k+Np|k)=A^{N_p}x(k)+A^{N_p-1}B\Delta u(k)+A^{N_p-2}B\Delta u(k+1)+…+A^{N_p-N_c}B\Delta u(k+N_c-1|k)

Seguidamente veremos como es que sucede la predición de salida:

y(k+1|k)=Cx(k+1)=CAx(k)+CB\Delta u(k)

y(k+2|k)=CAx(k+1|k)+CB\Delta u(k+1|k)=CA^2x(k)+CAB\Delta u(k)+CB\Delta u(k+1|k)

\vdots

y(k+Np|k)=CA^{N_p}x(k)+CA^{N_p-1}B\Delta u(k)+CA^{N_p-2}B\Delta u(k+1)+…+CA^{N_p-N_c}B\Delta u(k+N_c-1|k)

Todas las predicciones son formuladas en terminos de las variables de estado actuales $x(k)$ y del futuro en el movimiento de control $\Delta u(k+j)$ donde $(j=1,2,3,...,N_c-1)$ . Asi se definen los vectores como:

Y=\begin{bmatrix} y(k+1) &y(k+2) &\cdots & y(k+N_p) \end{bmatrix}^T

\Delta U=\begin{bmatrix} \Delta u(k) & \Delta u(k+1) & \cdots & \Delta u(k+N_c-1) \end{bmatrix}^T

Para un caso SISO la dimensión de Y será $N_p$ y la dimensión de $\Delta U$ será $N_c$ .

De forma compacta, las predicciones de salida las podemos represenar como:

$Y=Fx(k)+G\Delta u$

Donde:

$F=\begin{bmatrix} CA\\ CA^2\\ CA^3\\ \vdots \\ CA^{N_p} \end{bmatrix}$

$G=\begin{bmatrix} CB & 0 & 0 & \cdots & 0\\ CAB& CB & 0 & \cdots & 0\\ CA^2B& CAB & CB & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots &\vdots \\ CA^{N_p-1}B& CA^{N_p-2}B & CA^{N_p-3}B & \cdots & 0 \end{bmatrix}$

En este punto te podría interesar estos temas:

GPC con Perturbaciones

Restricciones en el MPC

MIMO MPC

FUNCIÓN OBJETIVO

Tenemos la siguiente función objetivo que deseamos minimizar. La función es similar a la vista en el DMC y GPC:

J=(R-y)^TQ_{\delta}(R-y)+\Delta u^TQ_{\lambda}\Delta u

Donde $R^T=[1,1,1,...,1]w(k)$ es el vector que contiene las referencias, $w(k)$ es la referencia actual, $Q_{\delta}$ es la matriz de ponderación del error de seguimiento y $Q_{\lambda}$ es la matriz de ponderación del incremento de control.

Para el caso sin restricción solo basta con cumplir con:

\dfrac{\partial J}{\partial \Delta u}=0

Realizando el mismo procedimiento empleado en la entrada del DMC y GPC, el incremento de control para el caso sin restricciones se reduce a:

\Delta U=(G^TQ_{\delta}G+Q_{\lambda})^{-1}G^TQ_{\delta}(R-Fx(k))

donde el termino $(G^TQ_{\delta}G+Q_{\lambda})^{-1}G^TQ_{\delta}$ puedo calcularlo una sola vez y reducirlo a una ganancia:

$K_1=(G^TQ_{\delta}G+Q_{\lambda})^{-1}G^TQ_{\delta}$

Recordando que solo me interesa la primera fila de esta matriz.

Asi mi incremento de control que debe ser calculado en cada interacción es:

$\Delta U=K_1(R-Fx(k))$

Caso con Restricciones

Para el caso restricto partimos de la función de costo, y aplicamos el mismo concepto visto en la clase de las restricciones del MPC.

$J=(R-y)^TQ_{\delta}(R-y)+\Delta u^TQ_{\lambda}\Delta u$

Sabemos que $Y=Fx(k)+G\Delta u$ , entonces reemplazamos

$J=(R-Fx(k)-G\Delta u)^TQ_{\delta}(R-Fx(k)-G\Delta u)+\Delta u^TQ_{\lambda}\Delta u$

Desarrollando tenemos que:

J=[(R-Fx(k))^T-\Delta u^TG^T]Q_{\delta}[(R-Fx(k))-G\Delta u]+\Delta u^TQ_{\lambda}\Delta u

J=(R-Fx(k))^TQ_{\delta}(R-Fx(k))-2(R-Fx(k))^TQ_{\delta}G\Delta u+\Delta u^TG^TQ_{\delta}G\Delta u+\Delta u^TQ_{\lambda}\Delta u

J=\underset{ParteQueOptimizo}{\underbrace{\Delta u^T(G^TQ_{\delta}G+Q_{\lambda})\Delta u+2(Fx(k)-R)^TQ_{\delta}G\Delta u}}+\underset{NoInfluenciaEnLaOptimizacion}{\underbrace{(R-Fx(k))^TQ_{\delta}(R-Fx(k))}}

Multiplico por (1/2)

J=\dfrac{1}{2}\Delta u^T\underset{H}{\underbrace{(G^TQ_{\delta}G+Q_{\lambda})}}\Delta u+\underset{F_o^T}{\underbrace{(Fx(k)-R)^TQ_{\delta}G}}\Delta u

Código de Implementación (MATLAB code)

A continuación se muestra un ejemplo en MATLAB de como implementar el MPC en espacio de estados, para poder ver el código debes compartir el contenido de este post y así ayudar a que la pagina siga creciendo. Gracias por contribuir con el sitio web.

[sociallocker id=948]

Diseñar un control predictivo por variables de estado para el siguiente sistema:

$x_m(k+1)= \begin{bmatrix} 1 & 1\\ 0 & 1 \end{bmatrix} x_m(k)+\begin{bmatrix} 0.5\\ 1 \end{bmatrix}u(k)$

$y(k)=\begin{bmatrix} 1 & 1 \end{bmatrix}x_m(k)$

% Sergio Andres Castaño Giraldo
% http:// controlautomaticoeducacion.com
% ____________________________________________________________________
%
clc
clear all
close all

p=menu('EEMPC','Con Restricciones','Sin Restricciones');

%Parametros de Sintonia del MPC
Np=20;              %Horizonte de prediccion
Nc=4;               %Horizonte de Control
lambda=1;           %Parametro de ponderacion
delta=1;            %Parametro de ponderacion
Ql=eye(Nc)*lambda;  %Matriz de Ponderacion
Qd=eye(Np)*delta;

%Variables de Estado caso Integrador
%X(k+1)=AX(k)+Bu(k)

Am=[1 1;0 1];
Bm=[0.5;1];
Cm=[1 1];

Na=size(Am);
Om=zeros(1,Na(1)); %Numero de zeros para completar la matriz aumentada

%Espacio de estados
A=[Am Om';Cm*Am 1];
B=[Bm;Cm*Bm];
C=[Om 1];

%Calculo de la matriz G y vector de respuesta libre F
h(1,:)=C;       %Auxiliar para calcular G
F(1,:)=C*A;     %Calculo del vetor F
for kk=2:Np
    h(kk,:)=h(kk-1,:)*A;  %Voy calculado C*A aumentando el indice de A
    F(kk,:)= F(kk-1,:)*A; %Voy aumentando el indice de A
end
v=h*B;              %Multiplico por B para obter la primera columna de G
G=zeros(Np,Nc);   %Creo la Matriz G
G(:,1)=v;         %Coloco los datos en la matriz G
for i=2:Nc
    G(:,i)=[zeros(i-1,1);v(1:Np-i+1,1)]; %Expando em toda la Matriz G
end

%% Loop de Controle

%% inicializa parametros de Simulacion
nit=160;    %Numero de interacciones
inc_u=0;
u_ant(1:10) = 0;
u(1:20) = 0; ym(1:20) = 0; r(1:20) = 0;

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

[n,n_in]=size(B);       %Salvo a dimenção de B
xm=[0;0];               %Estados iniciales del Processo
Xf=zeros(n,1);          %Realimentacion de estados
y=0;

for k=2:nit;
    ym(k)=y;
    %----------- Restriccion en la Señal de Control -----------------%
    Triang=tril(ones(Nc));  
    T=ones(Nc,1);
    u_max=0.3*ones(Nc,1);  
    u_min=-0.3*ones(Nc,1);
    %----------- Restriccion en la salida (estados) -----------------%
    x_max=5*ones(Np,1);  
    x_min=-5*ones(Np,1);
    
    
    a=[Triang; -Triang;G;-G];
    b=[u_max-T*u(k-1); T*u(k-1)-u_min;x_max-F*Xf; F*Xf-x_min];
    
    H=(G'*Qd*G+Ql);
    Fo=(F*Xf-r(k)*ones(Np,1))'*Qd*G;
    
    %Calculo del Control
    
    if p==1
    %---------------------------Con restriccion --------------------------
        options = optimset('LargeScale','off');
        [x,fval,exitflag] = quadprog(H,Fo,a,b,[],[],[],[],[],options);
        inc_u=x(1);
    end
    if p==2
    %---------------------------Sin restriccion --------------------------
        Mn=inv(H)*G'*Qd;
        K1=Mn(1,:);            %Tomo a primeira linha del incremento de control
        inc_u=K1*(r(k)-F*Xf);  %Incremento de control    
    end
    %-----------------------------------------------------------------------------------
    
    u(k)=u(k-1)+inc_u;            %Calculo de la ley de control
    xm_old=xm;                      %Guardo los estados anteriores
    xm=Am*xm+Bm*u(k);                  %Calculo los estados actuales
    y=Cm*xm;                   %Salida atual
    Xf=[xm-xm_old;y];              %Calculo los estados realimentados
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]

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

Sergio A. Castaño Giraldo

Mi nombre es Sergio Andres Castaño Giraldo, y en este sitio web voy a compartir una de las cosas que mas me gusta en la vida y es sobre la Ingeniería de Control y Automatización. El sitio web estará en constante crecimiento, voy a ir publicando material sobre el asunto desde temas básicos hasta temas un poco más complejos. Suscríbete al sitio web, dale me gusta a la página en Facebook y únete al canal de youtube. Espero de corazón que la información que comparto en este sitio, te pueda ser de utilidad. Y nuevamente te doy las gracias y la bienvenida a control automático educación.

MPC por Variables de Estado

FUNCIÓN OBJETIVO

Caso con Restricciones

Código de Implementación (MATLAB code)

Comparte esto: