Eigenvalue problem optimized with lqsnonlin

Question

Fabio Cavagna 2023-8-3

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2004382-eigenvalue-problem-optimized-with-lqsnonlin

评论： Matt J 2024-10-22

采纳的回答： Matt J

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Torsten 2023-8-3

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

Sam Chak 2023-8-4

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 2023-8-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2004382-eigenvalue-problem-optimized-with-lqsnonlin#answer_1283047

编辑：Matt J 2023-8-3

在 MATLAB Online 中打开

% 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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Marina 2024-10-22

编辑：Matt J 2024-10-22

在 MATLAB Online 中打开

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 

Marina 2024-10-22

编辑：Marina 2024-10-22

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 2024-10-22

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

请先登录，再进行评论。

Eigenvalue problem optimized with lqsnonlin

2 个评论
显示无隐藏无

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Eigenvalue problem optimized with lqsnonlin

2 个评论 显示 无隐藏 无

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

7 个评论
显示 5更早的评论隐藏 5更早的评论