Steady-State Initialization of Quarter-Car Simulink Model

R2026b

This example uses:

Open Live Script

This example demonstrates how to find and apply steady-state initialization to a Simulink® model representing a quarter-car, using Simscape™ Multibody™ MATLAB® classes. It shows how to eliminate initial transients by initializing the model at equilibrium.

The example shows how to:

  • Extract a CompiledMultibody object from an existing Simulink model

  • Find the steady-state equilibrium position using optimization

  • Create an operating point for model initialization

  • Apply the operating point to eliminate startup transients

  • Compare simulation results with and without proper initialization

For best interaction with this example, run it section by section.

Load and Compile the Simulink Model

First, load the existing quarter-car Simulink model and compile it to access the CompiledMultibody object:

% To simplify the process and avoid repeatedly typing the namespace name
% for the classes, you can use the import function.
import simscape.Value simscape.op.* simscape.multibody.*;
addpath(fullfile('SteadyStateInitializationOfQuarterCarSimulinkModelSupport','Scripts'));

% Load the Simulink model
modelName = "SteadyStateInitializationOfQuarterCarSimulinkModel";
load_system(modelName);

% Compile to CompiledMultibody object
cmb = simscape.multibody.compileBlockDiagram(modelName);

Find Steady-State Equilibrium

Find the configuration where all accelerations are zero, representing the steady-state position under gravity:

% Define equilibrium equations (accelerations = 0)
equilibriumEqs = @(x) computeAccelerations(cmb, x);

% Initial guess - assume some compression in both joints
x0 = [-0.05; -0.02];  % [suspension; tire] displacement in meters

% Solver options
options = optimoptions('fsolve', ...
    'Display', 'iter', ...
    'FunctionTolerance', 1e-10, ...
    'StepTolerance', 1e-10);

% Find equilibrium
[equilibrium, fval, exitflag] = fsolve(equilibriumEqs, x0, options);
                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0          3             3861.91                        2.82e+05              1
     1          6         2.77122e-15       0.119177         0.000219              1
     2          9         4.10208e-29    5.24078e-10         2.86e-11              1

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

fprintf('\n=== Steady-State Equilibrium Found ===\n');
=== Steady-State Equilibrium Found ===

fprintf('Suspension position: %.4f m (%.1f mm compression)\n', ...
    equilibrium(1), -equilibrium(1)*1000);
Suspension position: -0.1692 m (169.2 mm compression)

fprintf('Tire position: %.4f m (%.1f mm compression)\n', ...
    equilibrium(2), -equilibrium(2)*1000);
Tire position: -0.0180 m (18.0 mm compression)

fprintf('Residual accelerations: [%.2e, %.2e] m/s^2\n', fval(1), fval(2));
Residual accelerations: [5.33e-15, -3.55e-15] m/s^2

Create Operating Point for Initialization

Create an operating point with the equilibrium positions and zero velocities:

% Create operating point
steadyStateOP = OperatingPoint;

% Set positions at equilibrium
steadyStateOP("Suspension Joint/Pz/p") = Target(equilibrium(1), "m", "High");
steadyStateOP("Tire Joint/Pz/p") = Target(equilibrium(2), "m", "High");

% Set velocities to zero
steadyStateOP("Suspension Joint/Pz/v") = Target(0, "m/s", "High");
steadyStateOP("Tire Joint/Pz/v") = Target(0, "m/s", "High");
disp(steadyStateOP);
  OperatingPoint with children:

  OperatingPoints:

   ChildId             Size
   __________________  ____

   'Suspension Joint'   1x1
   'Tire Joint'         1x1

% Save operating point to base workspace for Simulink
assignin('base', 'quarterCarSteadyStateOP', steadyStateOP);

Simulate Without Initialization

First, run the model without steady-state initialization to see the transient behavior:

% Configure simulation without initialization
set_param(modelName, 'LoadInitialState', 'off');
simTime = 5;  % seconds
set_param(modelName, 'StopTime', num2str(simTime));

% Run simulation
simOut_noInit = sim(modelName);

% Extract results
logsOut_noInit = simOut_noInit.simlog;
time_noInit = simOut_noInit.tout;
suspension_pos_noInit = logsOut_noInit.get("Suspension Joint/Pz/p").series.values;
tire_pos_noInit = logsOut_noInit.get("Tire Joint/Pz/p").series.values;

Apply Steady-State Initialization

Configure the model to use the steady-state operating point:

% Set operating point in model configuration
set_param(modelName,'SimscapeUseOperatingPoints','on');
set_param(modelName,'SimscapeOperatingPoint','quarterCarSteadyStateOP');

% Alternatively, the operating point could be set using the Configuration
% Parameters of the QuarterCar Simulink model, as shown in the following image:

Simulate with Steady-State Initialization

Now run the model with proper initialization:

% Run simulation with initialization
simOut_withInit = sim(modelName);

% Extract results
logsOut_withInit = simOut_withInit.simlog;
time_withInit = simOut_withInit.tout;
suspension_pos_withInit = logsOut_withInit.get("Suspension Joint/Pz/p").series.values;
tire_pos_withInit = logsOut_withInit.get("Tire Joint/Pz/p").series.values;

Compare Results

Visualize the difference between initialized and non-initialized simulations:

% Create comparison plots
plotInitializationComparison(time_noInit, suspension_pos_noInit, tire_pos_noInit, ...
    time_withInit, suspension_pos_withInit, tire_pos_withInit, ...
    equilibrium);

Figure contains 3 axes objects and another object of type subplottext. Axes object 1 with title Suspension Response, xlabel Time (s), ylabel Suspension Position (mm) contains 3 objects of type line, constantline. These objects represent No Init (Transient), With Init (Steady), Equilibrium. Axes object 2 with title Tire Response, xlabel Time (s), ylabel Tire Position (mm) contains 3 objects of type line, constantline. These objects represent No Init (Transient), With Init (Steady), Equilibrium. Axes object 3 with title Transient Error (Without Initialization), xlabel Time (s), ylabel Distance from Equilibrium (mm) contains 2 objects of type line. These objects represent Suspension, Tire.

% Calculate settling metrics
fprintf('\n=== Transient Analysis ===\n');
=== Transient Analysis ===

[suspension_settling_time, tire_settling_time] = analyzeTransients(...
    time_noInit, suspension_pos_noInit, tire_pos_noInit, equilibrium);

fprintf('Without initialization:\n');
Without initialization:

fprintf('  Suspension settling time: %.3f s\n', suspension_settling_time);
  Suspension settling time: 2.723 s

fprintf('  Tire settling time: %.3f s\n', tire_settling_time);
  Tire settling time: 2.297 s

fprintf('\nWith initialization:\n');
With initialization:

fprintf('  Settling time: 0 s (starts at equilibrium)\n');
  Settling time: 0 s (starts at equilibrium)

Summary

This example demonstrated the power of combining Simulink models with the Simscape Multibody MATLAB classes for steady-state initialization. Key benefits include:

  1. Immediate steady state: No waiting for transients to settle

  2. Consistent initial conditions: Reproducible simulations

  3. Efficient workflow: Find equilibrium once, use many times

  4. Better for control: Start controller design from steady state

  5. Automated process: Can be scripted for parameter variations

The compileBlockDiagram function bridges the gap between Simulink graphical models and programmatic analysis, enabling powerful initialization and analysis workflows.

rmpath(fullfile('SteadyStateInitializationOfQuarterCarSimulinkModelSupport','Scripts'));

% Close the system
close_system(modelName, 0);

See Also