Please help for objectiveFunction command.

% Optimization of Spring Constants and Damping Coefficients for Minimum Acceleration
% Given upper and lower boundaries for spring constants and damping coefficients
lower_bound = [1000, 1000, 50, 50];
upper_bound = [5000, 5000, 200, 200];
% Objective function for optimization (minimize vertical acceleration)
objective = @(x) objectiveFunction(x, k1, k2, c1, c2);
% Initial guess for optimization variables [k1, k2, c1, c2]
x0 = [3000, 2500, 150, 100];
% Optimization options
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'interior-point');
% Perform optimization
[x_opt, fval] = fmincon(objective, x0, [], [], [], [], lower_bound, upper_bound, [], options);

Unrecognized function or variable 'k1'.

Error in solution>@(x)objectiveFunction(x,k1,k2,c1,c2) (line 6)
objective = @(x) objectiveFunction(x, k1, k2, c1, c2);

Error in fmincon (line 563)
initVals.f = feval(funfcn{3},X,varargin{:});

Caused by:
Failure in initial objective function evaluation. FMINCON cannot continue.

% Display optimized results
disp('Optimized Spring Constants and Damping Coefficients:');
disp(['k1: ', num2str(x_opt(1)), ' N/m']);
disp(['k2: ', num2str(x_opt(2)), ' N/m']);
disp(['c1: ', num2str(x_opt(3)), ' Ns/m']);
disp(['c2: ', num2str(x_opt(4)), ' Ns/m']);
disp(['Optimized Objective Value: ', num2str(fval)]);
% Simulate the system with the optimized parameters and obtain acceleration before optimization
[x_sim_before, v_sim_before, a_sim_before] = simulateQuarterCarModel(m1, m2, k1, k2, c1, c2);
% Simulate the system with the optimized parameters and obtain acceleration after optimization
[x_sim_after, v_sim_after, a_sim_after] = simulateQuarterCarModel(m1, m2, x_opt(1), x_opt(2), x_opt(3), x_opt(4));
% Plot results
figure;
% Plot acceleration before optimization
subplot(3, 1, 1);
plot(t, a_sim_before, 'b', 'LineWidth', 1.5);
title('Acceleration of Sprung Mass (Before Optimization)');
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
grid on;
% Plot acceleration after optimization
subplot(3, 1, 2);
plot(t, a_sim_after, 'r', 'LineWidth', 1.5);
title('Acceleration of Sprung Mass (After Optimization)');
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
grid on;
% Plot road profile
subplot(3, 1, 3);
plot(t, road_profile(t), 'g', 'LineWidth', 1.5);
title('Road Profile');
xlabel('Time (s)');
ylabel('Road Profile');
grid on;
function J = objectiveFunction(x, m1, m2, ~, ~)
% Simulate the system and compute the objective function
[~, ~, a_sim] = simulateQuarterCarModel(m1, m2, x(1), x(2), x(3), x(4));

% Objective: Minimize the vertical acceleration of the sprung mass
J = max(abs(a_sim));
end
function [x_sim, v_sim, a_sim] = simulateQuarterCarModel(m1, ~, k1, k2, c1, c2)
% Simulation parameters
dt = 0.01;  % time step (s)
t_end = 5;  % simulation time (s)
t = 0:dt:t_end;
% Preallocate arrays
x_sim = zeros(size(t));
v_sim = zeros(size(t));
a_sim = zeros(size(t));
% Initial conditions
x_sim(1) = 0.1;   % initial displacement of the sprung mass (m)
v_sim(1) = 10;   % initial velocity of the sprung mass (m/s)
% Simulation loop
for i = 1:length(t)-1
    % State-space equations
    A = [0, 1; -((k1 + k2)/m1), -(c1 + c2)/m1];
    B = [0; k2/m1];
    u = road_profile(t(i)); % road profile as an input
    
    % Update state using Euler integration
    x_dot = A * [x_sim(i); v_sim(i)] + B * u;
    x_sim(i+1) = x_sim(i) + v_sim(i) * dt;
    v_sim(i+1) = v_sim(i) + x_dot(2) * dt;
    a_sim(i+1) = x_dot(2);
end
end
function y = road_profile(t)
% Define a simple road profile function (you can replace this with your own)
y = 0.1 * sin(2 * pi * 0.5 * t);
end

I write this and end up to the following:

Error in newopt>@(x)objectiveFunction(x,k1,k2,c1,c2) (line 8)

