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:

MPC advanced control applied to a gas compression system of an offshore platform

Entradas Analógicas

Como Hacer un Control PID de Temperatura con Arduino

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.