Using Optimization to Tune a Suspension System

by Paul Goossens and Jason Schattman

WMI Marketing Dept.

Introduction

Suppose that we wish to design an automobile suspension system that exhibits some specified behavior in response to a bump in the road. As engineers, we have control over the spring constant k and the damper constant c . Given the mass m of the car and the expected amplitude of a typical bump, our task is to choose k and c so that the suspension system's response is as close as possible to the desired response.

Our approach to this engineering problem is to measure the discrepancy between the desired and actual motion of the system as a function of k and c. After deriving the actual response by solving the system's differential equation, we use an optimization algorithm to find the values of k and c that minimize this discrepancy. In this application, we will use Newton's Method for unconstrained nonlinear programming as the optimization algorithm.

Below is a diagram of the suspension system. x(t) denotes the amplitude of the automobile as a function of time.

An Example Desired Response

For illustration purposes, suppose our target response is a decaying exponential function. That is, when the automobile hits a bump of unit amplitude, we want it to spring up and down with an exponentially shrinking amplitude. (In practice, the desired response need not be an algebraic function. It could simply be a discrete set of ordered pairs of amplitudes and times.)

> restart;

> Target:=t->exp(-.2*t)*cos(t);

> plot(Target(t), t = 0..20, title=`Target Response`, thickness=3);

Deriving the Actual Response

To derive the actual response of the system to a unit bump, we must solve the differential equation governing the system's behavior with the initial condition x(0) = 1 . Below is a more general equation describing the system under a periodic disturbance with frequency and amplitude .

> System_Equation := diff(x(t),`$`(t,2))+2*zeta*omega[n]*diff(x(t),t)+omega[n]^2*x(t) = F[0]/m*sin(omega*t);

where:

and

To solve this differential equation, we've developed a package specifically for analyzing vibrating systems.

> libname:="c:/mylib/vibsys",libname:
with(plots):
with(VibratingSystem);

Warning, the name changecoords has been redefined

Our initial conditions are unit amplitude and zero velocity:

> Xinit:=1; Xdotinit:=0;

We now solve the system and define the response as a function of k , c and t, with m=100 and F0 = 0.

> System_Response:=unapply(Solve_System(Xinit,Xdotinit),(k,c,m,omega,F0,t));

> Actual:=(k, c, t)->System_Response(k, c, 100, 1, 0, t);

Let's now overlay plots of the actual response with k=10 and c=10 , and the target response. Notice that these choices of k and c do not provide a good match to the target. One approach to improving the match is to try different combinations of k and c until we're satisfied. However, we use a more systematic approach that utilizes optimization techniques.

> plot([ Target(t), Actual(10, 10, t) ], t=0..20, thickness=3,
title=`Target response vs. Actual response for k = 10 and c = 10`);

Measuring the Error between the Desired and Actual Responses

Squared differences are a common measure of error between two functions. Ideally, we would measure error as the integral of the squared difference between the desired and actual response. However, evaluating this integral in this application is extremely difficult, due to the high complexity of the actual response function. As an approximation, we sample the functions at key times and sum up the squared differences between the functions at these points. In this case as an illustration, we sample at the peaks, zeros and valleys of the target function over two periods.

> period:=2*Pi; num_periods:=2;

> time_sample:=seq( i * period / 4, i=1..4 * num_periods );

> Error:=add( (Actual(k, c, time_sample[i])
- Target(time_sample[i]))^2,
i = 1..4 * num_periods);

Minimizing the Error with Nonlinear Programming

We now minimize the function Error(k,c) using Newton's Method for unconstrained nonlinear programming. This algorithm is contained in our package NonlinearProgramming .

> Error:=subs( {k=x[1],c=x[2]} , Error):

> infolevel['UnconstrainedNewton']:=2;

The algorithm makes use of the LinearAlgebra and linalg packages, so we load them now.

> with(LinearAlgebra): with(linalg):

Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value

Warning, the protected names norm and trace have been redefined and unprotected

> libname:="C:/mylib/nlp", libname;

> with( NonlinearProgramming );

The parameters to Newton's Method are the objective function ( Error(k,c) ), the number of decision variables (two), the starting point for the search (<10, 10>), the convexity of the objective function (nonconvex, in this case), and the data type of the variables (complex floats). We use complex floats instead of real because Error(k,c) can take on complex values. In this case, the algorithm takes the real part of all function evaluations.

(Because Error(k,c) is very complicated, the algorithm takes about 10 minutes to run on a PC. The bottle neck is computing the Hessian matrix of second derivatives.)

> UnconstrainedNewton( Error, 2, <10,10>, 'nonconvex', 'complex[8]' );

Testing the Result

Newton's Method found a local minimum at k = 113.101, c = 40.013 . We now overlay the desired and actual plots with these values and see that the match is quite good.

> plot([Target(t), Actual(113, 40, t) ], t=0..20, thickness=3);

>