Simulating a Nonlinear model predictive controller for a nonlinear system

Question

naiva saeedia 2025-3-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175555-simulating-a-nonlinear-model-predictive-controller-for-a-nonlinear-system

编辑： naiva saeedia 2025-3-24

Hi guys. I want to simulate the work of the paper in the attachments. I have simulated the open-loop model in CHO function as follows:

function state_dot = CHO(u,x)
%Extracellular parameters
miu_max = 0.0522;
miu_death_max = 0.03;
K_Glc = 0.75;
K_Gln = 0.075;
K_i_Lac = 85.88;
K_i_Amm = 14.24;
K_death_Amm = 1.76;
n = 2;
Y_Xv_Glc = 2.6e+8;
Y_Xv_Gln = 8e+8;
m_Glc = 9.8e-14;
alpha_1 = 3.4e-13;
alpha_2 = 4;
Y_Lac_Glc = 2;
Y_Amm_Gln = 0.45;
K_deg_Gln = 9.6e-3;
%Reduced model parameters
Y_mAb = 4.1331583e-9;
m_mAb = 1.212663e-8;
%States
V = x(1);
X_v = x(2);
Glc = x(3);
Gln = x(4);
Lac = x(5);
Amm = x(6);
mAb = x(7);
%Input
Q_in =u(1); %2.34375e-4;
Glc_in = 500;
Gln_in = 100;
%Auxilary equations
miu = miu_max*Glc/(K_Glc+Glc)*Gln/(K_Gln+Gln)*K_i_Lac/(K_i_Lac+Lac)*K_i_Amm/(K_i_Amm+Amm);
miu_death = miu_death_max/(1+(K_death_Amm/Amm)^n);
q_Glc = -miu/Y_Xv_Glc-m_Glc;
m_Gln = alpha_1*Gln/(alpha_2+Gln);
q_Gln = -miu/Y_Xv_Gln-m_Gln;
q_Lac = -Y_Lac_Glc*q_Glc;
q_Amm = -Y_Amm_Gln*q_Gln;
q_mAb = Y_mAb*miu + m_mAb;
V_dot = Q_in;
Xv_dot = (miu-miu_death)*X_v-Q_in*X_v/V;
Glc_dot = q_Glc*X_v-Q_in*(Glc-Glc_in)/V;
Gln_dot = q_Gln*X_v-Q_in*(Gln-Gln_in)/V-K_deg_Gln*Gln;
Lac_dot = q_Lac*X_v-Q_in*Lac/V;
Amm_dot = q_Amm*X_v-Q_in*Amm/V+K_deg_Gln*Gln;
mAb_dot = q_mAb*X_v-Q_in*mAb/V;
state_dot = [V_dot, Xv_dot, Glc_dot, Gln_dot, Lac_dot, Amm_dot, mAb_dot]';

and I disceretizing it using newton direct method with the help of this function:

function xk1 = Disc_CHO(xk, uk, Ts)
    substeps = 100;
    dt = Ts/substeps;
    % Initial state
    x_temp = xk;
    % Euler integration with sub-steps
    for i = 1:substeps
        state_dot = CHO(uk,x_temp);
        x_temp = x_temp + dt*state_dot;
    end
    % Output the state at the end of the sampling interval
    xk1 = x_temp;
end

Next I created a NLMPC object for the controller as follows with the custom cost function that has mentioned in the paper and its constraint for single feed case:

clear
clc
%Initial conditions
V0 =0.2;
Xv0=2e8;
Glc0=25.1;
Gln0=5.01;
Lac0=0;
Amm0=0;
mAb0=100;
Q0 =0;
x0=[V0;Xv0;Glc0;Gln0;Lac0;Amm0;mAb0];
u0 = Q0;
nx = 7;   %Number of states
ny = 1;   %Number of outputs (mAb)
nu = 1;   %Input: Qin
Ts = 12; %hours
PredictionHorizon = 14;
ControlHorizon = 7;
%Scaling states
% Assuming you have 7 states
nlobj.States(1).ScaleFactor = 1;        % V
nlobj.States(2).ScaleFactor = 1e8;      % Xv (very large)
nlobj.States(3).ScaleFactor = 100;      % Glc
nlobj.States(4).ScaleFactor = 10;       % Gln
nlobj.States(5).ScaleFactor = 10;       % Lac
nlobj.States(6).ScaleFactor = 1;        % Amm
nlobj.States(7).ScaleFactor = 200;      % mAb
nlobj.MV(1).ScaleFactor = 1e-4;  % if Q_in is small
nlobj.OutputVariables(1).ScaleFactor = 100;  % if using mAb as output
%Input constraints
nlobj.MV.Min = 0;
nlobj.MV.Max =2e-3;
nlobj.Weights.ManipulatedVariablesRate=1e-5;
%Create NLMPC object
nlobj = nlmpc(nx, ny, nu);
nlobj.Ts = Ts;
nlobj.PredictionHorizon = PredictionHorizon;
nlobj.ControlHorizon = ControlHorizon;
nlobj.Optimization.ReplaceStandardCost = true;
%State and output function
nlobj.Model.StateFcn = @(x,u) Disc_CHO(x, u, Ts);
nlobj.Model.OutputFcn = @(x,u) x(7);
nlobj.Model.IsContinuousTime = false;
%Custom cost Function (Maximize mAb final production)
nlobj.Optimization.CustomCostFcn = @(X,U,e,data) 1/(X(end,1)*X(end,7));
%nlobj.Optimization.CustomCostFcn = @(X, U, e, data) 1/((X(end,1)/0.2) * (X(end,7)/200));
%Custom constraint function clearly defined for volume limit
V_max = 0.24;
nlobj.Optimization.CustomIneqConFcn = @(X,U,e,data) X(end,1)-V_max;
validateFcns(nlobj, x0, u0)

My problem is the mv of the controller is not changing at all it looks like the optimal solution is Q_in = 0 but in the paper for single case you will see clearly this is not the case in this plot:

Also my simulation in simulink looks like this:

I don't have any refrence signal for this controller so I assumed it as constant. I only want to meet the conditions and maximize the objective function that the paper has mentioned. Also my simulation for open_loop model is correct based on the open loop plots for the single feed case in the paper.