# Simulation of control system with only matlab script withou simulink

59 views (last 30 days)

Show older comments

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 Comments

### Accepted Answer

Sam Chak
on 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 Comments

Sam Chak
on 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.

Figure 1: Reactor's dynamic system model in Simulink

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

### More Answers (2)

Sam Chak
on 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 Comments

Sam Chak
on 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
on 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 Comments

Sam Chak
on 12 Mar 2024

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!