Simulation of control system with only matlab script withou simulink

I want to simulate coontrol system using PID for a start. I have sets of DAEs. I have tried to simulate closed loop control system without using simulink blocks. But I got an error. I think my script has structural errors and I feel i am missing some things. I have alwys been using Simulink but I this time, I need to do some things and my research group are not familiar with simulink. So I have been requested to write the code uisng only script so that they can understand it since they know python.

Please, I need help on fixing this script given below. I have tried but I am stucked. I want to have freedom in selecting solver choice like ode 45,ode23s etc, I also want to be able to use any reference input for transient simulation.

function solveReactorDAE()
%% Parameters
Sig_x=2.65e-22;
yi=0.061;
yx=0.002;
lamda_x = 2.09e-5;
lamda_I = 2.87e-5;
Sum_f = 0.3358;
%nv=2.2e+3; % average velocisty of neutrons (thermal)
%nu=2.43; % the total number of neutrons liberated per rx
%Sigf=1/(gen*nv*nu); % what is this 
l=1.68e-3;
beta = 0.0065;
beta_1 = 2.267510e-4;
beta_2 = 1.22023e-3;
beta_3 = 1.193061e-3;
beta_4=2.785813e-3;
beta_5=1.257395e-3;
beta_6=5.226188e-4;
Lamda_1 = 0.0124;
Lamda_2 = 0.0369;
Lamda_3 = 0.632;
Lamda_4=3.044508e-1;
Lamda_5=8.539933e-1;
Lamda_6=2.868585;
%Gr8   =-2350E-5/180*2; % All drums
Gr4   =-1450E-5/180*2; % half
%Gr2   =-660E-5/180*2; % 1/4
Gr1   =-250E-5/180*2; % one
cp_f=977;
cp_m=1697;
cp_c=5188.6;
M_f=2002;
M_m=11573;
M_c=500;

mu_f=M_f*cp_f;
mu_m=M_m*cp_m;
mu_c=M_c*cp_c;
f_f = 0.96;
P_0 = 20;
Tf0=1105;
Tm0=1087;
T_in=864;
T_out=1106;
Tc0=(T_in+T_out)/2;
K_fm=f_f*P_0/(Tf0-Tm0);
K_mc=P_0/(Tm0-Tc0);
M_dot=1.75E+01;
alpha_f=-2.875e-5;
alpha_m=-3.696e-5;
alpha_c=0.0;
X0=2.35496411413791e10;

