Optimization Eigenvalue problem Stiffness matrix

I have a eigenvalue problem of the form (K - w²*M)*R = 0, where:

Please see a reference question. They do not use non-diagonal stiffness matrix

Eigenvalue problem optimized with lqsnonlin - MATLAB Answers - MATLAB Central

• K is a 3x3 diagonal matrix with unknown values,
• M is a known 3x3 matrix,
• w² represents the eigenvalues (typically denoted as lambda) and is related to the known frequencies f (where w = 2pif),
• R are the corresponding eigenvectors.

My goal is to determine the diagonal values of K that minimize the difference between my exact frequency values (f_exact) and the calculated frequency values (f_calc), which are computed iteratively.

I’ve attempted to use lsqnonlin for optimization, but the calculated frequencies (f_calc) are significantly off from my target values (f_exact). Does anyone have suggestions on how to improve the optimization approach?

Here is the code I used:

% Marina Krauze Thomaz
% REF: https://www.mathworks.com/matlabcentral/answers/2004382-eigenvalue-problem-optimized-with-lqsnonlin
% Given data

m1_kip_s2_ft = 1035;
m2_kip_s2_ft = 628;
m3_kip_s2_ft = 591;

% m1 = 1035;
% m2 = 628;
% m3 = 591;

M = diag([m1_kip_s2_ft, m2_kip_s2_ft, m3_kip_s2_ft]);  % Mass matrix

T_wanted1_s = 5.2;
T_wanted2_s = 4.0;
T_wanted3_s = 1.17;
T_wanted_s  = [T_wanted1_s; T_wanted2_s; T_wanted3_s];  % Target periods

f1_Hz = 1/T_wanted1_s;
f2_Hz = 1/T_wanted2_s;
f3_Hz = 1/T_wanted3_s;

f_exact_Hz = [f1_Hz, f2_Hz, f3_Hz]; 

% Improved Initial Guess for Stiffness Matrix (Non-Diagonal Form)
% Use a heuristic guess for k1, k2, and k3 based on mass and frequency
k1_guess = (1);  % Guess 
k2_guess = (1);  % Guess 
k3_guess = (1);  % Guess 

% Combine into initial guess vector [k1, k2, k3]
K_guess_kip_ft = [k1_guess + k2_guess, -k2_guess,             0;
                 -k2_guess,             k2_guess + k3_guess, -k3_guess;
                  0,                   -k3_guess,             k3_guess];

% Bounds for K
LB = [1, 1, 1];
UB = [1e16, 1e16, 1e16];  

% Define bounds and flatten them for optimization (vectorize)
LB_vec = LB(:);
UB_vec = UB(:);
K_guess_vec = K_guess_kip_ft(:);  % Vectorize the initial guess matrix

% Optimization options with relaxed tolerances
options = optimoptions('lsqnonlin', 'FunctionTolerance', 1e-8, 'StepTolerance', 1e-8, ...
    'OptimalityTolerance', 1e-6, 'MaxFunEvals', inf, 'MaxIterations', 5000);

% Perform optimization
[K_opt_vec,~,res] = lsqnonlin(@(K) obj_function_reshaped_matrix(K, M, f_exact_Hz), K_guess_vec, LB_vec, UB_vec, options);

Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.

Warning: Length of upper bounds is < length(x); filling in missing upper bounds with +Inf.

Solver stopped prematurely.

lsqnonlin stopped because it exceeded the iteration limit,
options.MaxIterations = 5.000000e+03.

%% Reshape the optimized vector back into a matrix
K_opt =     [K_opt_vec(1) + K_opt_vec(2), -K_opt_vec(2),                 0;
            -K_opt_vec(2),                 K_opt_vec(2) + K_opt_vec(3), -K_opt_vec(3);
             0,                           -K_opt_vec(3),                 K_opt_vec(3)];

% Check the optimized stiffness matrix
disp('Optimized Stiffness Matrix (K):');

Optimized Stiffness Matrix (K):

disp(K_opt);

   1.0e+16 *

    0.5000   -0.5000         0
   -0.5000    1.0000   -0.5000
         0   -0.5000    0.5000

