PID controller for a car suspension system not functioning.

Hello everyone, I am currently working on a project and am near the end of it. It is comparing the performance of a passive suspension system to that of a suspension system system with an LQR controller and another with a PID controller. So far the suspension system with the LQR controller and the Passive suspension system with no controller are working as intended. However, the PID controller seems to not be working at all. When I remove the PID controller in Simulink and treat that segment of the model like it is a Passive suspension system with no controller it yields no results at all either. There isn't any error popping up for me so I'm not sure as to where I may have made my error. Is there something wrong with my Matlab code or did I fumble with the Simulink model? Any help would be greatly appreciated!

Here is my code for the project so far (one of bode_plot and one for the model) along with the file for the Simulink model, my appologies for the length of it.

clc

clear all

ms = 1000;

mus = 200;

ks = 100000;

bs = 6000;

kus = 400000;

bus = 0;

sim_time = 0.01:0.04:10;

road_surface = zeros(2,length(sim_time));

A = [ 0 1 0 -1;

-ks/ms -bs/ms 0 bs/ms ;

0 0 0 1;

ks/mus bs/mus -kus/mus -(bs+bus)/mus];

B = [ 0 0;

0 1/ms;

-1 0;

bus/mus -1/mus];

C = [ eye(4)

-ks/ms -bs/ms 0 bs/ms];

D = [zeros(4,2)

0 1/ms];

%% PID Tuning

rank(ctrb(A,B));

Q = diag([1e10, 1e8, 1, 1]);

R = 0.01;

K = lqr( A, B(:,2), Q, R );

fprintf( 'Starting Simulation\n' )

tic;

model_sim = sim('Quarter_Car_Suspension_LQR_PID3');

elapsed_time_simulation = toc; % stop measuring time

fprintf( 'Simulation completed\n' )

fprintf( 'The total simulation time was %.2f seconds\n', elapsed_time_simulation)

% state space variable

state_passive = model_sim.state.state_passive;

state_LQR = model_sim.state.state_LQR;

state_PID = model_sim.state.state_PID;

figure('Name','Suspension travel'); hold on;

plot(state_passive.xs_xus)

plot(state_LQR.xs_xus)

plot(state_PID.xs_xus)

legend({'passive','LQR','PID'})

title('Suspension travel')

ylabel('x_s - x_{us} [m]')

xlabel('Time [s]')

xlim([0 10])

figure('Name','Sprung mass velocity'); hold on;

plot(state_passive.xs_d)

plot(state_LQR.xs_d)

plot(state_PID.xs_d)

legend({'passive','LQR','PID'})

title('Sprung mass velocity')

ylabel('dot{x_s} [m/s]')

xlabel('Time [s]')

xlim([0 10])

figure('Name','Tire deflection'); hold on;

plot(state_passive.xus_xr)

plot(state_LQR.xus_xr)

plot(state_PID.xus_xr)

legend({'passive','LQR','PID'})

title('Tire deflection')

ylabel('x_{us} - x_r [m/s]')

xlabel('Time [s]')

xlim([0 10])

figure('Name','Unsprung mass velocity'); hold on;

plot(state_passive.xus_d)

plot(state_LQR.xus_d)

plot(state_PID.xus_d)

legend({'passive','LQR','PID'})

title('Unsprung mass velocity')

xlabel('Time (s)')

ylabel('dot{x_{us}} [m/s]')

xlabel('Time [s]')

xlim([0 10])

% sprung mass acceleration

sprung_mass_a = model_sim.sprung_mass_acc;

figure('Name','Sprung mass acceleration'); hold on;

plot(sprung_mass_a.xs_dd_passive)

plot(sprung_mass_a.xs_dd_LQR)

plot(sprung_mass_a.xs_dd_PID)

legend({'passive','LQR','PID'})

title('Sprung mass acceleration')

ylabel('ddot{x_s} [m/s^2]')

xlabel('Time [s]')

xlim([0 10])

% sprung mass motion

sprung_mass_motion = model_sim.ms_motion;

figure('Name','Sprung mass motion'); hold on;

plot(sprung_mass_motion.xs_passive)

plot(sprung_mass_motion.xs_LQR)

plot(sprung_mass_motion.xs_PID)

legend({'passive','LQR','PID'})

title('Sprung mass motion')

ylabel('x_s [m]')

xlabel('Time [s]')

xlim([0 10])

%% Evaluation param

sp_passive = state_passive.xs_xus.Data;

sp_LQR = state_LQR.xs_xus.Data;

sp_PID = state_PID.xs_xus.Data;

msd_passive = state_passive.xs_d.Data;

msd_LQR = state_LQR.xs_d.Data;

msd_PID = state_PID.xs_d.Data;

td_passive = state_passive.xus_xr.Data;

td_LQR = state_LQR.xus_xr.Data;

td_PID = state_PID.xus_xr.Data;

musd_passive = state_passive.xus_d.Data;

musd_LQR = state_LQR.xus_d.Data;

musd_PID = state_PID.xus_d.Data;

msdd_passive = sprung_mass_a.xs_dd_passive.Data;

msdd_LQR = sprung_mass_a.xs_dd_LQR.Data;