% Define your DAE system (example only)
%function [dx, algebraic] = reactorDAE(t, x, u, rho_inf)
    function [dx, algebraic] = reactorDAE(t, x, u1, u4)
        % Extract state variables
        %neutron_flux = x(1);
        %temperature = x(2);
        %xenon = x(3);
        n_r = x(1); Cr1 = x(2); Cr2 = x(3); Cr3 = x(4); Cr4 = x(5);Cr5=x(6);Cr6=x(7);
        X = x(8); I = x(9); Tf = x(10); Tm = x(11); Tc = x(12); Rho_d1 = x(13);Rho_d2 = x(14);
        
        % Define inputs (control signals) - here, for simplicity, assume constant
        %u1 = 0.01;          % Control signal 1 (e.g., PID output)
        %u2 = 0.005;         % Control signal 2 (e.g., temperature control)
        kp = 1;                 % Proportional gain
        ki = 0.1;               % Integral gain
        kd = 0.01;              % Derivative gain
        % u1=Kp * error + Ki * error_integrated + Kd * (error - error_prev) / dt;
        % %u2=Kp * error + Ki * error_integrated + Kd * (error - error_prev) / dt;
        % %u3=Kp * error + Ki * error_integrated + Kd * (error - error_prev) / dt;
        % u4=Kp * error + Ki * error_integrated + Kd * (error - error_prev) / dt;
        u1=pid(kp,ki,kd);
        %u2=pid(kp,ki,kd);
        %u3=pid(kp,ki,kd);
        u4=pid(kp,ki,kd);
        
        
        
        
        % Compute reactivity
        %rho = u(1) + u(2) - rho_inf * (neutron_flux - 1) + alpha_T * (temperature - 300) + alpha_Xe * xenon;
        
        
        rho = Rho_d1+Rho_d2+alpha_f*(Tf-Tf0)+alpha_c*(Tc-Tc0)+alpha_m*(Tm-Tm0)-Sig_x*(X-X0)/Sum_f;
        % Compute derivatives
        %dx = zeros(3, 1);
        %dx(1) = (beta / Lambda) * rho; % Neutron flux dynamics
        %dx(2) = 0.1 * (neutron_flux - 1); % Temperature dynamics (example only)
        %dx(3) = 0.0001 * (neutron_flux - 1); % Xenon dynamics (example only)
        
        %define the initial states
        dx = zeros(14,1) ;
        %%kinetics equations with six-delayed neutron groups
        dx(1) = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
        dx(2) = Lamda_1*n_r-Lamda_1*Cr1;
        dx(3) = Lamda_2*n_r-Lamda_2*Cr2;
        dx(4) = Lamda_3*n_r-Lamda_3*Cr3;
        dx(5) = Lamda_4*n_r-Lamda_4*Cr3;
        dx(6) = Lamda_5*n_r-Lamda_5*Cr3;
        dx(7) = Lamda_6*n_r-Lamda_6*Cr3;
        
        %% Xenon and Iodine dynamics
        G = 3.2e-11;
        V = 400*200;% i need to confirm the volume to be used 
        Pi = P_0/(G*Sum_f*V);
        % dx(8) = yx*Sum_f*Pi/nX0+lamda_I*I*nI0/nX0-Sig_x*X*Pi/nX0-lamda_x*X;
        % dx(9) = yi*Sum_f*Pi/nI0-lamda_I*I;
        dx(8) = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
        dx(9) = yi*Sum_f*Pi-lamda_I*I;
        %% thermal–hydraulics model of the reactor core
        dx(10)=f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
        dx(11)=(1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
        dx(12)=K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
        dx(13) = Gr1*u1;
        dx(14) = Gr4*u4;
        % Algebraic equations
        algebraic = [];
    end
% Define initial conditions
t_span = [0 100];
%x0 = [1; 300; 0];  % Initial neutron flux, temperature, xenon concentration
x0 = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032;
    0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 
    0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
% Define the setpoint
setpoint = 1.0;         % Desired neutron flux
neutron_flux_0=x0(1);
% Compute the error
%t_span=[t_start t_end];
%error = setpoint - neutron_flux_0;
%t_start = 0;
%t_end = 100;
%dt = 0.01;
%t = t_start:dt:t_end;
% Initialize arrays
%rho = zeros(size(t));
%neutron_flux = zeros(size(t));
%temperature = zeros(size(t));
%xenon = zeros(size(t));
%error_integrated = 0;
%error_prev = 0;

% Define function handle for DAE solver (e.g., ode15i for index-1 DAEs)
%dae_solver = @ode15i;
% Solve the DAE system
%[t, x] = dae_solver(@(t, x, u) reactorDAE(t, x, [u1; u4], rho_inf), t_span, x0, []);
[t, x] = ode15i(@(t, x, u1, u4) reactorDAE(t, x, u1, u4), t_span, x0, []);
% Plot results
figure;
subplot(4,1,1);
plot(t, x(1,:), 'r');
xlabel('Time (s)');
ylabel('Neutron Flux');
title('Neutron Flux vs Time');
% subplot(4,1,2);
% plot(t, x(2,:), 'g');
% xlabel('Time (s)');
% ylabel('Temperature (°C)');
% title('Temperature vs Time');
% 
% subplot(4,1,3);
% plot(t, x(3,:), 'b');
% xlabel('Time (s)');
% ylabel('Xenon concentration (ppm)');
% title('Xenon Concentration vs Time');
% 
% % Compute reactivity over time
% rho = zeros(size(t));
% for i = 1:length(t)
%     rho(i) = u1 + u2 - rho_inf * (x(1,i) - 1) + alpha_T * (x(2,i) - 300) + alpha_Xe * x(3,i);
% end
% 
% subplot(4,1,4);
% plot(t, rho, 'm');
% xlabel('Time (s)');
% ylabel('Reactivity (pcm)');
% title('Reactivity vs Time');

[Sam Chak]

Hi @Kamal,

I noticed that the script appears a bit untidy, making it difficult to understand your intentions regarding the solution of the primary DAE function named reactorDAE(t, x, u1, u4). I suggest editing the code to separate the two functions with proper indentation levels and adding 'end' at the end of each function. This will improve code readability and organization. Also try removing the 'commented-out' lines to further improve the readability. Later you can make wonderful plots once the code runs correctly.

The main issue with the code is likely an incompatibility problem. I noticed that you used the pid() command in u1 and u4 to generate the PID control signals. However, this approach is incorrect because the pid() command only produces a transfer function-like system as a function of a complex variable s, which operates in the complex-valued frequency domain. On the other hand, the DAE of your Reaction in the code is clearly described in the time domain, t.

Please request your control design team to convert the PID controller from the s-domain to an equivalent time-domain representation. This conversion is necessary to ensure consistency with the MATLAB code. Once the code is updated accordingly, please return here within the next few days.

kp  = 1;                % Proportional gain
ki  = 0.1;              % Integral gain
kd  = 0.01;             % Derivative gain
u1  = pid(kp, ki, kd)   % PID controller (it's a system, not a signal)

u1 =
 
             1          
  Kp + Ki * --- + Kd * s
             s          

  with Kp = 1, Ki = 0.1, Kd = 0.01
 
Continuous-time PID controller in parallel form.

[Kamal]

Hi @Sam Chak,

Thank you so much for the insight. I am the only one using MATLAB and doing this. The indentation issue is becasue I copied and pasted it. I have attached the file here instead of copying. Also, I had tried using time-domain PID if you observe that I commented it out. The issues are:

1. I don't know at what line to write the PID code, is it after the differentilal algebraic equation (DAEs) of my system or where?
2. I have two control inputs, that is why I wrote two PIDs.
3. States 13 and 14 are the inputs or used to input the control signals
4. The error is difference between the set point and output, initializing the first output to calculate the first error is also a challenge for me.

Pardon my ignorance, please because I am migrating from simulink to matlab. Once I get this, I believe I will be able to do subsequent ones

Thank for the inisght

[Sam Chak]

No worries, @Kamal. I suggest presenting the results to your team first and then inquire about their expectations, similar to the questions I raised in my Answer below.

Answer 1

[Sam Chak]

Hi @Kamal,

The Reactor system seems to exhibit stability even without any control inputs (u1 = u4 = 0). I'm curious about what specific aspect you intend to control. In other words, what are the desired output responses that you expect to observe in a practical scenario?

tspan   = [0 20];   % [0 500] to see the convergence of x2 and x3
x0      = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032; 0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
[t, x]  = ode45(@reactorDAE, tspan, x0);

figure(1)
tl1     = tiledlayout(3, 3, 'TileSpacing', 'Compact');
nexttile([1 3])
plot(t, x(:,1)), grid on, title('x_{1}')
nexttile
plot(t, x(:,2)), grid on, title('x_{2}')
nexttile
plot(t, x(:,3)), grid on, title('x_{3}')
nexttile
plot(t, x(:,4)), grid on, title('x_{4}')
nexttile
plot(t, x(:,5)), grid on, title('x_{5}')
nexttile
plot(t, x(:,6)), grid on, title('x_{6}')
nexttile
plot(t, x(:,7)), grid on, title('x_{7}')
title(tl1, 'Six neutron groups'), xlabel(tl1, 'Time (sec)')

figure(2)
tl2     = tiledlayout(2, 1, 'TileSpacing', 'Compact');
nexttile
plot(t, x(:,8)), grid on, title('x_{8}')
nexttile
plot(t, x(:,9)), grid on, title('x_{9}')
title(tl2, 'Xenon and Iodine'), xlabel(tl2, 'Time (sec)')

figure(3)
tl3     = tiledlayout(3, 3, 'TileSpacing', 'Compact');
nexttile
plot(t, x(:,10)),    grid on, title('x_{10}')
nexttile([2 2])
plot(t, x(:,10:14)), grid on, title('x_{10} to x_{14}')
nexttile
plot(t, x(:,11)),    grid on, title('x_{11}')
nexttile
plot(t, x(:,12)),    grid on, title('x_{12}')
nexttile
plot(t, x(:,13)),    grid on, title('x_{13}')
nexttile
plot(t, x(:,14)),    grid on, title('x_{14}')
title(tl3, 'Thermal–hydraulics states'), xlabel(tl3, 'Time (sec)')

function [dx, algebraic] = reactorDAE(t, x)
    %% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = -1450E-5/180*2; % half
    % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = -250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
    K_mc    = P_0/(Tm0-Tc0);
    M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;

    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
    
    G       = 3.2e-11;
    V       = 400*200;% i need to confirm the volume to be used 
    Pi      = P_0/(G*Sum_f*V);
    
    % Define inputs (control signals) - here, for simplicity, assume constant
    % u1 = 0.01;          % Control signal 1 (e.g., PID output)
    % u2 = 0.005;         % Control signal 2 (e.g., temperature control)
    % u1 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u2 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u3 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u4 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    
    kp      = 1;                 % Proportional gain
    ki      = 0.1;               % Integral gain
    kd      = 0.01;              % Derivative gain
    % u1    = pid(kp,ki,kd);
    % u2    = pid(kp,ki,kd);
    % u3    = pid(kp,ki,kd);
    % u4    = pid(kp,ki,kd);
    
    u1  = 0;    % stability test if the system is control-free
    u4  = 0;    % stability test if the system is control-free
        
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    
    %% ODEs
    dx    = zeros(14,1);
    
    %% kinetics equations with six-delayed neutron groups
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
    
    % Algebraic equations
    algebraic = [];
end

[Kamal]

%% Set reference
%dt=0.0001;
dt=0.1;
T = 2000;
n0=1;
time=0:dt:T;
nt=length(time);
ref=zeros(nt,1);
%time_point=[0	100	700	1000	1300	1600	1780	2080	2260	2560	2680	2980	3040	3340	3400	3700]
%pow_point=[1	1	0.9	0.9	0.8	0.8	0.65	0.65	0.95	0.95	0.55	0.55	0.15	0.15	0.95	0.95]
%time_point=[0 20 30 50 60 80 90 110 130 150];%work with only this time
%time_point=[0 10 40 70 80 100 110 130 140 150];
%time_point=[0 30 40 50 60 70 90 100 140 150];
time_point=[0 20 30 50 60 80 90 110 130 200];
%pow_point=[0.8 0.8 0.4 0.4 0.8 0.8 0.4 0.4 0.8 0.8];power Transient
%pow_point=[1 1 0.4 0.4 1 1 0.4 0.4 1 1];
%pow_point=[1 1 0.4 0.4 0.6 0.6 0.8 0.8 1 1];
%pow_point=[0.9 0.9 0.3 0.3 0.9 0.9 0.3 0.3 0.9 0.9];
%pow_point=[0.7 0.7 0.4 0.4 0.4 0.4 0.6 0.6 0.7 0.7];
pow_point=[0.5 0.5 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1];
ref_old=pow_point(1);
for it=1:nt %to vectorize nt or make each point of nt
    if(it>1)
        time(it)=time(it-1)+dt;%to be able to make equal time step of the time
        ref_old=ref(it-1); %for each time slot 
    end
    ref(it)=ref_old;
    i1=1;
    i2=1;
    for ii=1:length(time_point)-1
        if(time_point(ii)<=time(it) && time(it)<=time_point(ii+1))
            i1=ii;
            i2=ii+1;
            frac1=(time_point(ii+1)-time(it))/(time_point(ii+1)-time_point(ii));
            frac2=1.0-frac1;
            ref(it)=frac1*pow_point(i1)+frac2*pow_point(i2);
            break
        end
    end
end
ref(:)=ref(:);
yref=(ref.*n0);

@Sam Chak Thank you very much for the curiousity and the effort. I appreciate. Yes the reactor is marginally stable even without control. We have established that it has one pole in the negative side and one margnal pole at the origin. I want conduct transient scenarios as shown above. I have conducted the above transients with simulink already. In summary what I wnat to do is to incorporate the above reference inputs into the code and observe the transient response. I also want to input matched disturbance, conduct error performance of control algorithms I want to later develop for the same system and other tests. I had done all these with Simiulink already but my organization want the script because they dont use MATLAB. I am the only one using MATLAB.

Thank you and I appreciate.

[Sam Chak]

Hi @Kamal

I'm aware that you are the MATLAB coder in your team, while your other team members are the system designers. I plotted the reference input without interpreting the code. Is this what your Reactor team wants the control system to track?

Typically, I don't directly develop the code myself. Instead, I write out the governing mathematical equations for my budding coding scientists to grasp the underlying mathematics at a higher level. They can then learn coding independently while gaining a deeper understanding of the subject matter.

%% --- original code starts here ---
%% Set reference
dt=0.1;
T = 2000;
n0=1;
time=0:dt:T;
nt=length(time);
ref=zeros(nt,1);
time_point=[0 20 30 50 60 80 90 110 130 200];
pow_point=[0.5 0.5 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1];
ref_old=pow_point(1);
for it=1:nt %to vectorize nt or make each point of nt
    if(it>1)
        time(it)=time(it-1)+dt;%to be able to make equal time step of the time
        ref_old=ref(it-1); %for each time slot 
    end
    ref(it)=ref_old;
    i1=1;
    i2=1;
    for ii=1:length(time_point)-1
        if(time_point(ii)<=time(it) && time(it)<=time_point(ii+1))
            i1=ii;
            i2=ii+1;
            frac1=(time_point(ii+1)-time(it))/(time_point(ii+1)-time_point(ii));
            frac2=1.0-frac1;
            ref(it)=frac1*pow_point(i1)+frac2*pow_point(i2);
            break
        end
    end
end
ref(:)=ref(:);
yref=(ref.*n0);
%% --- original code ends here ---

plot(time, yref), ylim([-0.5, 1]), grid on
title('Reference input to be tracked by the control system')

[Kamal]

@Sam Chak Thank you. But the inputs above should give the figure belwo which is what I got using Simulink.

[Kamal]

@Sam Chak Yes the plot you got is one of the transients. Sorry I just looked again. Thank. I want to make sure other transienst above also works. It seems they will work once i change them. Thank you so much. Please, can I have the code so that I can run it and see some other transients. Also, If I want to calculate Integral absolute error, time absolute error, Integral time absolute, and absolute. And in case I want to pass disturbance into it, Please how do I go about this?

[Sam Chak]

Hi @Kamal

There are 14 states altogether as you can see. Please check with your Reactor team. Which state do they want to track the reference input?

Additionally, it is important to request your team to provide the error equation so that you can calculate the desired performance metrics such as ... "Integral absolute error (IAE), time absolute error (TAE), Integral time absolute (ITA?), and absolute (A?)" ...

Error,

[Kamal]

@Sam Chak I want to track the first state. The error equation is reference input -output. Error=nr-n, note that the first state is the nr so the difference is the error. Thank you

[Sam Chak]

Hi @Kamal,

Since I'm unfamiliar with your system, I took the liberty of "forcing" the 1st state (or ) to track the reference input, which I used as a square wave for testing purposes. However, in practice, your control design team should provide the appropriate control equation. They should already be well-versed in the system through their mathematical research before implementing any code for simulations.

tspan   = linspace(0, 20, 2001);
x0      = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032; 0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
[t, x]  = ode45(@reactorDAE, tspan, x0);

figure(1)
tl1     = tiledlayout(3, 3, 'TileSpacing', 'Compact');
nexttile([1 3])
plot(t, [x(:,1), 1/2*sign(sin(2*pi/10*t)) + 1/2]), grid on, ylim([-0.2, 1.2]), title('x_{1}'), legend('n_{r}', 'ref')
nexttile
plot(t, x(:,2)), grid on, title('x_{2}')
nexttile
plot(t, x(:,3)), grid on, title('x_{3}')
nexttile
plot(t, x(:,4)), grid on, title('x_{4}')
nexttile
plot(t, x(:,5)), grid on, title('x_{5}')
nexttile
plot(t, x(:,6)), grid on, title('x_{6}')
nexttile
plot(t, x(:,7)), grid on, title('x_{7}')
title(tl1, 'Six neutron groups'), xlabel(tl1, 'Time (sec)')

figure(2)
tl2     = tiledlayout(2, 1, 'TileSpacing', 'Compact');
nexttile
plot(t, x(:,8)), grid on, title('x_{8}')
nexttile
plot(t, x(:,9)), grid on, title('x_{9}')
title(tl2, 'Xenon and Iodine'), xlabel(tl2, 'Time (sec)')

figure(3)
tl3     = tiledlayout(3, 3, 'TileSpacing', 'Compact');
nexttile
plot(t, x(:,10)),    grid on, title('x_{10}')
nexttile([2 2])
plot(t, x(:,10:14)), grid on, title('x_{10} to x_{14}')
nexttile
plot(t, x(:,11)),    grid on, title('x_{11}')
nexttile
plot(t, x(:,12)),    grid on, title('x_{12}')
nexttile
plot(t, x(:,13)),    grid on, title('x_{13}')
nexttile
plot(t, x(:,14)),    grid on, title('x_{14}')
title(tl3, 'Thermal–hydraulics states'), xlabel(tl3, 'Time (sec)')

function [dx, algebraic] = reactorDAE(t, x)
    %% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = -1450E-5/180*2; % half
    % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = -250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
    K_mc    = P_0/(Tm0-Tc0);
    M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;

    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
    
    G       = 3.2e-11;
    V       = 400*200;% i need to confirm the volume to be used 
    Pi      = P_0/(G*Sum_f*V);
    
    % Define inputs (control signals) - here, for simplicity, assume constant
    % u1 = 0.01;          % Control signal 1 (e.g., PID output)
    % u2 = 0.005;         % Control signal 2 (e.g., temperature control)
    % u1 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u2 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u3 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u4 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    
    kp      = 1;                 % Proportional gain
    ki      = 0.1;               % Integral gain
    kd      = 0.01;              % Derivative gain
    % u1    = pid(kp,ki,kd);
    % u2    = pid(kp,ki,kd);
    % u3    = pid(kp,ki,kd);
    % u4    = pid(kp,ki,kd);
    
    u1  = 0;    % stability test if the system is control-free
    u4  = 0;    % stability test if the system is control-free
        
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    
    %% ODEs
    dx    = zeros(14,1);
    
    %% kinetics equations with six-delayed neutron groups

    n       = 1/2*sign(sin(2*pi/10*t)) + 1/2;   % square wave
    fx      = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    k       = 10;   % Proportional gain
    u       = - fx - k*(n_r - n);
    
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6 + u;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
    
    % Algebraic equations
    algebraic = [];
end

[Kamal]

@Sam Chak Thank you so much for going this far for me. I really appreciate. The code runs. I will try to apply the transients I supplied above. We actually understand the system and its mathematics. Please, let me know what you don dont understand in the system. I have worked with lots of transients with this system. I have designed PID, Super-twisting control and nonlinear MPC with the system just couple of weeks ago. Howevver, I used Simulink blocks except the function blocks that I wrote codes.

They want me to implement it with only scripts without blocks. You have saved me several man-hours. In fact let me attach all the transients I have plotted using Simulink. I need to replicate them again with scripts. So I will implement Sliding mode, PID and Nonlinear MPC now.

please

1. I noticed that the feedback control you used is 'u= - fx - k*(n_r - n)'; seems like only P and you feedback the state 1, which does not seem like PID. Is it that the PID is not working or cannot be impelemented as it is in Simulink?
2. what if I want to simulate for other control systems?
3. The fact that you were able to force it to follow the sqaurewave transinet indicates, I should be able to implement other transients. thank you

[Sam Chak]

Hi @Kamal

I'm a bit puzzled as to why your team requested you to convert the Simulink simulation into MATLAB code. The Simulink model was constructed based on the mathematical equations developed by the team, and the results appear reliable. It begs the question: why aren't they considering publishing their findings? As a coder, your role primarily involves validating that the control algorithms are correctly implemented based on the control equations derived by the design team on paper. Moreover, since they are proficient in Python, wouldn't it be more efficient for you as a coder to translate the Python code into MATLAB code?

Reply to Q1:

As clarified previously, it is a mathematical trick to enforce to track the reference input using a simple proportional gain structure. This is achieved by introducing an input signal u, where . The expression for u is given by , resulting in a closed-loop equation of . Although it is possible to introduce an additional input signal u directly into the dynamics, your control design team should determine the appropriate location for injecting u among the 14 differential equations. Alternatively, your Simulink model, in theory, can also provide insights into the suitable point for injecting the PID control signal.

Reply to Q2:

Simulating other control algorithms is not a problem once you have determined the appropriate location to inject the control input signal u. In the original code, you injected two control input signals, and , at the and dynamics, respectively. Since your team members have a good understanding of the system and have derived mathematical equations, they should be able to confirm whether these two injection points can enable to track the reference input or not.

Reply to Q3:

It depends. My approach incorporates a clever trick, perhaps an unconventional method to decouple the dynamics from the other 13 state dynamics. This results in , enabling the system to achieve any desired transient performance in theory. Based on the Simulink results displayed in the image (figure1.png), it appears that your designed controllers (PID, STSMC, NMPC) can rapidly achieve various transient performances.

[Kamal]

@Sam Chak You have really done justice to this system for me with explanation above. Than you so muc for the time and the effrot to clarify and convince me.

1. I designed the control algorithms and the system equation is from a system that is still under design verification and validation by a firm. Therefore, They want to undersatnd everything I am doing. I feel it is transparency. Since they dont know Simulink and they dont know the mathematics of Simulink. So they requested I write everything in Matlab script so that they can undertstand everything becasue they know Python. That is why. I believe it is good practice and in fact learning curve for me as well because I am very comforatble with Simulink lol.
2. States x13 and x14 are actually the control devices that recieve the control signal. I will say actuator in this case. Although we have lumped things together there such as motor, rotor and some other mechanisms
3. Please I want to verify some thing and if you have any paper to cite or any evidence you can please direct me. Is it reasonable to say that the higher the controller gain of aany control method, the larger the mechanisms required to implement such coontrol algoritm. For instance, a PID gain of Kp=1000, Ki=300, Kd=10 will require bigger motor and acuator than that Kp=100, Ki=30, Kd=1?
4. Thank you so much

[Sam Chak]

Hi @Kamal

First and foremost, it is crucial to verify whether the simulations of the control-free system yield the same output response in both the MATLAB code and Simulink model. Can you show that?

[Kamal]

@Sam Chak So sorry for my late response. I was away for couple of days. I will work on implementing the control algorithms I already did on Simulink through the weekend. I will update you here. I will also conduct the control-free system with simulink and show it here as requested. Thank you. The tasks will take my weekend.

[Kamal]

Hello @Sam Chak I have tried to apply control algorithms to the script as promised. However, It has not worked. I actually want the matlab file to contian

1. The paramters
2. The state declaration
3. reference input
4. The system DAEs
5. Then I want the control algoriths like PID, STC , MPC to be implemented
6. Finally I want the plots of closed loop system of each control algorithms plotted on the same graph

I have tried to work it out but I was getting error that 'a function defination must appear at the end of a file.......'' I have tried to create local functions for the DAEs, and control algorithms. I think I should create function for plotting the results as well.

The issue I have been having is implementing the PID for instance and passing the PID signal to the system to form a closed loop system that tracks the reference input I set? Below is the code. Could you please advise on what I can do to the code. I believe if I get the PID working, I can use the initiative to implement the STC and the MPC in the same script.

Thank you

%% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = -1450E-5/180*2; % half
     % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = -250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
     K_mc    = P_0/(Tm0-Tc0);
     M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;
     G       = 3.2e-11;
     V       = 400*200;
     Pi      = P_0/(G*Sum_f*V);
    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
  
    
    %n       = 1/2*sign(sin(2*pi/10*t)) + 1/2;   % square wave
    %% step response
    %yref = 0.8;
    %% transient on and off 
    % dt=0.0001; % Simulation time step
    % %T=24*3600*5;% Simulation time or period
    % T=150;
    % time=0:dt:T; % Simulation time slots
    % nt=length(time);
    % n0=1;
    % %hr=3600; % hour in seconds
    % zp=1e-30; % Almost zero to indicate zero power level 
    % time_point=0;
    % pow_point =1;
    % for i=1:30 % I think this is to obtain the reference power points with time
    %     %Overnight
    %     it=length(time_point);
    %     %time_point(it+1)=time_point(it)+300; %note the negligible rise in power at this time point
    %     time_point(it+1)=time_point(it)+0.003; 
    %     pow_point (it+1)=zp;% zero power point
    %     %time_point(it+2)=time_point(it)+8*hr;
    %     time_point(it+2)=time_point(it)+1;
    %     %pow_point (it+2)=zp;
    %     %Daytime
    %     it=length(time_point);
    %     %time_point(it+1)=time_point(it)+300;
    %     time_point(it+1)=time_point(it);
    %     %time_point(it+1)=time_point(it);
    %     pow_point (it+1)=1;
    %     %time_point(it+2)=time_point(it)+16*hr;
    %     time_point(it+2)=time_point(it)+2;
    %     pow_point (it+2)=1;
    % end
    % 
    % ref_old=pow_point(1); 
    % for it=1:nt % nt is length of the simulation time slots
    %     if(it>1)
    %         time(it)=time(it-1)+dt; % Simulation time step and putting it in  time slots
    %         ref_old=ref(it-1);
    %     end
    %     ref(it)=ref_old;
    %     i1=1;
    %     i2=1;
    %     for ii=1:length(time_point)-1
    %         if(time_point(ii)<=time(it) && time(it)<=time_point(ii+1)) %why this condition?
    %             i1=ii;
    %             i2=ii+1;
    %             frac1=(time_point(ii+1)-time(it))/(time_point(ii+1)-time_point(ii));
    %             frac2=1.0-frac1;
    %             ref(it)=frac1*pow_point(i1)+frac2*pow_point(i2);
    %             break
    %         end
    %     end
    % end
    % yref=(ref.*n0); 
    %% load following
    dt=1;
    T = 200;
    n0=1;
    time=0:dt:T;
    nt=length(time);
    ref=zeros(nt,1);
    %time_point=[0	100	700	1000	1300	1600	1780	2080	2260	2560	2680	2980	3040	3340	3400	3700]
    %pow_point=[1	1	0.9	0.9	0.8	0.8	0.65	0.65	0.95	0.95	0.55	0.55	0.15	0.15	0.95	0.95]
    %time_point=[0 20 30 50 60 80 90 110 130 150];%work with only this time
    %time_point=[0 10 40 70 80 100 110 130 140 150];
    %time_point=[0 30 40 50 60 70 90 100 140 150];
    time_point=[0 20 30 50 60 80 90 110 130 200];
    %pow_point=[0.8 0.8 0.4 0.4 0.8 0.8 0.4 0.4 0.8 0.8];power Transient
    %pow_point=[1 1 0.4 0.4 1 1 0.4 0.4 1 1];
    %pow_point=[1 1 0.4 0.4 0.6 0.6 0.8 0.8 1 1];
    %pow_point=[0.9 0.9 0.3 0.3 0.9 0.9 0.3 0.3 0.9 0.9];
    pow_point=[0.7 0.7 0.4 0.4 0.4 0.4 0.6 0.6 0.7 0.7];
    %pow_point=[0.5 0.5 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1];
    ref_old=pow_point(1);
    for it=1:nt %to vectorize nt or make each point of nt
        if(it>1)
            time(it)=time(it-1)+dt;%to be able to make equal time step of the time
            ref_old=ref(it-1); %for each time slot 
        end
        ref(it)=ref_old;
        i1=1;
        i2=1;
        for ii=1:length(time_point)-1
            if(time_point(ii)<=time(it) && time(it)<=time_point(ii+1))
                i1=ii;
                i2=ii+1;
                frac1=(time_point(ii+1)-time(it))/(time_point(ii+1)-time_point(ii));
                frac2=1.0-frac1;
                ref(it)=frac1*pow_point(i1)+frac2*pow_point(i2);
                break
            end
        end
    end
    ref(:)=ref(:);
    yref=(ref.*n0);
    %% kinetics equations with six-delayed neutron groups
    %ODEs
    dx    = zeros(14,1);
    function dx = reactorDAE(t, x)
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6 + u;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
    
    end
    % fx      = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    % k       = 10;   % Proportional gain
    % u       = - fx - k*(n_r - yref);
    
    %% PID controller 
    function u_pid= PID_reactor()
    
    % Define PID gains
    kp      = 1;                 % Proportional gain
    ki      = 0.1;               % Integral gain
    kd      = 0.01;              % Derivative gain
    
    
    % Initialize simulation parameters
    % dt = 0.01;  % Time step (s)
    % endTime = 10;  % Simulation end time (s)
    % time = 0:dt:endTime;  % Time vector
    % 
    % % Initialize variables for simulation
    % output = zeros(size(time));  % Output of the PID controller
    % input = sin(time);  % Example input signal (replace with actual sensor reading)
    previousError = 0;  % Previous error for derivative term
    integral = 0;  % Integral term
    % Simulation loop
    for i = 1:length(time)
        % Calculate the error signal
        setpoint = yref;  % Desired setpoint (replace with actual setpoint if variable)
        error = setpoint - n_r(i);
        
        % Update integral term
        integral = integral + error * dt;
        
        % Calculate derivative term
        if i == 1
            derivative = 0;
        else
            derivative = (error - previousError) / dt;
        end
        
        % Compute the PID controller output
        u_pid(i) = kp * error + ki * integral + kd * derivative;
        
        % Send the output to the actuator (replace with actual actuator control code)
        % actuatorControl(u_pid(i));
        
        % Update previous error
        previousError = error;
        
        % Wait for next time step (optional, for real-time simulation)
        pause(dt);
    end 
    end
    
    x0      = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032; 0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
    [t, x]  = ode45(@reactorDAE, tspan, x0);
    figure(1)
    tl1     = tiledlayout(3, 3, 'TileSpacing', 'Compact');
    nexttile([1 3])
    %plot(t, [x(:,1), 1/2*sign(sin(2*pi/10*t)) + 1/2]), grid on, ylim([-0.2, 1.2]), title('x_{1}'), legend('n_{r}', 'ref')
    plot(t, [x(:,1), yref]), grid on, ylim([-0.2, 1.2]), title('x_{1}'), legend('n_{r}', 'ref')
    nexttile
    plot(t, x(:,2)), grid on, title('x_{2}')
    nexttile
    plot(t, x(:,3)), grid on, title('x_{3}')
    nexttile
    plot(t, x(:,4)), grid on, title('x_{4}')
    nexttile
    plot(t, x(:,5)), grid on, title('x_{5}')
    nexttile
    plot(t, x(:,6)), grid on, title('x_{6}')
    nexttile
    plot(t, x(:,7)), grid on, title('x_{7}')
    title(tl1, 'Six neutron groups'), xlabel(tl1, 'Time (sec)')
    figure(2)
    tl2     = tiledlayout(2, 1, 'TileSpacing', 'Compact');
    nexttile
    plot(t, x(:,8)), grid on, title('x_{8}')
    nexttile
    plot(t, x(:,9)), grid on, title('x_{9}')
    title(tl2, 'Xenon and Iodine'), xlabel(tl2, 'Time (sec)')
    figure(3)
    tl3     = tiledlayout(3, 3, 'TileSpacing', 'Compact');
    nexttile
    plot(t, x(:,10)),    grid on, title('x_{10}')
    nexttile([2 2])
    plot(t, x(:,10:14)), grid on, title('x_{10} to x_{14}')
    nexttile
    plot(t, x(:,11)),    grid on, title('x_{11}')
    nexttile
    plot(t, x(:,12)),    grid on, title('x_{12}')
    nexttile
    plot(t, x(:,13)),    grid on, title('x_{13}')
    nexttile
    plot(t, x(:,14)),    grid on, title('x_{14}')
    title(tl3, 'Thermal–hydraulics states'), xlabel(tl3, 'Time (sec)')

[Sam Chak]

Hi @Kamal

It seems that you have made changes to the guidance plan by implementing the PID control scheme and then rearranged the MATLAB coding structure, resulting in the error message.

However, it was previously requested that you verify whether the system in MATLAB code produces the same result as the system in Simulink when no control input is applied. This step is crucial because if the two systems are not exactly the same, any control schemes implemented in the MATLAB code would yield invalid results.

To address this, you should take the following steps as shown in this example: disconnect the input signal of the controller and directly feed in the reference input, preferably using a Unit Step signal.

Result from Simulink model

Result from system in MATLAB

%% ----- Beginning of Script -----
tspan   = [0 10];
x0      = [0; 0];
[t, x]  = ode45(@mySystem, tspan, x0);
plot(t, x(:, 1)), grid on, xlabel('t'), ylabel('x(t)')
title('System''s Output Response')

%% ----- End of Script -----

%% System's function file placed at the end of the script
function dxdt = mySystem(t, x)
    ref       = 1;
    dxdt(1,1) = x(2);
    dxdt(2,1) = 1*( ref - x(1)) - 2*x(2);
end

[Kamal]

@Sam Chak yeah thank you. I modified it by testing If I could simulate closed loop system. Just to try my knowledge becasue I read some stuffs on MATLAB central and I was trying them out. Meanwhile I have implemented the system using DAEs without control and here is what I got. I commneted out the line where you implemented the feedback control and apply the step input as my u into the system. I used both ODE45 and ODE23s. Both gave simuiar results attached. I output the state 1 x(1), which is the output power of my system. Figure attacted below. Thank you.

[Sam Chak]

Hi @Kamal

Click on the indentation icon to copy/paste the MATLAB code into the gray field. You can run it online in MATLAB to verify if it produces the same output response.

%% --- Paste the code here ---

[Kamal]

%% --- added by Sam ---
tspan   = [0 10];
x0      = zeros(14, 1);
[t, x]  = ode45(@reactorDAE, tspan, x0);
plot(t, x(:,1)), grid on, 
xlabel('t'), ylabel('x_{1}'), title('Response of state x_{1}')

%% --- added by Sam ---

function [dx, n_r]= reactorDAE(t, x)    % <-- edited by Sam

%% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = -1450E-5/180*2; % half
    % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = -250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
    K_mc    = P_0/(Tm0-Tc0);
    M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;

    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
    
    G       = 3.2e-11;
    V       = 400*200;% i need to confirm the volume to be used 
    Pi      = P_0/(G*Sum_f*V);
    
    % Define inputs (control signals) - here, for simplicity, assume constant
    % u1 = 0.01;          % Control signal 1 (e.g., PID output)
    % u2 = 0.005;         % Control signal 2 (e.g., temperature control)
    % u1 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u2 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u3 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u4 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    
    kp      = 1;                 % Proportional gain
    ki      = 0.1;               % Integral gain
    kd      = 0.01;              % Derivative gain
    % u1    = pid(kp,ki,kd);
    % u2    = pid(kp,ki,kd);
    % u3    = pid(kp,ki,kd);
    % u4    = pid(kp,ki,kd);
    
    u1  = 0;    % stability test if the system is control-free
    u4  = 0;    % stability test if the system is control-free
        
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    
    %% ODEs
    dx    = zeros(14,1);
    
    %% kinetics equations with six-delayed neutron groups

    % n       = 1/2*sign(sin(2*pi/10*t)) + 1/2;   % square wave
    % fx      = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    % k       = 10;   % Proportional gain
    % u       = - fx - k*(n_r - n);
    u       = heaviside(t);
    
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6 + u;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
end

[Kamal]

@Sam Chak I have attached the simulink file and apsted the matlab code as instructed. please note that I applied step input as u. I dont know how to implement the step input without simulinlk

[Sam Chak]

Hi @Kamal,

I made some edits to your post to enable running the simulation directly from the provided code. This helps to keep the thread concise and avoid excessive code length. I added a few lines to ensure the code runs smoothly, without altering the system dynamics. The results obtained from both approaches appear identical. However, I noticed that you directly fed the signal u into the dynamics. I would suggest consulting your system modeling team to confirm whether this is accurate or not.

Else, any control schemes implemented later in the MATLAB code would yield invalid results!

[Kamal]

@Sam Chak Thank you I have run it and I confirm it gives the same result. Yeah I fed the u directly to the x1 dynamics as a test. However, in reality, the control input u is fed to x(13) and x(14) dynamics. The two states accept the control signal I am trying to design for the system. For instance, the PID signal or any other control signal will be fed to the two states to control the system to achieve the desired transient response.

[Sam Chak]

Hi @Kamal

I observed that the reactorDAE block in Figure 1 implements the exact dynamics as coded in MATLAB. This likely explains why the results matched in both MATLAB's ode45 and the Simulink's reactorDAE model. Could you please demonstrate how you and your research team generated the results shown in Figure 2 using two PID controllers?

If you are confident that the system model is accurate, you can try using two PID controller blocks to connect the reactorDAE system using the and channels in the Simulink model. By using the Desired Power profile as the Reference Input signal, you should obtain the same result as successfully simulated in Figure 2.

Figure 1: Reactor's dynamic system model in Simulink

Figure 2: Result of x1 under PID control simulated in Simulink.

[Kamal]

@Sam Chak Yes this exactly what I have been trying to explain. The simulink model has no issue for me. I do not want to use the simulink at all. I only want to use only MATLAB script with no simulink block to generate figure 2 above. My questions are:

1. can I wriite Matlab code for the PID algorthim to generate signals that I will feed into states x13 and x14 within the same matlab code above ?. It could be a single signal u_pid that I will use to multiply state x13 and x14.
2. Concerning how I generated figure 2, I used the PID block and tuned it to generate the signal that is fed into the system. The signal is used to multiply the state x13 and x14.

Now how do this affect the output, look into this equation

rho = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;

Note the first two terms on the right handside, you will notice that they are state x13 and x14. The state x13 and x14 are generated by multiplying the input signal with a constant as given here dx(13) = Gr1*u1; dx(14) = Gr4*u4 and integral gives x13 and x14, if you check the simulink implementation. The remaining terms are inherent feedbacks into the system. Now note the term on the left handside, you will see that it is the one in the state x(1); dx(1) = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6; That is how we influence the output state x1 with control input.

So it is the control input i want to generate with matlab script only WITHOUT SIMULINK BLOCKS.

In summary, what I want is to add control algrithm such as PID to the above MATLAB script that is working.

Note all the bolded

Answer 2

[Sam Chak]

Hi @Kamal

I've developed a custom function called 'pidController()' to mimic the PID controller based on the control equation provided by your Control Design team. This function takes {setpoint, Kp, Ki, Kd, dt} as additional parameters in the input argument. Scroll to the bottom of the script.

However, it appears that the chosen gain values in the PID controller have destabilized the Reactor system. It might be wise to double-check with your Control Design team to validate the PID control design and compare the results against the original Simulink model that you previously executed successfully for PID, SMC, and MPC, as illustrated in Figure 2.

%% Parameters
setpoint = 1;       % Reference input (actual one to be supplied by the user)
Kp       = 1;       % Proportional gain
Ki       = 0.1;     % Integral gain
Kd       = 0.01;    % Derivative gain
dt       = 1e-4;    % ideally it should be the time elapsed between the current and previous time steps

%% Call ode45
tspan    = [0 20];
x0       = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032; 0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
[t, x]   = ode45(@(t, x) reactorDAE(t, x, pidController(t, x, setpoint, Kp, Ki, Kd, dt)), tspan, x0);

Warning: Failure at t=3.009812e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t.

%% Plot result
figure(1)
plot(t, x(:,1)), grid on, xlabel('Time'), title('x_{1}')

%% Reactor dynamics
function [dx, algebraic] = reactorDAE(t, x, u)
    %% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = -1450E-5/180*2; % half
    % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = -250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
    K_mc    = P_0/(Tm0-Tc0);
    M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;

    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
    
    G       = 3.2e-11;
    V       = 400*200;% i need to confirm the volume to be used 
    Pi      = P_0/(G*Sum_f*V);
    
    %% Define inputs (control signals) - here, for simplicity, assume constant
    % u1 = 0.01;          % Control signal 1 (e.g., PID output)
    % u2 = 0.005;         % Control signal 2 (e.g., temperature control)
    % u1 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u2 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u3 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u4 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u1 = pid(kp,ki,kd);
    % u2 = pid(kp,ki,kd);
    % u3 = pid(kp,ki,kd);
    % u4 = pid(kp,ki,kd);
    
    %% The extra parameter u in reactorDAE(t, x, u) comes from the pidController()
    u1  = u;    % actuation thru channel x13, both PID controllers are the same 
    u4  = u;    % actuation thru channel x14
    
    %% ODEs
    dx      = zeros(14,1);
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    
    %% kinetics equations with six-delayed neutron groups
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
    
    %% Algebraic equations
    algebraic = [];
end

%% PID controller
function u = pidController(t, x, setpoint, Kp, Ki, Kd, dt)
    persistent integral_error prev_error prev_time;
    if isempty(integral_error)
        integral_error  = 0;
        prev_error      = 0;
        prev_time       = 0;
    end
    
    error               = setpoint - x(1);
    integral_error      = integral_error + error*(t - prev_time);
    derivative_error    = (error - prev_error)/dt;
    prev_error          = error;
    prev_time           = t;
    
    % user-supplied PID control equation
    u = Kp*error + Ki*integral_error + Kd*derivative_error;
end

[Kamal]

@Sam Chak; Thank you so much for this great effort. This is now taking the form of what I am looking for. I have attached the plot I got after changing the G parameter in state x13 and x14 from negative to positive. It has to do with poistive and negative feedback. But I need to remove the overshhoot and also implement the transient i have in the simulink. I will try it now.

%% Parameters
setpoint = 1;       % Reference input (actual one to be supplied by the user)
Kp       = 3000;       % Proportional gain
Ki       = 1000;     % Integral gain
Kd       = 50;    % Derivative gain
dt       = 1e-4;    % ideally it should be the time elapsed between the current and previous time steps

%% Call ode45
tspan    = [0 20];
x0       = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032; 0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
[t, x]   = ode45(@(t, x) reactorDAE(t, x, pidController(t, x, setpoint, Kp, Ki, Kd, dt)), tspan, x0);

%% Plot result
figure(1)
plot(t, x(:,1)), grid on, xlabel('Time'), title('x_{1}')

%% Reactor dynamics
function [dx, algebraic] = reactorDAE(t, x, u)
    %% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = 1450E-5/180*2; % half
    % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = 250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
    K_mc    = P_0/(Tm0-Tc0);
    M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;

    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
    
    G       = 3.2e-11;
    V       = 400*200;% i need to confirm the volume to be used 
    Pi      = P_0/(G*Sum_f*V);
    
    %% Define inputs (control signals) - here, for simplicity, assume constant
    % u1 = 0.01;          % Control signal 1 (e.g., PID output)
    % u2 = 0.005;         % Control signal 2 (e.g., temperature control)
    % u1 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u2 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u3 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u4 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u1 = pid(kp,ki,kd);
    % u2 = pid(kp,ki,kd);
    % u3 = pid(kp,ki,kd);
    % u4 = pid(kp,ki,kd);
    
    %% The extra parameter u in reactorDAE(t, x, u) comes from the pidController()
    u1  = u;    % actuation thru channel x13, both PID controllers are the same 
    u4  = u;    % actuation thru channel x14
    
    %% ODEs
    dx      = zeros(14,1);
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    
    %% kinetics equations with six-delayed neutron groups
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
    
    %% Algebraic equations
    algebraic = [];
end

%% PID controller
function u = pidController(t, x, setpoint, Kp, Ki, Kd, dt)
    persistent integral_error prev_error prev_time;
    if isempty(integral_error)
        integral_error  = 0;
        prev_error      = 0;
        prev_time       = 0;
    end
    
    error               = setpoint - x(1);
    integral_error      = integral_error + error*(t - prev_time);
    derivative_error    = (error - prev_error)/dt;
    prev_error          = error;
    prev_time           = t;
    
    % user-supplied PID control equation
    u = Kp*error + Ki*integral_error + Kd*derivative_error;
end

[Kamal]

Also, please I am trying to run the transients I set as my new setpoints as shown at the begining of the code below to get the same outputs like the one in my simulink but it seems I am missing soemthing in the 'for loop' I created. It is only capturing the first point of the setpoint. Please see the below and note the for loop I created. Observe that it captures only the initial point of 0.7 instead of the entire points pow_point=[0.7 0.7 0.4 0.4 0.4 0.4 0.6 0.6 0.7 0.7];

%% Parameters
 %% step response
     %yref = 0.8;
     %% transient on and off 

     % dt=0.0001; % Simulation time step
     % %T=24*3600*5;% Simulation time or period
     % T=150;
     % time=0:dt:T; % Simulation time slots
     % nt=length(time);
     % n0=1;
     % %hr=3600; % hour in seconds
     % zp=1e-30; % Almost zero to indicate zero power level 
     % time_point=0;
     % pow_point =1;
     % for i=1:30 % I think this is to obtain the reference power points with time
     %     %Overnight
     %     it=length(time_point);
     %     %time_point(it+1)=time_point(it)+300; %note the negligible rise in power at this time point
     %     time_point(it+1)=time_point(it)+0.003; 
     %     pow_point (it+1)=zp;% zero power point
     %     %time_point(it+2)=time_point(it)+8*hr;
     %     time_point(it+2)=time_point(it)+1;
     %     %pow_point (it+2)=zp;
     %     %Daytime
     %     it=length(time_point);
     %     %time_point(it+1)=time_point(it)+300;
     %     time_point(it+1)=time_point(it);
     %     %time_point(it+1)=time_point(it);
     %     pow_point (it+1)=1;
     %     %time_point(it+2)=time_point(it)+16*hr;
     %     time_point(it+2)=time_point(it)+2;
     %     pow_point (it+2)=1;
     % end
     % 
     % ref_old=pow_point(1); 
     % for it=1:nt % nt is length of the simulation time slots
     %     if(it>1)
     %         time(it)=time(it-1)+dt; % Simulation time step and putting it in  time slots
     %         ref_old=ref(it-1);
     %     end
     %     ref(it)=ref_old;
     %     i1=1;
     %     i2=1;
     %     for ii=1:length(time_point)-1
     %         if(time_point(ii)<=time(it) && time(it)<=time_point(ii+1)) %why this condition?
     %             i1=ii;
     %             i2=ii+1;
     %             frac1=(time_point(ii+1)-time(it))/(time_point(ii+1)-time_point(ii));
     %             frac2=1.0-frac1;
     %             ref(it)=frac1*pow_point(i1)+frac2*pow_point(i2);
     %             break
     %         end
     %     end
     % end
     % yref=(ref.*n0); 
    %% load following

    dt=1e-4;
    T = 200;
    n0=1;
    time=0:dt:T;
    nt=length(time);
    ref=zeros(nt,1);



    %time_point=[0	100	700	1000	1300	1600	1780	2080	2260	2560	2680	2980	3040	3340	3400	3700]
    %pow_point=[1	1	0.9	0.9	0.8	0.8	0.65	0.65	0.95	0.95	0.55	0.55	0.15	0.15	0.95	0.95]

    %time_point=[0 20 30 50 60 80 90 110 130 150];%work with only this time
    %time_point=[0 10 40 70 80 100 110 130 140 150];
    %time_point=[0 30 40 50 60 70 90 100 140 150];
    time_point=[0 20 30 50 60 80 90 110 130 200];
    %pow_point=[0.8 0.8 0.4 0.4 0.8 0.8 0.4 0.4 0.8 0.8];power Transient
    %pow_point=[1 1 0.4 0.4 1 1 0.4 0.4 1 1];
    %pow_point=[1 1 0.4 0.4 0.6 0.6 0.8 0.8 1 1];
    %pow_point=[0.9 0.9 0.3 0.3 0.9 0.9 0.3 0.3 0.9 0.9];
    pow_point=[0.7 0.7 0.4 0.4 0.4 0.4 0.6 0.6 0.7 0.7];
    %pow_point=[0.5 0.5 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1];

ref_old=pow_point(1);
for it=1:nt %to vectorize nt or make each point of nt
    if(it>1)
        time(it)=time(it-1)+dt;%to be able to make equal time step of the time
        ref_old=ref(it-1); %for each time slot 
    end
    ref(it)=ref_old;
    i1=1;
    i2=1;
    for ii=1:length(time_point)-1
        if(time_point(ii)<=time(it) && time(it)<=time_point(ii+1))
            i1=ii;
            i2=ii+1;
            frac1=(time_point(ii+1)-time(it))/(time_point(ii+1)-time_point(ii));
            frac2=1.0-frac1;
            ref(it)=frac1*pow_point(i1)+frac2*pow_point(i2);
            break
        end
    end

end
ref(:)=ref(:);
yref=(ref.*n0);
for s =1:nt
    setpoint=yref(s);
end
%setpoint = 1;       % Reference input (actual one to be supplied by the user)
Kp       = 3000;       % Proportional gain
Ki       = 1000;     % Integral gain
Kd       = 50;    % Derivative gain
dt       = 1e-4;    % ideally it should be the time elapsed between the current and previous time steps

%% Call ode45
tspan    = [0 20];
x0       = [0.121886811335360; 0.121886893341315; 0.121885271145068; 0.121886810283437; 0.0651359556098032; 0.0704553771394268; -0.0244616314872170; 23549641141.3791; 16604965156.7948; 0.000506372787486842; 0.00153255417325857; 863.999999999640; -0.000221924659921813; -0.000221924659921813];
[t, x]   = ode45(@(t, x) reactorDAE(t, x, pidController(t, x, setpoint, Kp, Ki, Kd, dt)), tspan, x0);

%% Plot result
figure(1)
plot(t, x(:,1)), grid on, xlabel('Time'), title('x_{1}')

%% Reactor dynamics
function [dx, algebraic] = reactorDAE(t, x, u)
    %% Parameters
    Sig_x   = 2.65e-22;
    yi      = 0.061;
    yx      = 0.002;
    lamda_x = 2.09e-5;
    lamda_I = 2.87e-5;
    Sum_f   = 0.3358;
    % nv    = 2.2e+3;           % average velocisty of neutrons (thermal)
    % nu    = 2.43;             % the total number of neutrons liberated per rx
    % Sigf  = 1/(gen*nv*nu);    % what is this 
    l       = 1.68e-3;
    beta    = 0.0065;
    beta_1  = 2.267510e-4;
    beta_2  = 1.22023e-3;
    beta_3  = 1.193061e-3;
    beta_4  = 2.785813e-3;
    beta_5  = 1.257395e-3;
    beta_6  = 5.226188e-4;
    Lamda_1 = 0.0124;
    Lamda_2 = 0.0369;
    Lamda_3 = 0.632;
    Lamda_4 = 3.044508e-1;
    Lamda_5 = 8.539933e-1;
    Lamda_6 = 2.868585;
    % Gr8   = -2350E-5/180*2; % All drums
    Gr4     = 1450E-5/180*2; % half
    % Gr2   = -660E-5/180*2; % 1/4
    Gr1     = 250E-5/180*2; % one
    cp_f    = 977;
    cp_m    = 1697;
    cp_c    = 5188.6;
    M_f     = 2002;
    M_m     = 11573;
    M_c     = 500;
    mu_f    = M_f*cp_f;
    mu_m    = M_m*cp_m;
    mu_c    = M_c*cp_c;
    f_f     = 0.96;
    P_0     = 20;
    Tf0     = 1105;
    Tm0     = 1087;
    T_in    = 864;
    T_out   = 1106;
    Tc0     = (T_in+T_out)/2;
    K_fm    = f_f*P_0/(Tf0-Tm0);
    K_mc    = P_0/(Tm0-Tc0);
    M_dot   = 1.75E+01;
    alpha_f = -2.875e-5;
    alpha_m = -3.696e-5;
    alpha_c = 0.0;
    X0      = 2.35496411413791e10;

    %% Declaration of state variables, x(i), where i = 1 to 14
    n_r     = x(1); 
    Cr1     = x(2); 
    Cr2     = x(3); 
    Cr3     = x(4); 
    Cr4     = x(5);
    Cr5     = x(6);
    Cr6     = x(7);
    X       = x(8); 
    I       = x(9); 
    Tf      = x(10); 
    Tm      = x(11); 
    Tc      = x(12); 
    Rho_d1  = x(13);
    Rho_d2  = x(14);
    
    G       = 3.2e-11;
    V       = 400*200;% i need to confirm the volume to be used 
    Pi      = P_0/(G*Sum_f*V);
    
    %% Define inputs (control signals) - here, for simplicity, assume constant
    % u1 = 0.01;          % Control signal 1 (e.g., PID output)
    % u2 = 0.005;         % Control signal 2 (e.g., temperature control)
    % u1 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u2 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u3 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u4 = Kp*error + Ki*error_integrated + Kd*(error - error_prev)/dt;
    % u1 = pid(kp,ki,kd);
    % u2 = pid(kp,ki,kd);
    % u3 = pid(kp,ki,kd);
    % u4 = pid(kp,ki,kd);
    
    %% The extra parameter u in reactorDAE(t, x, u) comes from the pidController()
    u1  = u;    % actuation thru channel x13, both PID controllers are the same 
    u4  = u;    % actuation thru channel x14
    
    %% ODEs
    dx      = zeros(14,1);
    rho     = Rho_d1 + Rho_d2 + alpha_f*(Tf - Tf0) + alpha_c*(Tc - Tc0) + alpha_m*(Tm - Tm0) - Sig_x*(X - X0)/Sum_f;
    
    %% kinetics equations with six-delayed neutron groups
    dx(1)   = (rho-beta)/l*n_r+beta_1/l*Cr1+beta_2/l*Cr2+beta_3/l*Cr3+beta_4/l*Cr4+beta_5/l*Cr5+beta_6/l*Cr6;
    dx(2)   = Lamda_1*n_r-Lamda_1*Cr1;
    dx(3)   = Lamda_2*n_r-Lamda_2*Cr2;
    dx(4)   = Lamda_3*n_r-Lamda_3*Cr3;
    dx(5)   = Lamda_4*n_r-Lamda_4*Cr3;
    dx(6)   = Lamda_5*n_r-Lamda_5*Cr3;
    dx(7)   = Lamda_6*n_r-Lamda_6*Cr3;
    
    %% Xenon and Iodine dynamics
    dx(8)   = yx*Sum_f*Pi+lamda_I*I-Sig_x*X*Pi-lamda_x*X;
    dx(9)   = yi*Sum_f*Pi-lamda_I*I;
    
    %% thermal–hydraulics model of the reactor core
    dx(10)  = f_f*P_0/mu_f*n_r-K_fm/mu_f*(Tf-Tc);
    dx(11)  = (1-f_f)*P_0/mu_m*n_r+(K_fm*(Tf-Tm)-K_mc*(Tm-Tc))/mu_m;
    dx(12)  = K_mc*(Tm-Tc)/mu_c-2*M_dot*cp_c*(Tc-T_in)/mu_c;
    dx(13)  = Gr1*u1;
    dx(14)  = Gr4*u4;
    
    %% Algebraic equations
    algebraic = [];
end

%% PID controller
function u = pidController(t, x, setpoint, Kp, Ki, Kd, dt)
    persistent integral_error prev_error prev_time;
    if isempty(integral_error)
        integral_error  = 0;
        prev_error      = 0;
        prev_time       = 0;
    end
    
    error               = setpoint - x(1);
    integral_error      = integral_error + error*(t - prev_time);
    derivative_error    = (error - prev_error)/dt;
    prev_error          = error;
    prev_time           = t;
    
    % user-supplied PID control equation
    u = Kp*error + Ki*integral_error + Kd*derivative_error;
end

[Sam Chak]

@Kamal, It's great to hear that the Reactor system can be stabilized using the PID controller code proposed in MATLAB. Regarding the over 300% overshoot issue, it would be wise to consult your Control Design team, as they possess significant knowledge in this research area.

However, I'm a bit confused. Previously, your Control Design team successfully designed a PID controller with the values Kp = 1, Ki = 0.1, and Kd = 0.01. You then implemented this controller using the PID block in Simulink, resulting in the Figure 2 outcome. It's quite perplexing that the MATLAB result turned out to be completely different from the Simulink result. I expected the MATLAB result to exhibit behavior similar to the Simulink result shown in Figure 2.

[Kamal]

@Sam Chak, the gains you mentioned above were just my own trial while trying my hand with MATLAB scripts. I used about 3000, 1000 abd 100 for the three gains respectively. I will send the simulink to you once I get hold of my system. I am in transit now. Thank you so much for your help.

[Kamal]

Please help check the concern I raised earlier about the ‘for loop’ I am trying to simulate repeating sequence input

[Sam Chak]

Hi @Kamal

The generation of the reference input signal seems to be a separate issue. Instead of treating it as a control problem, I recommend posting it as a new problem focusing on the MATLAB coding aspect. Users who excel at loop structures would likely provide faster responses. Additionally, please include the expected shape or form of the reference input signal in your post.

[Kamal]

@Sam Chak Alright I will do that. Thank you so much for your effort, time, patience. I really appreciate.

Answer 3

[Sam Chak]

Hi @Kamal

Below, you will notice a slight perceptible difference in the results between the Ideal PID controller and the PID Control Emulator when simulating a stably-damped mass-spring system. I should emphasize that the code for the PID Control Emulator, which utilizes the 'persistent' variables technique, only approximates the Ideal PID controller with limited accuracy.

The reason the PID Control Emulator cannot precisely replicate the Ideal PID controller's results lies in the fact that declaring the 'persistent' variables is merely a programming technique, without taking into consideration for the true mathematics governing the dynamics of the Ideal PID controller:

Ideal PID controller: ,

where e (error signal) is the difference between the desired (setpoint) and measured outputs.

Additionally, while the ode45 solver employs an adaptive timestep integration method, the 'integral_error' term in the PID Control Emulator relies on the less accurate forward Euler method for numerical integration, known for its first-order accuracy. Furthermore, as mentioned by @Torsten and @Jan in this thread, the concept of a "previous timestep" in ODE45's adaptive timestep execution holds no meaningful relevance.

For now, you can enhance the accuracy of the results from the PID Control Emulator by adjusting the value of the delta time 'dt'. Determine the delta time so that it evenly divides the simulation time by the number of time intervals. The number of time intervals equals the number of elements in the computed time vector minus one, which can be checked by examining the size of the time vector, 't'.

%% Parameters
setpoint = 1;                   % Reference input
Kp       = 0.75;                % Proportional gain
Ki       = 0.5;                 % Integral gain
Kd       = 0.125;               % Derivative gain
dt       = 10/1080;             % delta time (simulation time / number of time spacing)

%% Generate result from Ideal PID controller (requires Control System Toolbox)
Gp       = tf(1, [1, 2, 1]);    % plant
Gc       = pid(Kp, Ki, Kd);     % ideal PID controller
Gcl      = feedback(Gc*Gp, 1);  % closed-loop system
tt       = linspace(0, 10, 1001);
y        = step(Gcl, tt);       % step response

%% Generate result from PID Control Emulator using ode45
tspan    = [0 10];              % simulation time span
x0       = [0; 0];              % initial values of the state variables
[t, x]   = ode45(@(t, x) systemDyn(t, x, pidEmulator(t, x, setpoint, Kp, Ki, Kd, dt)), tspan, x0);
size_t   = size(t)

size_t = 1x2

1081 1

%% Plot results for comparison
plot(tt, y),      hold on
plot(t,  x(:,1)), grid on, 
legend('Ideal PID controller', 'PID control emulator')
xlabel('Time'), ylabel('x_{1}'), title('Comparison')

%% System Dynamics
function [dx, y] = systemDyn(t, x, u)
    % 'x' is the state variables
    % 'u' is the control input signal
    
    A   = [0,  1;   % state  matrix
          -1, -2];
    B   = [0;       % input  matrix
           1];
    C   = [1,  0];  % output matrix
    D   = [0];      % direct matrix
    
    % State-space model
    dx  = A*x + B*u;
    y   = C*x + D*u;
end

%% PID Control Emulator
function u = pidEmulator(t, x, setpoint, Kp, Ki, Kd, dt)
    persistent integral_error prev_error prev_time;
    if isempty(integral_error)
        integral_error  = 0;
        prev_error      = 0;
        prev_time       = 0;
    end
    
    error               = setpoint - x(1);                          % feedback loop
    integral_error      = integral_error + error*(t - prev_time);   % Euler method
    derivative_error    = (error - prev_error)/dt;                  % de/dt
    prev_error          = error;                                    % previous error value
    prev_time           = t;                                        % previous time value
    
    % PID control equation
    u = Kp*error + Ki*integral_error + Kd*derivative_error;
end

[Kamal]

@Sam Chak Thnak you so much for this wonderful thread, which has actually improved my knowledge. I understand your point here. I will adjust the dt as advised.

[Kamal]

@Sam Chak I simulated for 2000s and the initial oscilation was not captured. Is there anyway in the code to select desired decimation point of my data just like in the ToFile block of Simulink that I can select decimation point. Also, is there any means to save data generated to file instead of workspace or both with the code?

[Sam Chak]

Hi @Kamal, I'm not familiar with the process of selecting the decimation point in the To File block in Simulink. I suggest copying and pasting the MATLAB code and creating a new question to address this issue. By doing so, you can seek assistance from experts in ode solvers who can provide insights and help resolve the problem.