% Solve eigenvalue problem with optimized K
[R, L] = eig(K_opt, M);
f_calc_Hz = (diag(sqrt(L) ./ (2 * pi)).');  % Do not sort to preserve order
T_calc_s = 1 ./f_calc_Hz;                   % Periods as a row vector

%% Error
% Calculate the squared error between target and calculated periods
error_squared = (T_wanted_s - T_calc_s.').^2; %Transpose T_calc_s

% Display calculated periods and error
disp('Calculated periods (T_calc_s):');

Calculated periods (T_calc_s):

disp(T_calc_s);

  127.2559    0.0000    0.0000

disp('Squared error between target and calculated periods:');

Squared error between target and calculated periods:

disp(error_squared);

   1.0e+04 *

    1.4898
    0.0016
    0.0001

% Check the optimized stiffness matrix
disp('Optimized Stiffness Matrix (K):');

Optimized Stiffness Matrix (K):

disp(K_opt);

   1.0e+16 *

    0.5000   -0.5000         0
   -0.5000    1.0000   -0.5000
         0   -0.5000    0.5000

%% Objective function definition
function diff = obj_function_reshaped_matrix(K_vec, M, f_exact)
    % Stiffness values
    k1 = K_vec(1);
    k2 = K_vec(2);
    k3 = K_vec(3);
    
    % Construct stiffness matrix
    K_mat = [k1 + k2, -k2,       0;
            -k2,       k2 + k3, -k3;
             0,       -k3,       k3];
    
    % Solve eigenvalue problem
    [~, L] = eig(K_mat, M);
    
    % Calculate natural frequencies from eigenvalues
    omega_calc = sqrt(diag(L));             % Rad/s
    f_calc_Hz  = omega_calc / (2 * pi);     % Convert to Hz
    
    % Return squared frequency differences (to minimize)
    diff = (f_calc_Hz - f_exact).^2;        % Use squared difference
end

Answer 1

[Torsten]

Should be

% Perform optimization
[K_opt_vec,~,res] = lsqnonlin(@(K) obj_function_reshaped_matrix(K, M, f_exact_Hz), [k1_guess,k2_guess,k3_guess], LB_vec, UB_vec, options);

instead of

% Perform optimization
[K_opt_vec,~,res] = lsqnonlin(@(K) obj_function_reshaped_matrix(K, M, f_exact_Hz), K_guess_vec, LB_vec, UB_vec, options);

and

% Return squared frequency differences (to minimize)
diff = f_calc_Hz - f_exact;        % Use squared difference

instead of

% Return squared frequency differences (to minimize)
diff = (f_calc_Hz - f_exact).^2;        % Use squared difference

(squaring is done internally by "lsqnonlin").

And I think it's necessary to compare sort(f_calc_Hz) with sort( f_exact), thus

diff = sort(f_calc_Hz) - sort(f_exact);        % Use squared difference

Answer 2

[Sam Chak]

Hi @Marina

The problem has two parts. I have solved only the first part for you, which is to determine the characteristic polynomial that yields the desired natural frequencies. You need to solve the second part of the problem, which involves finding the stiffness matrix that produces the target characteristic polynomial. However, my preliminary analysis indicates that there is no real-valued solution. My analysis could be inaccurate; please double-check this.

format long g

%% Desired natural frequencies
m1_kip_s2_ft= 1035;
m2_kip_s2_ft= 628;
m3_kip_s2_ft= 591;
M           = diag([m1_kip_s2_ft, m2_kip_s2_ft, m3_kip_s2_ft]);  % Mass matrix
T_wanted1_s = 5.2;
T_wanted2_s = 4.0;
T_wanted3_s = 1.17;
T_wanted_s  = [T_wanted1_s; T_wanted2_s; T_wanted3_s];  % Target periods
f1_Hz       = 1/T_wanted1_s;
f2_Hz       = 1/T_wanted2_s;
f3_Hz       = 1/T_wanted3_s;
f           = [f1_Hz; f2_Hz; f3_Hz] % <-- Desired natural frequencies

f = 3×1

0.192307692307692 0.25 0.854700854700855

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Desired Characteristic Polynomial, CP
omega       = 2*pi*f;
myRoots     = omega.^2;
% r1        = 1.46000065104872
% r2        = 2.46740110027234
% r3        = 28.8395190330612
CP          = poly(myRoots)     % <-- target Characteristic Polynomial

CP = 1×4

1 -32.7669207843822 116.866784770497 -103.891691378266

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Testing results (according to OP's code)
sol         = roots(CP);
omega_calc  = sqrt(sol);             % rad/s
f_calc_Hz   = sort(omega_calc/(2*pi))

f_calc_Hz = 3×1

0.192307692307692 0.25 0.854700854700855

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Part 2: Need to solve for {k1, k2, k3} to satisfy CP
syms k1 k2 k3 real positive

K   = [k1 + k2, -k2,       0;
      -k2,       k2 + k3, -k3;
       0,       -k3,       k3];
M   = diag([m1_kip_s2_ft, m2_kip_s2_ft, m3_kip_s2_ft]);  % Mass matrix
MK  = M\K;
cp  = charpoly(MK)

%% One of six complex-valued solution sets
k1  = 2031.99 - 2061.74i;
k2  = 1089.46 + 1192.27i;
k3  = 8530.08 - 322.27i;
K_mat   = [k1 + k2, -k2,       0;
          -k2,       k2 + k3, -k3;
           0,       -k3,       k3];
[~, L]  = eig(K_mat, M);

%% Calculate natural frequencies from eigenvalues
omega_calc = sqrt(diag(L));     % rad/s
f_calc_Hz  = omega_calc/(2*pi)  % Convert to Hz

f_calc_Hz = 
          0.854701005410721 -  6.07592924005276e-08i
          0.250000256584473 -  8.05443027787641e-07i
          0.192307431901452 +  4.33118420843562e-07i