ADMM-Based MPC Control for Quarter-Car Suspension Dynamics

Quarter-Car Suspension Model

The Quarter-Car Suspension Model is a two degrees of freedom system that simplifies the vehicle suspension system by focusing on the vertical dynamics of one wheel and a quarter of the car body mass. This model comprises: two masses: the sprung mass $$m_s$$, which corresponds to the vehicle body, and the unsprung mass $$m_u$$, which represents the wheel assembly. The spring and damper of the suspension, with constants $$k_s$$ and $$c_s$$ respectively, link the two masses, while the tire is modeled as a spring with constant $$k_t$$. Active suspension forces are introduced through the control input $$u(t)$$, and the road effect is modeled by the excitation $$w(t)$$. The equations governing the system dynamics are:

$\{\begin{array}{l}{\mathit{m}}_{\mathit{s}}{\ddot{\mathit{x}}}_{\mathit{s}}=\mathit{f}\\ {\mathit{m}}_{\mathit{u}}{\ddot{\mathit{x}}}_{\mathit{u}}=-\mathit{f}-{\mathit{k}}_{\mathit{t}}\left({\mathit{z}}_{\mathit{u}}-{\mathit{z}}_{\mathit{r}}\right)\end{array}$, where$$f=-k_s(z_s-z_u)-b_s({\underset{}{\overset{\dot{}}{z}}}_s-{\underset{}{\overset{\dot{}}{z}}}_u)+u$$

Defining a state vector as:

$\mathit{z}\left(\mathit{t}\right)=\left[\begin{array}{c}\mathbf{x}\\ \dot{\mathbf{x}}\end{array}\right]$, where $\mathbf{x}=\left[\begin{array}{c}{\mathit{x}}_{\mathit{s}}\\ {\mathit{x}}_{\mathit{u}}\end{array}\right]$and $\dot{\mathbf{x}}=\left[\begin{array}{c}{\dot{\mathit{x}}}_{\mathit{s}}\\ {\dot{\mathit{x}}}_{\mathit{u}}\end{array}\right]$

The Quarter-Car Suspension Model is implemented using a linear state-space system.

% Physical parameters
ms = 300;    % kg
mu = 60;     % kg
bs = 1000;   % N/m/s
ks = 16000 ; % N/m
kt = 190000; % N/m

% Define the mass matrix M.
M = [ms 0; 0 mu];
% Define the stiffness matrix F.
F = [-ks ks; ks -(ks + kt)];
% Define the damping matrix B.
B = [-bs bs; bs -bs];
% Define the input matrix q.
q = [1; -1];
% Define the disturbance matrix d.
d = [0; kt];

% Define the State Space Matrices.
A = [zeros(2, 2), eye(2, 2); M\F, M\B];
B = [zeros(2,2); M\q M\d];
C = [1 -1 0 0];
D = 0;
plant = ss(A,B,C,D);

Road Disturbance Profile

For this example, the road disturbance is modeled as a sinusoidal bump with a frequency of $$f$$ Hz. The disturbance is active only between $$2/f$$ and $$3/f$$ seconds to simulate a single bump encounter by the vehicle. The amplitude of the disturbance is scaled by a factor to represent the severity of the bump.

% Define the frequency of the road disturbance.
f = 4; % Hz
Tstop = 2;

% Define the sampling time.
Ts = 0.01;

% Create a time vector from 0 to Tstop.
% with increments of Ts.
T = (0:Ts:Tstop)';

% Define the road disturbance profile.
v = 0.025*(1-cos(2*pi*f*T));

% Activate the disturbance only between 2/f and 3/f seconds.
v(T < 2*1/f) = 0;
v(T > 3*1/f) = 0;

Simulate the suspension system without the active control input to provide a baseline response of the system when subjected to the defined road disturbance.

% Assign the road disturbance profile 
% as a column vector [0 v] for simulation.
Unn = zeros(length(T),2);
Unn(:,2) = v;

% Simulate the uncontrolled system response 
% using the lsim function.
[Yuc,Tuc,Xuc] = lsim(plant,Unn,T);

Plot the impact of the road disturbance on the quarter-car suspension system along with the suspension deflection over time.

figure;
title('Response of Uncontrolled Suspension to Road Disturbance');
plot(T, v, 'DisplayName', 'Road Disturbance'); hold on;
plot(Tuc, Yuc, 'DisplayName', 'Suspension Deflection'); hold off;
legend;
ylabel('Meters');
xlabel('Time (seconds)');

Design Model Predictive Controller

By default all input signals of the plant state-space model are manipulated variables and all outputs are measured outputs. Define the second input as a measured disturbances correspondent to the road's effect.

plant = setmpcsignals(plant,'MD',2);

-->Assuming unspecified input signals are manipulated variables.

Create an MPC controller for a prediction horizon of 40 steps, control horizon of 4 steps and a sample time of 0.01 seconds.

Ts = 1e-2;
p = 40;
m = 4;
mpcobj = mpc(plant,Ts,p,m);

-->"Weights.ManipulatedVariables" is empty. Assuming default 0.00000.
-->"Weights.ManipulatedVariablesRate" is empty. Assuming default 0.10000.
-->"Weights.OutputVariables" is empty. Assuming default 1.00000.

mpcobj.OutputVariables(1).Min = -.05;
mpcobj.OutputVariables(1).Max = .05;

Set MPC Controller parameters to define the scale factors and weights. The scale factor is used to normalize the output variables, and the weights are used to specify the relative importance of each variable in the cost function. In this case, more emphasis is placed on the control input by setting a higher weight for the manipulated variable rate.

mpcobj.OutputVariables(1).ScaleFactor = 1e-5;
mpcobj.Weights.ManipulatedVariables = 1;
mpcobj.Weights.ManipulatedVariablesRate = 10;
mpcobj.Weights.OutputVariables = 5;

Review the MPC controller design to ensure it is configured as intended.

review(mpcobj)

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.
Test report has been saved to:
 C:\Users\gcampa\AppData\Local\Temp\index.html

% Define the reference signal 'r' as a vector of zeros, 
% which represents the desired output 
% (i.e., no suspension deflection).
r = zeros(length(T),1);

% Define the road disturbance 'v' for the simulation, 
% scaled to represent the severity of the bump 
% and active only within a specific time window.
v = 0.025 * (1 - cos(2 * pi * f * T));
v(T < 2 * 1/f) = 0;
v(T > 3 * 1/f) = 0;

% Set simulation options for the MPC controller. 
% The 'MDLookAhead' option is turned off, 
% indicating that the controller does not 
% anticipate future measured disturbances.
params = mpcsimopt();
params.MDLookAhead = 'off';

Setting the Optimization Solver

The 'active-set' algorithm is the default optimization method for MPC problems that contain only continuous manipulated variables. When you inspect the Optimizer property of an mpc object, you will typically see a structure similar to the following:

mpcobj.Optimizer

ans = struct with fields:
         OptimizationType: 'QP'
                   Solver: 'active-set'
            SolverOptions: [1×1 mpc.solver.options.qp.ActiveSet]
             MinOutputECR: 0
    UseSuboptimalSolution: 0
             CustomSolver: 0
      CustomSolverCodeGen: 0

Note that since the manipulated variable is continuous, the optimization type that needs to be solved is a Quadratic Programming (QP) problem. If you have finite (discrete) control sets, the problems become Mixed-Integer Quadratic Programming (MIQP) problems. For more information, see Finite Control Set MPC.

For this example, set the MPC object to use the ADMM solver:

mpcobj.Optimizer.Solver = "admm";

By changing the solver to admm, the MPC controller is now configured to use the Alternating Direction Method of Multipliers (ADMM), which is well-suited for large-scale and distributed optimization problems.

To modify the solver parameters and tailor the optimization process to specific requirements, you can access and adjust the settings through the SolverOptions property of the Optimizer attribute within the MPC object. This property allows you to fine-tune various aspects of the solver's behavior, such as convergence tolerances and iteration limit options.

mpcobj.Optimizer.SolverOptions

ans = 
  ADMM with properties:

      AbsoluteTolerance: 1.0000e-04
      RelativeTolerance: 1.0000e-04
          MaxIterations: 1000
    RelaxationParameter: 1
       PenaltyParameter: 1.6000
               StepSize: 1.0000e-06

[Yad, Tad, Uad, XPad, XCad] = sim(mpcobj,Tstop/Ts+1,r,v,params);

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.

Simulate the system using the 'active-set' optimization algorithm. This algorithm is one of the options available for solving the quadratic programming problem at each control interval.

mpcobj.Optimizer.Solver = 'active-set';
[Yas, Tas, Uas, XPas, XCas] = sim(mpcobj, Tstop/Ts + 1, r, v, params);

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.

Simulate the system using the 'interior-point' optimization algorithm. This algorithm is another option for solving the optimization problem and may offer different performance characteristics compared to the 'active-set' algorithm.

mpcobj.Optimizer.Solver = 'interior-point';
[Yip, Tip, Uip, XPip, XCip] = sim(mpcobj, Tstop/Ts + 1, r, v, params);

-->Converting model to discrete time.
-->Assuming output disturbance added to measured output #1 is integrated white noise.
-->"Model.Noise" is empty. Assuming white noise on each measured output.

Plot results.

TT = {Tas, Tip, Tad, Tuc};
XX = {Yas, Yip, Yad, Yuc};
UU = {Uas, Uip, Uad, Tuc*0};

plotResults(TT,XX,UU)
subplot(211); 
legend('active-set','interior-point','admm','Uncontrolled')

function plotResults(TT,XX,UU)
% Plot the results of the simulations to compare 
% the performance of the two optimization algorithms. 
% The plots will show the system output (suspension deflection) 
% and control input (active suspension force) over time.
subplot(211); 
for i = 1:length(TT)
    plot(TT{i}, XX{i})
    hold on;
end
ylabel('Meters')
xlabel('Time, s')

subplot(212); 
for i = 1:length(TT)
    plot(TT{i}, UU{i})
    hold on;
end
ylabel('N')
xlabel('Time, s')
end