Cylinder Design Parameter Estimation

Open Model

This example shows how to parameterize and test a tandem primary cylinder starting from a manufacturer datasheet. This example uses optimization to determine remaining unknown parameters, given the numerical data extracted from the datasheet. After the model is simulated, it compares the resulting push rod force - brake pressure relationship curve to the curve provided on the manufacturer datasheet. Understanding the behavior of the tandem primary cylinder is important when selecting other braking system components.

The live script requires the Optimization Toolbox™ and Global Optimization Toolbox to be installed and licensed to run.

Model

The following figure shows the tandem primary cylinder test harness. The model uses spring-loaded accumulators for loading in the $1^{st} $ and the $2^{nd} $ circuits in place of disk or drum brakes. The input to the model is the push rod force coming from the lever mechanism or the brake booster. The output of the model is the pressures achieved in the brake circuits 1 and 2 respectively.

open_system('CylinderDesignParameterEstimation')

Manufacturer Datasheet Data

The supplier has provided the following data in the datasheet.

  1. tpcPistonDiameter : first and second circuit piston diameter.

  2. tpcPistonStroke1 : first circuit piston stroke.

  3. tpcPistonStroke2 : second circuit piston stroke.

  4. tpcTotalStroke : overall stroke.

  5. tpcCircuit1Disp : first circuit displacement.

  6. tpcCircuit2Disp : second circuit displacement.

  7. tpcTotalDisp : net displacement.

  8. _tpcMaxPressure _ : max pressure.

Pressure Development Characteristics

The supplier datasheet gives the following function diagram for the push rod force-brake pressure curve for the tandem primary cylinder which should be followed by the brake circuit 1 and 2 as per supplier catalog. It is required for a tandem primary cylinder to have same brake pressures created in brake circuits for uniform braking.

Equations of Motion (EoM) of System

The following figure shows the tandem primary cylinder. Mass 1 and mass 2 use the coordinate frames shown below. This description only includes the mechanical EoM of the system. The fluid dynamics equations are not shown.

$m_1\frac{d^2x_1}{dt^2}+c_1(\frac{dx_1}{dt}-\frac{dx_2}{dt})+k_1(x_1-x_2)=F_p -p_1 A -Preload_1+R_1$, when $x_1=0$

$m_2\frac{d^2x_2}{dt^2}+c_1(\frac{dx_2}{dt}-\frac{dx_1}{dt})+c_2\frac{dx_2}{dt}+k_1(x_2-x_1)+k_2 x_2 = p_1 A-p_2 A +Preload_1-Preload_2+R_2$, when $x_2=0$

$m_1\frac{d^2x_1}{dt^2}+c_1(\frac{dx_1}{dt}-\frac{dx_2}{dt})+k_1(x_1-x_2)=F_p -p_1 A -Preload_1$, when $x_1>0$

$m_2\frac{d^2x_2}{dt^2}+c_1(\frac{dx_2}{dt}-\frac{dx_1}{dt})+c_2\frac{dx_2}{dt}+k_1(x_2-x_1)+k_2 x_2 = p_1 A-p_2 A +Preload_1-Preload_2$, when $x_2>0$

where,

$x_1 $ and $x_2 $ are positions of piston mass 1 and mass 2 from the left side assumed hard stop.

$m_1 $ and $m_2 $ are masses of piston 1 and 2 respectively.

$c_1 $ and $c_2 $ are damping coefficients of piston 1 and 2 respectively.

$k_1 $ and $k_2 $ are stiffness coefficients of left (Spring 1) and right (spring 2) respectively.

$F_p $ is force applied on push rod of primary cylinder.

$P_1 $ and $P_2 $ are pressures in brake circuit 1 and 2 respectively.

$R_1 $ and $R_2 $ are reaction forces on mass 1 and 2 respectively when both masses are at the leftmost position.

$Preload_1 $ and $Preload_2 $ are preloads on left and right side springs respectively.

Parameter Estimation Scheme

The following steps outline how to estimate unknown feasible parameters of the tandem primary cylinder.

1. Identify Parameters to be Estimated

Some parameters, which are required for component model simulation, are not provided in the manufacturer's datasheet. Therefore, in order to validate the model with supplier function diagram it is necessary to estimate these parameters. The supplier function diagram gives quasi- steady state behavior of the tandem primary cylinder.

Where,

Data1, Data2, and Data3 are 3 points on the manufacturer function diagram.

Data1 point represents the leftmost points of the pistons in their respective coordinate frames.

Data3 point represent the right most points of the pistons in their respective coordinate frame, which are given in the datasheet.

Data2 point represent pistons movement at which brake circuits 1 and 2 get disconnected from the brake oil reservoir by the closing of the compensating ports. Movement beyond Data2 point in the right direction results in pressure building in the brake circuits. These points are not given in the datasheet and have to be estimated.

