MPC en Variables de Estado

Sergio A. Castaño Giraldo

hace 10 años

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

1. MPC por Variables de Estado

1.1. Control Predictivo Basado en Modelo – DMC

1.2. NMPC Usando Modelos de Hammerstein

1.3. MPC Toolbox usando Matlab y Simulink

2. FUNCIÓN OBJETIVO

2.1. Caso con Restricciones

3. Código de Implementación (MATLAB code)

4. Citar este Artículo

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 es la entrada o variable manipulada, es la salida del proceso y es el vector de estados de dimensión .

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 a la salida como:

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

Como siempre asumimos que la entrada no afecta la salida en el mismo tiempo por el principio de horizonte deslizante. Con esto se dice que .

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 debemos relacionar con la salida . 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

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 de dimensión 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 y del futuro en el movimiento de control donde . 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á y la dimensión de será .

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

Donde:

En este punto te podría interesar estos temas:

Control Predictivo Basado en Modelo – DMC

NMPC Usando Modelos de Hammerstein

MPC Toolbox usando Matlab y Simulink

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 es el vector que contiene las referencias, es la referencia actual, es la matriz de ponderación del error de seguimiento y 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 puedo calcularlo una sola vez y reducirlo a una ganancia:

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:

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.

Sabemos que , entonces reemplazamos

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:

% 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

Citar este Artículo

📌 Selecciona el formato:

Giraldo, S. A. C. (2016). MPC en Variables de Estado. Control Automático Educación. https://controlautomaticoeducacion.com/control-predictivo/mpc-en-variables-de-estado/

Sergio A. Castaño Giraldo

¡Hola! Soy Sergio Andrés Castaño Giraldo. Bienvenidos a Control Automático Educación, un espacio dedicado a una de mis grandes pasiones: la Ingeniería de Control y la Automatización.

Mi objetivo aquí es tender un puente entre la teoría académica y la práctica profesional. Como profesor e investigador, he tenido la oportunidad de trabajar en proyectos complejos —desde sistemas de compresión en plataformas offshore hasta modelos de control avanzado— y todo ese aprendizaje lo vuelco en este sitio.

Aquí encontrarás material que crece constantemente: desde fundamentos de instrumentación hasta estrategias de control avanzado y programación profesional. Mi meta es que la información que comparto te sea de utilidad real en tu carrera o formación.

Te invito a que te suscribas, te unas a nuestro canal de YouTube y formemos una comunidad de aprendizaje continuo. ¡Gracias por estar aquí!

MPC por Variables de Estado

FUNCIÓN OBJETIVO

Caso con Restricciones

Código de Implementación (MATLAB code)

Citar este Artículo

Comparte esto: