Control NMPC – Matlab – Simulink

En esta entrada vamos a diseñar e implementar un controlador predictivo no lineal basado en modelo NMPC para el laboratório de control de temperatura TCLAB utilizando el toolbox de MATLAB/Simulink.

Antes de comenzar, te hago la invitación para que te suscribas al canal de YouTube y que mires también nuestras entradas del Curso de Control Predictivo Basado en Modelo.

Proceso No Lineal

El proceso no lineal que deseamos controlar en esta entrada corresponde a un laboratório de control de temperatura (TCLAB) mostrado a continuación:

El TCLAB consta de dos calentadores (Transistores) los cuales se aumentan su temperatura en función a la corriente que fluye por ellos. Esta tarjeta junto con sus modelos lineales y no lineales y códigos de manipulación de la misma en matlab, simulink y python fueron propuestos por el sitio web APMonitor.

Las dos ecuaciones diferenciales que describen el proceso a través de la convección, radiación y conducción vienen dados por:

Q_{C12}=UA_s(T_2-T_1)

Q_{R12}=\epsilon \sigma A (T_{2}^4-T_1^4)

mc_p\dfrac{d T_1}{d t}=\alpha_1 Q_1+UA(T_{\infty}-T_1)+\epsilon \sigma A (T_{\infty}^4-T_1^4)+Q_{C12}+Q_{R12}

mc_p\dfrac{d T_2}{d t}=\alpha_2 Q_2+UA(T_{\infty}-T_2)+\epsilon \sigma A (T_{\infty}^4-T_2^4)+Q_{C12}+Q_{R12}

Los parámetros del modelo fueron tomados de la pagina de APMonitor, modificando el coeficiente de transferencia de calor y el factor del calentador para Mi caso específico:

Sintonía del Controlador MPC (MPC Tuning)

MPC Toolbox usando Matlab y Simulink

Escalado de Sistemas vía Condicionamiento Mínimo

NMPC Toolbox Matlab – Simulink

Vamos a utilizar la propia herramienta de Matlab para crear e implementar fisicamente el controlador NMPC con el TCLAB.

Toda la información sobre el uso de esta funcionalidad se encunetra en el sitio oficial de Mathworks.

Inicialmente comenzamos creando el objeto nlmpc (Nonlinear model predictive controller)

Un controlador predictivo de modelo no lineal calcula los movimientos de control óptimos a través de un horizonte de predicción utilizando un modelo de predicción no lineal, una función de costo no lineal y unas restricciones no lineales.

Una vez es creado el archivo nlmpc como mostrado en el video al comienzo del post, procedemos a simular inicialmente el controlador directamente en el modelo del proceso (caso nominal sin error de modelado) para ver el comportamiento del sistema con el controlador NMPC.

Posteriormente se hace la implementación, modificando la S-Function del diagrama anterior (S_TCLAB) por la planta real.

Para eso usamos los propios bloque sumunistrados por APMonitor para conseguir comunicación entre Simulink y el TCLAB. El esquema es el siguiente:

Finalmente, hacemos una comparación entre la respuesta del sistema controlando el modelo (caso nominal) y la respuesta del sistema controlando el TCLAB.

El caso nominal (rojo) como el caso en el sistema real (negro) son bastante próximos y muestra como usando el toolbox de matlab, conseguimos realizar una implementación del controlador NMPC en el sistema TCLAB.

Códigos

Los códigos y todos los archivos de simulink y simulación puedes descargarlo a continuación

[sociallocker id=»948″]

clc
clear all
close all

%% Nonlinear model
 % Parameters
Ta = 26 + 273.15;   % K
U = 2.3367;           % W/m^2-K
m = 4.0/1000.0;     % kg
Cp = 0.5 * 1000.0;  % J/kg-K    
A = 10.0 / 100.0^2; % Area in m^2
As = 2.0 / 100.0^2; % Area in m^2
alpha1 = 0.010364;    % W / % heater 1
alpha2 = 0.007126;    % W / % heater 2
eps = 0.9;          % Emissivity
sigma = 5.67e-8;    % Stefan-Boltzman

Par = [Ta;U;m;Cp;A;As;alpha1;alpha2;eps;sigma];
    
%% Crear objeto NMPC
nx = 2; %Estados
ny = 2;  %Salidas
nu = 2; %entradas
nlobj = nlmpc(nx,ny,nu);
nlobj.States(1).Name = 'T1';
nlobj.States(2).Name = 'T2';
nlobj.MV(1).Name = 'H1';
nlobj.MV(2).Name = 'H2';

nlobj.Model.StateFcn = @(x,u,Par) prediction_model_tclab(x,u,Par);
% No output function specified. Assuming "y = x" in the prediction model.
% nlobj.Model.OutputFcn  = 'arduino_ode';
nlobj.Model.IsContinuousTime = true;
nlobj.Model.NumberOfParameters=1;

%Jacobian
% syms T1 T2 Q1 Q2
% [aj,bj] = tclab_Jacobian(Par);
% 
% Ja = @(x,u,Par)sym_Jacobian(x,u,aj,bj);
% nlobj.Jacobian.StateFcn = @(x,u,Par)tclab_Jacobian(x,u,Par);


%% NMPC Tuning
N=[14  14];
Nu=[6  6];
delta=[0.69408     0.49601];
lambda=[0.24209    0.016057];
Ts = 8;
nlobj.Ts = Ts;                                  % Sample time
nlobj.PredictionHorizon = N(1);                    % Prediction horizon
nlobj.ControlHorizon =Nu(1);                       % Control horizon
nlobj.Weights.OutputVariables = delta;
nlobj.Weights.ManipulatedVariablesRate =lambda;

%% NMPC Constrains
nlobj.States(1).Min = 0;
nlobj.States(2).Min = 0;

nlobj.MV(1).Min = 0;    nlobj.MV(2).Min = 0;
nlobj.MV(1).Max = 100;  nlobj.MV(2).Max = 100;

nlobj.MV(1).RateMin = -100;   nlobj.MV(2).RateMin = -100;
nlobj.MV(1).RateMax = 100;    nlobj.MV(2).RateMax = 100;

%% Y ref
zeta1 = 0.5;
Wn1 = 4/(zeta1*180);
zeta2 = 0.5;
Wn2 = 4/(zeta2*180);
Pref=[tf(Wn1^2,[1 2*zeta1*Wn1 Wn1^2]),0;0,tf(Wn2^2,[1 2*zeta2*Wn2 Wn2^2])];
Prefz=c2d(Pref,Ts,'zoh'); 

%% Validate Custom FunctionsTi1 = 23 + 273.15;
Ti1 = 23 + 273.15;
Ti2 = 23 + 273.15;
x0 = [Ti1;Ti2];

time = [0,1000];
Q1 = 40;
Q2 = 0;
u0 = [Q1;Q2];

[time,T] = ode45(@(t,x)tclab_ode(t,x,u0,Par),time,x0);

validateFcns(nlobj,x0,u0,[],{Par});

% mdl = 'NMPC_TCLAB_1';
mdl = 'NMPC';
open_system(mdl)
createParameterBus(nlobj,[mdl '/Nonlinear MPC Controller'],'myBusObject',{Par});
% open_system(mdl)

% mv = nlmpcmove(nlobj,x0,u0,zeros(1,4),[Tc C2H4Avalability]);

%% Result
% 

load('NMPC_Tuning2');

Yref=lsim(Pref,out.Wt-(Ta - 273.15),out.tt,'zoh')+ones(1201,2).*(Ta-273.15);


%Tamanho da letra
ftsz1 = 24;
ftsz2 = 24;

figure
subplot(221)
plot(out.tt,Yref(:,1),'-r',out.tt,out.Yt(:,1),'-k',out.tt,out.Wt(:,1),'--','Linewidth',2);
xlabel('Time (s)','fontsize',ftsz2);
ylabel('°C','fontsize',ftsz2);
title('T_1','fontsize',ftsz2)
set(gca,'fontsize',ftsz1)

subplot(223)
stairs(out.tt,out.Qt(:,2),'k','Linewidth',2); 
xlabel('Time (s)','fontsize',ftsz2);
ylabel('%','fontsize',ftsz2);
title('H_1','fontsize',ftsz2)
set(gca,'fontsize',ftsz1)

subplot(222)
plot(out.tt,Yref(:,2),'-r',out.tt,out.Yt(:,2),'-k',out.tt,out.Wt(:,2),'--','Linewidth',2);
xlabel('Time (s)','fontsize',ftsz2);
ylabel('°C','fontsize',ftsz2);
title('T_2','fontsize',ftsz2)
set(gca,'fontsize',ftsz1)
 
subplot(224)
stairs(out.tt,out.Qt(:,2),'k','Linewidth',2); hold on
xlabel('Time (s)','fontsize',ftsz2);
ylabel('%','fontsize',ftsz2);
title('H_2','fontsize',ftsz2)
set(gca,'fontsize',ftsz1)

[/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.

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í!