Left and right spring stiffness coefficients and pre-compressions are unknown.

Therefore there are 6 unknown parameters.

  1. tpcPosOpen1, which is piston 1 movement up to which compensating port remains connected with the brake circuit 1.

  2. tpcPosOpen2, which is piston 2 movement up to which compensating port remains connected with the brake circuit 2.

  3. tpcStiffness1Coeff, which is first spring stiffness coefficient

  4. tpcStiffness2Coeff, which is second spring stiffness coefficient

  5. tpcSpring1Comp, which is first spring pre-compression displacement

  6. tpcSpring2Comp, which is second spring pre-compression displacement. This parameter can be computed from the constraint that the Preload1 and Preload2 are equal at rest in the steady state of the system.

The datasheet does not provide the step test data, therefore the overshoot, rise time, settling time, and recovery time are not known. Therefore, the inertia and damping coefficient can't be estimated.

2. Provide Lower, Upper Bound, Scaling Factor and Initialization Point for the Parameter Estimation through Optimization 

% tpcPosOpen1Bound        = [0.5 1.5]*1e-3; % (m) Lower and upper bound for compensating port 1 parameter
% tpcPosOpen2Bound        = [0.5 1.5]*1e-3; % (m) Lower and upper bound for compensating port 2 parameter
% tpcStiffness1CoeffBound = [1e4 6e4];      % (N/m^2) Lower and upper bound for first spring stiffness coefficient
% tpcStiffness2CoeffBound = [1e4 6e4];      % (N/m^2) Lower and upper bound for second spring stiffness coefficient
% tpcSpring1CompBound     = [1e-3 1e-2];    % (m) Lower and upper bound for first spring pre-compression

Defining scaling factor for better numerical performance of the optimization algorithm

% scalingFactor = [1/tpcPosOpen1Bound(2) 1/tpcPosOpen2Bound(2) 1/tpcStiffness1CoeffBound(2), 1/tpcStiffness2CoeffBound(2) 1/tpcSpring1CompBound(2)];

Defining lower bound vector for optimization

% lowerBound = scalingFactor.*[tpcPosOpen1Bound(1) tpcPosOpen2Bound(1) tpcStiffness1CoeffBound(1) tpcStiffness2CoeffBound(1) tpcSpring1CompBound(1)];

Defining upper bound vector for optimization

% upperBound              = scalingFactor.*[tpcPosOpen1Bound(2) tpcPosOpen2Bound(2) tpcStiffness1CoeffBound(2) tpcStiffness2CoeffBound(2) tpcSpring1CompBound(2)];

Computing initialization point between lower and upper bound for the optimization algorithm for parameter estimation.

% rng default % Use the same random number generator seed
% tpcPosOpen1Rand           = tpcPosOpen1Bound(1)+(tpcPosOpen1Bound(2)-tpcPosOpen1Bound(1))*rand(1);
% tpcPosOpen2Rand           = tpcPosOpen2Bound(1)+(tpcPosOpen2Bound(2)-tpcPosOpen2Bound(1))*rand(1);
% tpcStiffness1CoeffRand    = tpcStiffness1CoeffBound(1)+(tpcStiffness1CoeffBound(2)-tpcStiffness1CoeffBound(1))*rand(1);
% tpcStiffness2CoeffRand    = tpcStiffness2CoeffBound(1)+(tpcStiffness2CoeffBound(2)-tpcStiffness2CoeffBound(1))*rand(1);
% tpcSpring1CompRand        = tpcSpring1CompBound(1)+(tpcSpring1CompBound(2)-tpcSpring1CompBound(1))*rand(1);

Defining random initialization point for the optimization

% initParam                 = scalingFactor.*[tpcPosOpen1Rand tpcPosOpen2Rand tpcStiffness1CoeffRand tpcStiffness2CoeffRand tpcSpring1CompRand];

3. Define Objective Function

The optimization algorithm minimizes the objective function to estimate the required unknown parameters. The function used here defines the least squared error between the simulation output and the expected output from the datasheet. The Error Objective subsystem computes the least squared error vector.

set_param('CylinderDesignParameterEstimation/Error Objective','LinkStatus','none')
open_system('CylinderDesignParameterEstimation/Error Objective','force')

Open objective function definition file to see the cumulative error estimation code.

The optimization algorithm assumes that the objective function has one input x, where x has as many elements as the number of unknown variables in the problem. The objective function computes the scalar value of the objective function and returns it in its single output argument y.

4. Setting up model for Fast Restart

You can save computational time by setting up the model in Fast Restart mode. Set parameter configurability to Run-Time so that optimization algorithm can be executed in Fast Restart mode. This requires modeling spring pre-compressions and preload through Ideal force Source blocks in the model.

