Eigenvalue problem optimized with lqsnonlin

Hello everyone,

I have the generalized eigenvalue problem (K-w^2*M)*R = 0 where

K is a 6x6 diagonal matrix whose values are unknown,

M is a 6x6 matrix whose values are known,

w^2 are the eigenvalues (typically called lambda) and they are known through f ( w = 2*pi*f )

R are the eigenvectors

I would like to obtain the diagonal values of K that minimize the difference between f_exact (my values) and f_calc evaluated through iterations.

I tried with lsqnonlin but the f_calc is far from my f_exact.

Does anyone know how to optimize in a better way?

Here's my code:

% M matrix
M = diag([1e3, 1e3, 1e3, 1e5, 1e5, 1e5]);
M(1,4) = 100; M(4,1) = 100; M(1,5) = 150; M(5,1) = 150;
M(2,6) = 10; M(6,2) = 10; M(3,6) = 1; M(6,3) = 1;
M(2,3) = 30; M(3,2) = 30; M(3,5) = 500; M(5,3) = 500;

f_exact = [30, 30, 100, 10, 10, 20]; %lambda = (2*pi*f)^2
K_guess = [1e8, 1e8, 1e8, 1e8, 1e8, 1e8]; %Initial guess
LB = [1, 1, 1, 1, 1, 1];
UB = [1e15, 1e15, 1e15, 1e15, 1e15, 1e15];

options = optimoptions('lsqnonlin','FunctionTolerance',1e-16);
K_opt = lsqnonlin(@(K) obj_function(K, M, f_exact),K_guess,LB,UB,options);

%% Check 
K_check = diag(K_opt);
[R,L] = eig(K_check,M);
f_calc = diag(sqrt(L)./(2*pi));


function diff = obj_function(K, M, f_exact)

K_mat = diag(K);
[R,L] = eig(K_mat,M);
eig_calc = diag(L);
f_calc = sqrt(eig_calc)./(2*pi)
diff = f_calc - f_exact;

end

[Torsten]

At least you should sort given and calculated eigenvalues before comparing them.

[Sam Chak]

It's good to learn optimization with lqsnonlin. However, in this eigenvalue problem, the algorithm is kind of defeating the purpose because your objective function depends eig() computation itself.

Might as well use eig() directly to find the eigenvalues. You can use eig() for checking the accuracy, but not relying on it to solve the problem.

The logical question is how can the objective function be formulated without eig()?

Answer 1

[Matt J]

% M matrix
M = diag([1e3, 1e3, 1e3, 1e5, 1e5, 1e5]);
M(1,4) = 100; M(4,1) = 100; M(1,5) = 150; M(5,1) = 150;
M(2,6) = 10; M(6,2) = 10; M(3,6) = 1; M(6,3) = 1;
M(2,3) = 30; M(3,2) = 30; M(3,5) = 500; M(5,3) = 500;

f_exact = [30, 30, 100, 10, 10, 20]; 
K_guess = ((2*pi*f_exact(:)).^2.*diag(M))'; %<---- better initial guess
LB = [1, 1, 1, 1, 1, 1];
UB = [1e15, 1e15, 1e15, 1e15, 1e15, 1e15];

options = optimoptions('lsqnonlin','FunctionTolerance',1e-16,'StepTol',1e-16, 'OptimalityTol',1e-12,'MaxFunEvals',inf);
[K_opt,~,res] = lsqnonlin(@(K) obj_function(K, M, f_exact),K_guess,LB,UB,options);

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

%% Check 
K_check = diag(K_opt);
[R,L] = eig(K_check,M);
f_calc = diag(sqrt(L)./(2*pi)).';

sort(f_calc)  %<---compare sorted

ans = 1×6

9.9997 10.0000 20.0000 30.0000 30.0000 100.0000

sort(f_exact)

ans = 1×6

10 10 20 30 30 100

function diff = obj_function(K, M, f_exact)

K_mat = diag(K);
[R,L] = eig(K_mat,M);
eig_calc = diag(L).';  %<---transpose
f_calc = sqrt(eig_calc)./(2*pi); 
diff = sort(f_calc) - sort(f_exact); %<---compare sorted

end

[Fabio Cavagna]

Thank you very much!

[Marina]

Hi,

Could you please let me know how to solve that for a non-diagonal K matrix?

where:

K = [k1 + k2, -k2, 0;

-k2, k2 + k3, -k3;

0, -k3, k3];

[Matt J]

Just modify the way K_mat is constructed,

function diff = obj_function(k, M, f_exact)

k1=k(1); k2=k(2);

K_mat = [k1 + k2, -k2, 0;
         -k2, k2 + k3, -k3;
         0, -k3, k3];

 [R,L] = eig(K_mat,M);
 eig_calc = diag(L);
 f_calc = sqrt(eig_calc)./(2*pi)
 diff = f_calc - f_exact;
 
end

[Marina]

See my code below, I am getting f values very distant from my target f values. Do you know why?

% Given data
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_exact_Hz = [f1_Hz, f2_Hz, f3_Hz]; 

% Initial Guess for Stiffness Matrix (Non-Diagonal Form)
k1_guess = ((2 * pi * f1_Hz)^2) * m1_kip_s2_ft; 
k2_guess = ((2 * pi * f2_Hz)^2) * m2_kip_s2_ft;  
k3_guess = ((2 * pi * f3_Hz)^2) * m3_kip_s2_ft; 

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

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

UB = [1e20, 1e20, 1e20];

% LB = eps*[2, -1, 0;
%      -1, 2, -1;
%       0, -1, 1];
% 
% max_val = 1e15;
% UB = max_val*[2, -1, 0;
%      -1, 2, -1;
%       0, -1, 1];
% 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
options = optimoptions('lsqnonlin','FunctionTolerance',1e-16,'StepTol',1e-16, 'OptimalityTol',1e-12,'MaxFunEvals',inf);

% Perform optimization (objective function should accept vectorized K)
% [K_opt_vec,~,res] = lsqnonlin(@(K) obj_function_reshaped(K, M, f_exact_Hz), K_guess_vec, LB_vec, UB_vec, options);
[K_opt_vec,~,res] = lsqnonlin(@(K) obj_function_reshaped_matrix(K, M, f_exact_Hz), K_guess_vec, LB_vec, UB_vec, options);

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):');
disp(K_opt);

% 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):');
disp(T_calc_s);

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

%% Objective function definition for reshaped K
function diff = obj_function_reshaped_matrix(K_vec, M, f_exact)
k1 = K_vec(1);
k2 = K_vec(2);
k3 = K_vec(3);
% Reshape the vector back to a 3x3 matrix
K_mat = [k1 + k2, -k2, 0;
    -k2, k2 + k3, -k3;
    0, -k3, k3];

% Solve eigenvalue problem
[R, L] = eig(K_mat, M);

% Extract eigenvalues and calculate frequencies
eig_calc = diag(L).';  % Transpose to row vector
f_calc = sqrt(eig_calc) ./ (2 * pi); 

% Return the difference between sorted calculated and exact frequencies
diff = (f_calc) - (f_exact);
end 

[Matt J]

No, I do not know why.

[Marina]

Do you have any suggestion on how I can get my K stiffness matrix? with a better exact vs calculated values. Thank you!

[Matt J]

I'm afraid not. I have no reason to believe it should work.