The control objective is to control the output temperature T by manipulating the input temperature Tc, the code uploaded below comprises of the NMPC code which is for the controller and Moving Horizon Estimation which is for the estimator that is estimating the state, the main issue here is that the code is working but not tracking setpoint. the response i got will be uploaded below to see that the response was actually tracking setpoint before it went off radar at t = 20 as since from the response. in correcting this i have try alot of tunning in both the controller and the estimator t remove the offset but all effort proves abortive.
NMPC_CSTR CODE
clc, clear all
%% Simulink Model
% To run this example, Simulink(R) is required.
if ~mpcchecktoolboxinstalled('simulink')
disp('Simulink is required to run this example.')
return
end
%%initialize NMPC object and set properties
nx = 4;
ny = 1;
%'UD', [5,6,7,8,9,10]
nlobj = nlmpc(nx, ny, "MV", 1, 'UD', [2,3]);
%The prediction model sample time is the same as the controller sample time
Ts = 0.05; %Ts = 2
nlobj.Ts = Ts;
%Prediction and control horizon
nlobj.PredictionHorizon = 5; % Increased for better long-term prediction
nlobj.ControlHorizon = 1;
%weights on output and manipulated variables
nlobj.Weights.OutputVariables = 65; % Balanced weight for faster tracking %0.1
nlobj.Weights.ManipulatedVariables = 0.1; % Increased to reduce control effort penalty %0.05;
nlobj.Weights.ManipulatedVariablesRate = 10; % Reduced to improve responsiveness
%constraints on manipulated variables
nlobj.MV(1).RateMin = -10;
nlobj.MV(1).RateMax = 10;
nlobj.MV(1).Min = 280;
nlobj.MV(1).Max = 350;
%constraint on output variables to keep them within reasonable bounds
nlobj.OV(1).Min = 300;
nlobj.OV(1).Max = 360;
%% specify the nonlinear state and output functions
nlobj.Model.StateFcn = 'cstrStateFcnCT';
nlobj.Model.OutputFcn = 'cstrOutputFcn';
%% validate the nonlinear MPC object with initial conditions
y0 = [0.877;324.5;0;0]; %[0.0033;920.86; 0.0104; 0.0085; 0.0118; 798];
u0 = 300;
validateFcns(nlobj, y0, u0);
MHE CODE
function [Xobs, d_opt, dist_output, fcost, Ecount, exitflag] = MHE_compute_testing(uo, um, ymeas, yplant, hstep, Nsteps, Nmhe, Rp, Rx, Ry, p_initial, p_MIN_colvector, p_MAX_colvector)
%disp('p_initial'); disp(p_initial);
%disp('p_MIN_colvector'); disp(p_MIN_colvector);
%disp('p_MAX_colvector'); disp(p_MAX_colvector);
% Persistent variable for xest_prev to retain value between calls
persistent PIo xest_prev
if isempty(xest_prev)
xest_prev = ymeas; % Use the measured value as the initial guess dynamically % Initial guess for CA and T
end
if isempty(PIo)
% Define system parameters
PIo = 1e-4 * eye(2); % Increased process noise covariance for better disturbance handling % Define PIo here as well
end
% Optimization options
OPTIONS = optimset('Algorithm', 'active-set', 'Display', 'iter');
tic
% Run fmincon with the step-by-step obj_est
[p_opt, fval, exitflag] = fmincon(@(p) obj_est(p, um, hstep, Nsteps, Nmhe, ymeas, xest_prev, PIo, Rp, Ry), ...
p_initial, [], [], [], [], p_MIN_colvector, p_MAX_colvector, [], OPTIONS);
% Update xest_prev for next function call
% xest_prev = p_opt(1:2); % Set xest_prev to the optimized state estimate
%disp('p_opt'); disp(p_opt)
fcost = fval;
Ecount = toc;
[Xobs] = state_cal(um,hstep, Nsteps, p_opt,Nmhe);
%disp('xobs'); disp(Xobs)
ypredicted = Xobs(2:2);
% Set outputs based on p_opt results
%Xobs = p_opt(1:2); % Assign first two elements of p_opt as Xobs
%d_opt = p_opt(3:4); % Assign next two elements as d_opt
%disp('yplant'); disp(yplant);
%disp('ypredicted'); disp(ypredicted)
dist_output = (yplant - ypredicted); % Dummy calculation for testing
dx1 = p_opt(Nmhe);
dx2 = p_opt(2*Nmhe);
d_opt = [dx1;dx2];
% Confirm p_opt result
%disp("Result from fmincon, p_opt:");
%disp(p_opt);
PIo = co_var(Xobs,uo,hstep,PIo,Rx,Ry);
end
%function stop = outputFcn(x, optimValues, state)
%stop = false; % Setting to false allows optimization to continue
%disp('Current value of p (optimization variable):');
%disp(x);
% Optionally, you can log values of x to a file or workspace
% You could use assignin('base', 'p_values', x) to track in base workspace
%end
function [Xobs] = state_cal(um,hstep,Nsteps,p_opt,Nmhe)
xpred0 = p_opt(2 * Nmhe + 1:2 * Nmhe + 2);
%disp('xpred');disp(xpred0)
%disp(um);disp(hstep);disp(Nsteps)
%disp('xpred0'); disp(xpred0)
%disp(size(p_opt))
xpred(:,1) = RK4(xpred0,um(:,1),hstep,Nsteps,p_opt(1),p_opt(Nmhe+1));
for kk = 2:Nmhe
xpred(:,kk) = RK4(xpred(:,kk-1),um(:,kk),hstep,Nsteps,p_opt(kk),p_opt(Nmhe+kk));
end
%disp('xpred'); disp(xpred)
[Xobs] = xpred(:,Nmhe);
%disp('xobs'); disp(Xobs)
end
function PIk = co_var(x, u, hstep, PIo, Rx, Ry)
% Define system parameters
F = 100;
V = 100;
k0 = 7.2e10;
E_R = 8750;
DeltaH = 5e4;
rho = 1000;
Cp = 0.239;
UA = 5e4;
CAF = 1.0;
T_f = 350;
% Linearized matrix A
CA = x(1);
T = x(2);
%input parameter
Tc = u;
Ac = [(-F/V - k0 * exp(-E_R / T) - CA * k0 * (E_R / T^2) * exp(-E_R / T)), (CA * k0 * (E_R / T^2) * exp(-E_R / T));
(-DeltaH * k0 * exp(-E_R / T)) / (rho * Cp), (-F/V - (DeltaH * CA * k0 * E_R / (rho * Cp * T^2)) * exp(-E_R / T) - UA / (V * rho * Cp))];
A = expm(Ac*hstep);
C = [0, 1]; % Measurement matrix, assuming direct state observation
PIk = Rx + A * PIo * A' - A * PIo * C' * (Ry + C * PIo * C')^(-1) * C * PIo * A';
end
function [xpredode] = RK4(xo, uo, hstep, Nsteps, dx1,dx2)
%disp('xo');disp(xo)
xk1 = xo;
%disp('xk1');disp()
for k = 1:Nsteps
k1 = hstep * odefun(xk1, uo, dx1,dx2);
k2 = hstep * odefun(xk1 + k1 / 2, uo, dx1,dx2);
k3 = hstep * odefun(xk1 + k2 / 2, uo, dx1,dx2);
k4 = hstep * odefun(xk1 + k3, uo, dx1,dx2);
xk1 = xk1 + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
end
xpredode = xk1;
disp('size of xpredode'); disp(size(xpredode))
%disp('xpredode');disp(xpredode)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xdot = odefun(x, u, dx1, dx2)
%values of parameters
%disp('x');disp(x)
F = 100; % Flow rate [L/min]
V = 100; % Volume [L]
k0 = 7.2e10; % Pre-exponential factor [min^-1]
E_R = 8750; % Activation energy divided by R [K]
DeltaH = 5e4; % Heat of reaction [J/mol]
rho = 1000; % Density [g/L]
Cp = 0.239; % Heat capacity [J/g/K]
UA = 5e4; % Heat transfer coefficient [J/min/K]
CAF = 1.0; % Feed concentration [mol/L]
T_f = 350; % Feed temperature [K]
% States
CA = x(1);
%T = x(2);
T = max(x(2), 1e-3);
%input
Tc = u;
% ODEs based on provided equations
dCAdt = F * (CAF - CA) / V - CA * k0 * exp(-E_R / T) + dx1;
dTdt = F * (T_f - T) / V + (-DeltaH) * CA * k0 * exp(-E_R / T) / (rho * Cp) + UA * (Tc - T) / (V * rho * Cp) + dx2;
xdot = [dCAdt; dTdt];
%disp('xdot');disp(xdot)
end
% Step-by-step obj_est for debugging
function er_fun = obj_est(p, um, hstep, Nsteps, Nmhe, ymeas, xest_prev, PIo, Rp, Ry)
%disp('p values in obj_est'); disp(p);
[ypred] = est_fun(um,p,hstep,Nsteps,Nmhe);
%disp(size(ypred(end,:)))
%disp('ymeas');disp(size(ymeas))
ypred = ypred(end,:)';
ymeas = ymeas';
%disp('ymeas'); disp(ymeas);
%disp('ypred'); disp(ypred);
%disp('ymeas-ypred'); disp(ymeas-ypred);
PI_p = Rp^-1;
PI_y = Ry^-1;
%disp(PI_y)
%PIo = diag([1e-6, 1e-6]);
%disp('PIo');disp(PIo)
%disp('Ry');disp(Ry)
%disp('xest_prev'); disp(xest_prev);
c_arrival = [(p(2*Nmhe + 1)-xest_prev(1)); (p(2*Nmhe + 2)-xest_prev(2))];
dest = p(1:2*Nmhe);
pred_cost = (ymeas-ypred)'*PI_y*(ymeas-ypred);
er_fun = pred_cost+c_arrival'*PIo^(-1)*c_arrival + 1e6*(dest)'*(dest);
end
function [ypred2] = est_fun(um,p,hstep,Nsteps,Nmhe)
%disp('p'); disp(p)
xm = p(2*Nmhe+1:2*Nmhe+2);
%disp('Input to RK4:'); disp(xm); disp(um(:,1));
xpred(:,1) = RK4(xm,um(:,1),hstep,Nsteps,p(1),p(Nmhe+1));
for kk=2:Nmhe
xpred(:,kk) = RK4(xpred(:,kk-1),um(:,kk),hstep,Nsteps,p(kk),p(Nmhe+kk));
end
%disp('Output from RK4:'); disp(xpred(:,1));
ypred2 = xpred;
end
CSTR_STATE_FUNCTION
function dydt = cstrStateFcnCT(y, u)
%definig the parameters
F = 100; % Flow rate [L/min]
V = 100; % Volume [L]
k0 = 7.2e10; % Pre-exponential factor [min^-1]
E_R = 8750; % Activation energy divided by R [K]
delta_H = 5e4; % Heat of reaction [J/mol]
rho = 1000; % Density [g/L]
Cp = 0.239; % Heat capacity [J/g/K]
U_A = 5e4; % Heat transfer coefficient [J/min/K]
Caf = 1.0; % Feed concentration [mol/L]
Tf = 350; % Feed temperature [K]
%%% state parameters
Ca = y(1);
T = y(2);
%%% input parameters
Tc = u(1);
dy1 = y(3);
dy2 = y(4);
dydt = zeros(4,1);
% Rate constant (Arrhenius equation)
k = k0 * exp(-E_R / T);
% Mass balance for Ca
dydt(1) = F/V * (Caf - Ca) - k * Ca + dy1;
% Energy balance for T
dydt(2) = F/V * (Tf - T) + (delta_H / (rho * Cp)) * k * Ca + (U_A / (rho * Cp * V)) * (Tc - T) + dy2;
dydt(3) = u(2);
dydt(4) = u(3);
CSTR_OUTPUT_FUNCTION
function k = cstrOutputFcn(y, u)
k = y(2);
