Simulation of control system with only matlab script withou simulink
Afficher commentaires plus anciens
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');
3 commentaires
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)
Réponse acceptée
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
22 commentaires
Kamal
le 4 Mar 2024
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
le 4 Mar 2024
Kamal
le 4 Mar 2024
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
le 4 Mar 2024
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.
) 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
Sam Chak
le 5 Mar 2024
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.
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.
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.
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
le 5 Mar 2024
Kamal
le 8 Mar 2024
Kamal
le 9 Mar 2024
Modifié(e) : Walter Roberson
le 9 Mar 2024
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
%% --- 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
le 9 Mar 2024
Sam Chak
le 10 Mar 2024
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.
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!
Sam Chak
le 10 Mar 2024
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.
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
le 10 Mar 2024
Plus de réponses (2)
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);
%% 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
7 commentaires
@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
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
le 11 Mar 2024
@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
le 11 Mar 2024
Kamal
le 11 Mar 2024
Sam Chak
le 11 Mar 2024
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
le 12 Mar 2024
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)
%% 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
3 commentaires
Kamal
le 12 Mar 2024
Kamal
le 12 Mar 2024
@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
le 12 Mar 2024
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.
Catégories
En savoir plus sur Programming dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
