Using ode45 to develop a PID controller to reduce transverse vibrations in a rotating beam

i have the following code but Im struggling to add a PID tuner.

% Angular velocity (specify your desired value)
w   = 10;

% Given parameters
L   = 16;           % Length of shaft
d   = 0.5;          % Diameter of shaft
p   = 7900;         % Density of steel
E   = 209e9;        % Young's modulus of material
dr  = 0.1;          % Damping ratio
ec  = d/2;          % Eccentricity
G   = 75e9;         % Shear Modulus of stainless steel

% Calculated parameters
A   = d^2;          % Cross-sectional area of shaft
V   = A * L;        % Volume of shaft
m   = p * V;        % Mass of shaft
I   = d^4/12;       % Moment of inertia of beam
k   = G*I/L;        % Stiffness of shaft
c2  = E*I/(p*A);    % Corrected damping coefficient
I2  = pi^2/2;
I1  = 1/2;
I4  = -pi^2/2;
Po  = (m*ec*w^2)/(p*A);

% Constants (same as before)
a   = 2*sqrt(m*k);      % Coefficient of z'(t)
b   = c2/L^4*I2/I1;     % Coefficient of z(t)
c   = 1/L^2* I4/I1;     % Coefficient of Q
d_coefficient = 1/I1 * w^2 * sin(1/2) * Po; % Coefficient of w^2*sin(wt)

% Desired value (setpoint)
desired_value = 0;

% Time span for simulation
tspan = [0, 10];

% Initialize variables for PID gains
Kp = 1;     % Proportional gain (adjust as needed)
Ki = 0.1;   % Integral gain (adjust as needed)
Kd = 0.01;  % Derivative gain (adjust as needed)

% Initialize variables for integral control
integral_term = 0;
prev_error = 0;

% Define the transfer function of the rotating system
sys = tf([1], [m, a, b, c])

sys =
 
                      1
  -----------------------------------------
  31600 s^3 + 1.757e06 s^2 + 83 s - 0.03855
 
Continuous-time transfer function.

% Define the differential equation function with PID controller
ode_pid = @(t, Z) [Z(2); 
               - a*Z(1) - b*Z(1) + c*(Kp*(desired_value - Z(1)) + Ki*integral_term + Kd*(desired_value - Z(1)) - prev_error) + d_coefficient*exp(1i*w*t)];

% Initial conditions [z(0), z'(0)]
initial_conditions = [0, 0];

% Solve the differential equation with PID controller
[t, Z] = ode45(ode_pid, tspan, initial_conditions);

% Extracting real part of position from solution
real_part = real(Z(:, 1));

% Initialize array to store Q values
Q_values = zeros(size(t));

% Calculate Q values
for i = 1:length(t)
    Q_values(i) = c * (Kp*(desired_value - real_part(i)) + Ki*integral_term + Kd*(desired_value - real_part(i)) - prev_error);
end

% Plotting
figure;

subplot(2,1,1);
plot(t, real_part, 'b');
xlabel('Time');
ylabel('Real Part of Position (z)');
title(['Real Part of Position (z) vs Time for w = ', num2str(w)]);
grid on;

subplot(2,1,2);
plot(t, Q_values, 'r');
xlabel('Time');
ylabel('Control Signal (Q)');
title('Control Signal (Q) vs Time');
grid on;

Any help would be appreciated

[Sam Chak]

The transfer function 'sys' in your code corresponds to a 3rd-order system. However, the ODE system described in the anonymous function 'ode_pid()' represents a 2nd-order system. The formulation for the ODE of the system can be expressed as follows:

,

To convert it into state-space form, we have:

% Define the transfer function of the rotating beam system
sys = tf([1], [m, a, b, c])

sys =
 
                      1
  -----------------------------------------
  31600 s^3 + 1.757e06 s^2 + 83 s - 0.03855
 
Continuous-time transfer function.

% Checking the order of the system
systemOrder = order(sys)

systemOrder = 3

Therefore, I suggest modifying the 'ode_pid()' function as follows, incorporating the PID control function call:

% PID controller function
pidCtrl = @(t, Z) ...;

% ODEs of rotating beam system with PID controller function call
ode_pid = @(t, Z) [Z(2); 
                   Z(3); 
             - c/m*Z(1) - b/m*Z(2) - a/m*Z(3) + 1/m*pidCtrl(t, Z)];

Answer 1

[Sam Chak]

Hi @Samuel Adshead

With this approach, you can utilize ode45 solver to simulate a PID-controlled rotating beam and plot the step response.

%% Angular velocity (specify your desired value)
w   = 10;

%% Given parameters
L   = 16;           % Length of shaft
d   = 0.5;          % Diameter of shaft
p   = 7900;         % Density of steel
E   = 209e9;        % Young's modulus of material
dr  = 0.1;          % Damping ratio
ec  = d/2;          % Eccentricity
G   = 75e9;         % Shear Modulus of stainless steel

%% Calculated parameters
A   = d^2;              % Cross-sectional area of shaft
V   = A * L;            % Volume of shaft
m   = p * V;            % Mass of shaft
I   = d^4/12;           % Moment of inertia of beam
k   = G*I/L;            % Stiffness of shaft
c2  = E*I/(p*A);        % Corrected damping coefficient
I2  = pi^2/2;
I1  = 1/2;
I4  = -pi^2/2;
Po  = (m*ec*w^2)/(p*A);

%% Constants (same as before)
a   = 2*sqrt(m*k);      % Coefficient of z'(t)
b   = c2/L^4*I2/I1;     % Coefficient of z(t)
c   = 1/L^2* I4/I1;     % Coefficient of Q

%% Transfer function of the rotating system
Gp  = tf([1], [m, a, b, c]);
Gp  = minreal(Gp)

Gp =
 
                 3.165e-05
  ---------------------------------------
  s^3 + 55.59 s^2 + 0.002627 s - 1.22e-06
 
Continuous-time transfer function.

%% PID Controller
Gc  = pidtune(Gp, 'PIDF', w)

Gc =
 
             1            s    
  Kp + Ki * --- + Kd * --------
             s          Tf*s+1 

  with Kp = 3.19e+07, Ki = 1.44e+07, Kd = 1.77e+07, Tf = 0.000875
 
Continuous-time PIDF controller in parallel form.

%% Closed-loop system
Gcl = feedback(Gc*Gp, 1);
sys = ss(Gcl)

sys =
 
  A = 
           x1      x2      x3      x4      x5
   x1   -1198  -248.2  -39.08  -8.811   -1.98
   x2     256       0       0       0       0
   x3       0      64       0       0       0
   x4       0       0       8       0       0
   x5       0       0       0       2       0
 
  B = 
       u1
   x1   8
   x2   0
   x3   0
   x4   0
   x5   0
 
  C = 
           x1      x2      x3      x4      x5
   y1       0       0   4.885   1.101  0.2475
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

%% State-space system
function [dxdt, y] = ode(t, x, A, B, C, D);
    u    = 1;           % unit step input
    dxdt = A*x + B*u;
    y    = C*x + D*u;
end

%% Call ode45 solver
tspan   = [0 5];
x0      = zeros(5, 1);
[t, x]  = ode45(@(t, x) ode(t, x, sys.A, sys.B, sys.C, sys.D), tspan, x0);
[~, y]  = ode(t', x', sys.A, sys.B, sys.C, sys.D);

%% Plot Step response
plot(t, y), grid on
title('Step response under PID controller')
xlabel('Time')
ylabel('Amplitude')