set_param('CylinderDesignParameterEstimation/Tandem Primary Cylinder','LinkStatus','none')
open_system('CylinderDesignParameterEstimation/Tandem Primary Cylinder','force')

% load_system('CylinderDesignParameterEstimation')
% set_param('CylinderDesignParameterEstimation','FastRestart','on')

5. Parameter Estimation Using Optimization Algorithm

It is important to choose an algorithm which matches the problem complexity and type. No single optimization algorithm can perform the same for all types of problems.

Defining objective function

% ObjectiveFunction = @(x) CylinderDesignParameterEstimationObjectiveFunction(x,scalingFactor);

Defining lower bound

% lb = lowerBound ;

Defining upper bound

% ub = upperBound ;

Defining optimization initialization point

% x0 = initParam ;
% rng default      % For reproducibility in case of Heuristic method

You can use a hybrid algorithm to estimate unknown parameters. In a hybrid algorithm, multiple algorithms are re used together. For example, first a Global optimization algorithm, such as Simulated Annealing, could perform the optimization iterations, and then feed the output to a gradient based optimization algorithm such as fmincon to refine the optimization output.

% Simulated annealing
% Configuring options for Simulated annealing
% options = optimoptions('simulannealbnd','Display','iter','ObjectiveLimit',0, ...
%     'MaxIterations', 250);
% [x,fval,exitFlag,output] = simulannealbnd(ObjectiveFunction,x0,lb,ub,options)

% fmincon
% Configuring options for fmincon
% options = optimoptions('fmincon','Display','iter','Algorithm','interior-point',...
%     'StepTolerance',1e-20,'OptimalityTolerance',1e-20,...
%     'ConstraintTolerance',1e-15,'MaxIterations',200);
% x = fmincon(ObjectiveFunction,x,[],[],[],[],lb,ub,[],options)

You can also use other existing optimization algorithm functions such as Pattern Search, lsqnonlin, Genetic Algorithm, or Particle Swarm. However, you need to reconfigure optimization options for better results.

% Pattern Search
% Configuring options for Pattern Search
% options = optimoptions('patternsearch','Display','iter');
% x = patternsearch(ObjectiveFunction,x0,[],[],[],[],lb,ub,options)

% lsqnonlin
% Configuring options for lsqnonlin
% options = optimoptions('lsqnonlin','Display','iter',...
%     'Algorithm','levenberg-marquardt',...
%     'FunctionTolerance',1e-20,'StepTolerance',1e-20);
% [x,fval] = lsqnonlin(ObjectiveFunction,x0,lb,ub,options);

% Genetic Algorithm
% Configuring options for Genetic Algorithm
% options = optimoptions('ga','Display','iter','UseParallel',true,'UseVectorized', false,...
%     'FunctionTolerance',1e-9);
% x = ga(ObjectiveFunction,5,[],[],[],[],lb,ub)

% Particle Swarm
% Configuring options for Particle Swarm
% options = optimoptions('particleswarm','Display','iter','UseParallel',true,'MaxIterations', ...
%     25,'FunctionTolerance',1e-10,'MaxStallIterations',20);
% x = particleswarm(ObjectiveFunction,5,lb,ub,options)

You can use Surrogate optimization to estimate the unknown parameters. Here optimization is performed using Surrogate optimization.

% Surrogate optimization
% Configuring options for Surrogate optimization
% options = optimoptions('surrogateopt','MaxTime',600, ...
%     'MaxFunctionEvaluations',200);
% x = surrogateopt(ObjectiveFunction,lb,ub);

The estimated parameters are converted to their respective values through the defined scale.

Converting scaled parameters to actual values.

% tpcPosOpen1        = x(1)/scalingFactor(1);
% tpcPosOpen2        = x(2)/scalingFactor(2);
% tpcStiffness1Coeff = x(3)/scalingFactor(3);
% tpcStiffness2Coeff = x(4)/scalingFactor(4);
% tpcSpring1Comp     = x(5)/scalingFactor(5);
% tpcSpring2Comp     = x(3)*x(5)/(x(4)*scalingFactor(5));

% Utilizing constraint that Preload1 and Preload2
% are equal for spring-mass damper system at rest in
% steady state for estimating the last unknown parameter

Switching off Fast Restart

% set_param('CylinderDesignParameterEstimation','FastRestart','off')

Simulation

The model generates push rod force-brake pressure plots for the selected manufacturer design based on the optimization algorithm parameter estimation. The applied push rod force to the tandem primary cylinder is ramped at 25 N/sec rate in the simulation.

sim('CylinderDesignParameterEstimation',simT);

Plotting output from Simscape™ tandem primary cylinder model and comparing with function diagram

CylinderDesignParameterEstimationPlotResponse


clear

See Also

Topics