objective = @(x) objectiveFunction(x, k1, k2, c1, c2);

Error in fmincon (line 573)

initVals.f = feval(funfcn{3},X,varargin{:});

Error in newopt (line 17)

[x_opt, fval] = fmincon(objective, x0, [], [], [], [], lower_bound, upper_bound, [], options);

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

[Dyuman Joshi]

The error is clear, k1 is undefined. So are k2, c1, and c2.

lower_bound = [1000, 1000, 50, 50];
upper_bound = [5000, 5000, 200, 200];
% Objective function for optimization (minimize vertical acceleration)

%                                     vv  vv  vv  vv
objective = @(x) objectiveFunction(x, k1, k2, c1, c2);

It's not clear to me what is being optimized.

[SOUVIK]

I want to minimize "Acceleration of Sprung Mass()" for given boundary values of k1, k2, c1, c2 and for given masses (m1 & m2). Please help me with the code. Thanks in advance.

Answer 1

[Sam Chak]

Hi @SOUVIK

The input argument to call fmincon() is incomplete. I prefer the ode45 solver over the Euler integration method. However, you can still use your own user-defined models to simulate before and after optimization.

%% Optimization of Spring Constants and Damping Coefficients for Minimum Acceleration

%% Settings
objfun  = @costfun;                 % tag the objective function
x0      = [3000, 2500, 150, 100];   % Initial guess for [k1, k2, c1, c2]
nvars   = 4;                        % number of optimization variables
A       = - eye(nvars);
b       = zeros(nvars, 1);          % A*x <= b
Aeq     = [];                       % nil
beq     = [];                       % nil
lb      = [1000, 1000,  50,  50];   % lower bound
ub      = [5000, 5000, 200, 200];   % upper bound
nonlcon = [];                       % nil
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'interior-point');

%% Call fmincon to solve the optimization problem
[x_opt, fval] = fmincon(objfun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options);

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       5    1.792113e+01    0.000e+00    3.267e-03
    1      10    1.792111e+01    0.000e+00    3.267e-03    5.262e-03
    2      15    1.792097e+01    0.000e+00    3.267e-03    2.638e-02
    3      20    1.792027e+01    0.000e+00    3.267e-03    1.319e-01
    4      25    1.791680e+01    0.000e+00    3.268e-03    6.596e-01
    5      30    1.789943e+01    0.000e+00    3.270e-03    3.297e+00
    6      35    1.781259e+01    0.000e+00    3.282e-03    1.647e+01
    7      40    1.749494e+01    0.000e+00    3.323e-03    5.999e+01
    8      46    1.749414e+01    0.000e+00    3.323e-03    1.501e-01
    9      51    1.749450e+01    0.000e+00    3.323e-03    6.671e-02
   10      56    1.749407e+01    0.000e+00    3.323e-03    8.818e-02
   11      61    1.749293e+01    0.000e+00    3.323e-03    2.733e-01
   12      66    1.748728e+01    0.000e+00    3.324e-03    1.366e+00
   13      71    1.745901e+01    0.000e+00    3.326e-03    6.826e+00
   14      76    1.731811e+01    0.000e+00    3.339e-03    3.401e+01
   15      81    1.723596e+01    0.000e+00    3.348e-03    1.982e+01
   16      89    1.723703e+01    0.000e+00    3.348e-03    2.519e-01
   17      95    1.723638e+01    0.000e+00    2.193e-03    1.768e-01
   18     100    1.723708e+01    0.000e+00    1.844e-03    1.323e-01
   19     105    1.723684e+01    0.000e+00    1.844e-03    1.468e-01
   20     110    1.723557e+01    0.000e+00    1.844e-03    5.221e-01
   21     115    1.722922e+01    0.000e+00    1.846e-03    2.611e+00
   22     120    1.719740e+01    0.000e+00    1.850e-03    1.306e+01
   23     125    1.703701e+01    0.000e+00    1.869e-03    6.543e+01
   24     130    1.620010e+01    0.000e+00    1.977e-03    3.306e+02
   25     135    1.137025e+01    0.000e+00    2.480e-03    1.635e+03
   26     140    1.134177e+01    0.000e+00    2.486e-03    8.173e+00
   27     145    1.134346e+01    0.000e+00    2.485e-03    8.318e-01
   28     150    1.134219e+01    0.000e+00    2.463e-03    5.505e-01
   29     155    1.134207e+01    0.000e+00    2.463e-03    4.958e-02
   30     160    1.134147e+01    0.000e+00    2.463e-03    2.452e-01

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31     165    1.133847e+01    0.000e+00    2.462e-03    1.226e+00
   32     170    1.132346e+01    0.000e+00    2.459e-03    6.132e+00
   33     175    1.124802e+01    0.000e+00    2.463e-03    3.070e+01
   34     180    1.086103e+01    0.000e+00    2.546e-03    1.544e+02
   35     185    1.009254e+01    0.000e+00    2.722e-03    2.918e+02
   36     191    1.009055e+01    0.000e+00    1.840e-03    7.295e-01
   37     197    1.008970e+01    0.000e+00    9.479e-04    3.700e-01
   38     202    1.008687e+01    0.000e+00    3.395e-04    5.256e-01
   39     207    1.008607e+01    0.000e+00    2.850e-06    1.634e-01
   40     212    1.008606e+01    0.000e+00    6.106e-08    1.474e-03

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

