Simulation of control system with only matlab script withou simulink
58 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposta accettata
Sam Chak
il 4 Mar 2024
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 Commenti
Sam Chak
il 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1639736/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1639741/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1639746/image.png)
Figure 1: Reactor's dynamic system model in Simulink
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1639751/image.png)
Figure 2: Result of x1 under PID control simulated in Simulink.
Più risposte (2)
Sam Chak
il 11 Mar 2024
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 Commenti
Sam Chak
il 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.
Sam Chak
il 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:
,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1640446/image.png)
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 Commenti
Sam Chak
il 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.
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!