msdd_PID = sprung_mass_a.xs_dd_PID.Data;

ms_motion_passive = sprung_mass_motion.xs_passive.Data;

ms_motion_LQR = sprung_mass_motion.xs_LQR.Data;

ms_motion_PID = sprung_mass_motion.xs_PID.Data;

Answer 1

[Sam Chak]

Hi @Bryan Shivram

I am not familiar with your MIMO quarter car suspension system, but I'd advise conducting a test in MATLAB before implementing the controller in Simulink. If Output 1 of the system is to be controlled, I'd likely mathematically design the controller for Input 2 in this manner, although I do not know the specific roles of Inputs 1 and 2 in your case. Logically, I'd suppose that one of the inputs must represent the road surface profile that will disturb the car suspension system.

%% parameters
ms  = 1000;
mus = 200;   
ks  = 100000; 
bs  = 6000;  
kus = 400000; 
bus = 0;     

%% system matrices
A   = [ 0        1        0       -1; 
       -ks/ms   -bs/ms    0        bs/ms ;
        0        0        0        1; 
        ks/mus   bs/mus  -kus/mus -(bs+bus)/mus];
B   = [ 0        0; 
        0        1/ms; 
       -1        0; 
        bus/mus -1/mus];
C   = [ eye(4)
       -ks/ms   -bs/ms    0        bs/ms];
D   = [zeros(4,2)
        0        1/ms];

%% MIMO system
sys = compreal(ss(A, B, C, D), 'o');
G   = tf(sys)

G =
 
  From input 1 to output...
                  -2000 s + 9.107e-12
   1:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
                    1.2e04 s + 2e05
   2:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
           -s^3 - 36 s^2 - 600 s - 2.636e-15
   3:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
              2000 s^2 + 1.2e04 s + 2e05
   4:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
            1.2e04 s^2 + 2e05 s - 8.601e-10
   5:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
  From input 2 to output...
              0.006 s^2 + 2.655e-16 s + 2
   1:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
       0.001 s^3 - 2.2e-17 s^2 + 2 s - 9.661e-15
   2:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
          -0.005 s^2 + 3.457e-17 s - 8.18e-16
   3:  -----------------------------------------
       s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
       -0.005 s^3 - 9.643e-17 s^2 - 5.978e-16 s + 1.853e-16
   4:  ----------------------------------------------------
            s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
       0.001 s^4 + 6.867e-18 s^3 + 2 s^2 + 3.509e-15 s - 2.348e-30
   5:  -----------------------------------------------------------
                s^4 + 36 s^3 + 2600 s^2 + 1.2e04 s + 2e05
 
Continuous-time transfer function.

%% PID controller (SISO)
ctrl            = pid(0, 2e5)   % a pure integral action seems to be good enough

ctrl =
 
        1 
  Ki * ---
        s 

  with Ki = 2e+05
 
Continuous-time I-only controller.

%% connect the SISO PID controller to the 2nd input of the MIMO plant
feedforwPath    = ctrl*sys; % the feedforward path
feedbackPath    = 1;        % the feedback path has a gain of 1
feedin          = 2;        % the feedback output is connected to input "2" of the plant
feedout         = 1;        % the output 1 of the plant G is connected to the feedback
closedLoopSys   = feedback(feedforwPath, feedbackPath, feedin, feedout, -1);

% simulation time vector
t   = 0:0.01:5;

% define a two-reference input signal
r1  = zeros(size(t));       % 1st ref input (no input)
r2  = ones(size(t));        % 2nd ref input (step input)
Ref = [r1', r2'];           % Create input matrix (2 columns for 2 inputs)

% Simulate the response of the closed-loop system
[y, tOut]       = lsim(closedLoopSys, Ref, t);

% Plot the responses
figure
plot(tOut, y(:,1)), grid on
xlabel('Time (s)')
ylabel('Output 1')
title('Response of MIMO Plant with a SISO PID Controller')

figure
tL  = tiledlayout(2, 2, 'TileSpacing', 'Compact');
out = 2:5;
for i = 1:numel(out)
    nexttile
    plot(tOut, y(:,out(i))), grid on
    xlim([0, 5])
    title("Output "+string(out(i)))
end
xlabel(tL, 'Time (s)')
ylabel(tL, 'Amplitude')
title(tL, 'Responses of other states')

[Sam Chak]

Hi @Bryan Shivram

By the way, the so-called LQR controller is not truly a controller in the mathematical sense. The LQR is a computational algorithm designed to determine the stabilizing gain matrix for the full-state feedback controller, expressed as .

The effectiveness of this approach compared to PID controllers lies in the assumption that a full-state feedback controller can measure all output states of the system accurately, a condition that is rarely met in real-world applications.

In contrast, the PID controller is effectively a simplistic 2nd-order system. For plants of order greater than two, we can only hope for favorable outcomes when employing a PID controller 'alone'. Additionally, your MIMO quarter car suspension system is a 4th-order system, as indicated by the term in the denominators of the MIMO transfer function . Thus, controlling a 4th-order system with a 2nd-order controller may not yield the desired performances that we hope to see.