%% Display optimized results
disp('Optimized Spring Constants and Damping Coefficients:');

Optimized Spring Constants and Damping Coefficients:

disp(['k1: ', num2str(x_opt(1)), ' N/m']);

k1: 1000 N/m

disp(['k2: ', num2str(x_opt(2)), ' N/m']);

k2: 1000 N/m

disp(['c1: ', num2str(x_opt(3)), ' Ns/m']);

c1: 200 Ns/m

disp(['c2: ', num2str(x_opt(4)), ' Ns/m']);

c2: 200 Ns/m

disp(['Optimized Objective Value: ', num2str(fval)]);

Optimized Objective Value: 10.0861

%% Objective Function
function J = costfun(param) % cost function must contain the parameters to be optimized
    % first things first is to define the parameters
    k1     = param(1);  % spring constant of shock absorber 1
    k2     = param(2);  % spring constant of shock absorber 2
    c1     = param(3);  % damping coeff   of shock absorber 1
    c2     = param(4);  % damping coeff   of shock absorber 2

    tspan  = [0 10];
    x0     = [0.1; 10];
    [t, x] = ode45(@(t, x) odefcn(t, x, k1, k2, c1, c2), tspan, x0);

    for j  = 1:length(t)
        dxdt(j,:) = odefcn(t(j), x(j,:).', k1, k2, c1, c2);
    end
    
    % double check if this maximum absolute acceleration is to minimize it
    J      = max(abs(dxdt(:,2)));
end

%% Quarter Car Model
function dxdt = odefcn(t, x, k1, k2, c1, c2)    % mass m2 is missing!!
    m1   = 1500;    % [kg] average car weight
    A    = [ 0,             1; 
          -(k1 + k2)/m1, -(c1 + c2)/m1];
    B    = [ 0; 
            k2/m1];
    u    = 0.1*sin(2*pi*0.5*t);    % road profile as an input
    dxdt = A*x + B*u;
end

[SOUVIK]

Thank you, but "Optimized Objective Value: 10.0861" -- whose value is it?

[Sam Chak]

Hi @SOUVIK

That value comes from the Objective Function:

J = max(abs(dxdt(:,2)));

which should be equivalent to yours:

J = max(abs(a_sim));

Basically, fmincon found that , , , return the minimum maximal value of the absolute acceleration, .

However, since the optimal values of the parameters lie on the boundaries, you may need to either shift the lower and upper bounds, revisit/redesign the objective function, or do both.

[SOUVIK]

Thank you!

Answer 2

[Torsten]

You don't supply k1, k2, c1 and c2 to the objective function.

If these are the parameters to be optimized, you don't need to pass them since they are stored in the vector x.

[SOUVIK]

I omitted it, still getting a problem. Please tell me all the possible solutions.

% Optimization of Spring Constants and Damping Coefficients for Minimum Acceleration

% Given upper and lower boundaries for spring constants and damping coefficients
lower_bound = [1000, 1000, 50, 50];
upper_bound = [5000, 5000, 200, 200];

% Objective function for optimization (minimize vertical acceleration)
objective = @(x) objectiveFunction(x);

% Initial guess for optimization variables [k1, k2, c1, c2]
x0 = [3000, 2500, 150, 100];

% Optimization options
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'interior-point');

% Perform optimization
[x_opt, fval] = fmincon(objective, x0, lower_bound, upper_bound, [], options);

Error using fmincon
FMINCON requires the following inputs to be of data type double: 'Beq'.

% Display optimized results
disp('Optimized Spring Constants and Damping Coefficients:');
disp(['k1: ', num2str(x_opt(1)), ' N/m']);
disp(['k2: ', num2str(x_opt(2)), ' N/m']);
disp(['c1: ', num2str(x_opt(3)), ' Ns/m']);
disp(['c2: ', num2str(x_opt(4)), ' Ns/m']);
disp(['Optimized Objective Value: ', num2str(fval)]);

% Simulate the system with the optimized parameters and obtain acceleration before optimization
[x_sim_before, v_sim_before, a_sim_before] = simulateQuarterCarModel(m1, m2, k1, k2, c1, c2);

% Simulate the system with the optimized parameters and obtain acceleration after optimization
[x_sim_after, v_sim_after, a_sim_after] = simulateQuarterCarModel(m1, m2, x_opt(1), x_opt(2), x_opt(3), x_opt(4));

% Plot results
figure;

% Plot acceleration before optimization
subplot(3, 1, 1);
plot(t, a_sim_before, 'b', 'LineWidth', 1.5);
title('Acceleration of Sprung Mass (Before Optimization)');
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
grid on;

% Plot acceleration after optimization
subplot(3, 1, 2);
plot(t, a_sim_after, 'r', 'LineWidth', 1.5);
title('Acceleration of Sprung Mass (After Optimization)');
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
grid on;

% Plot road profile
subplot(3, 1, 3);
plot(t, road_profile(t), 'g', 'LineWidth', 1.5);
title('Road Profile');
xlabel('Time (s)');
ylabel('Road Profile');
grid on;

function J = objectiveFunction(x, m1, m2, ~, ~)
    % Simulate the system and compute the objective function
    [~, ~, a_sim] = simulateQuarterCarModel(m1, m2, x(1), x(2), x(3), x(4));
    
    % Objective: Minimize the vertical acceleration of the sprung mass
    J = max(abs(a_sim));
end

function [x_sim, v_sim, a_sim] = simulateQuarterCarModel(m1, ~, k1, k2, c1, c2)
    % Simulation parameters
    dt = 0.01;  % time step (s)
    t_end = 5;  % simulation time (s)
    t = 0:dt:t_end;

    % Preallocate arrays
    x_sim = zeros(size(t));
    v_sim = zeros(size(t));
    a_sim = zeros(size(t));

    % Initial conditions
    x_sim(1) = 0.1;   % initial displacement of the sprung mass (m)
    v_sim(1) = 10;   % initial velocity of the sprung mass (m/s)

    % Simulation loop
    for i = 1:length(t)-1
        % State-space equations
        A = [0, 1; -((k1 + k2)/m1), -(c1 + c2)/m1];
        B = [0; k2/m1];
        u = road_profile(t(i)); % road profile as an input
        
        % Update state using Euler integration
        x_dot = A * [x_sim(i); v_sim(i)] + B * u;
        x_sim(i+1) = x_sim(i) + v_sim(i) * dt;
        v_sim(i+1) = v_sim(i) + x_dot(2) * dt;
        a_sim(i+1) = x_dot(2);
    end
end

function y = road_profile(t)
    % Define a simple road profile function (you can replace this with your own)
    y = 0.1 * sin(2 * pi * 0.5 * t);
end

[Steven Lord]

You can't simply leave out inputs you don't want to specify or pass the inputs in an arbitrary order when you call fmincon. The way you're calling it:

[x_opt, fval] = fmincon(objective, x0, lower_bound, upper_bound, [], options);

attempts to pass the variable lower_bound in the place fmincon expects the left-hand side of the linear inequality constraints, upper_bound in the place fmincon expects the right-hand side of the linear inequality constraints, and options in the place where fmincon expects the right-hand side of the linear equality constraints. fmincon doesn't "look at the names" or anything like that and try to guess what you meant, it does what you told it to do.

If you don't want to specify one or more of those inputs, you must include a [] in place of the input you don't want to specify.

[SOUVIK]

Okay, but please give a perfect solution from your point of view. Please give me the rectified codes (Only then I can understand my fault).

I want to minimize "Acceleration of Sprung Mass()" for given boundary values of k1, k2, c1, c2 and for given masses (m1 & m2). Please help me with the code. Thanks in advance.