Optimization Eigenvalue problem Stiffness matrix

10 Ansichten (letzte 30 Tage)
Marina
Marina am 22 Okt. 2024
Beantwortet: Sam Chak am 22 Okt. 2024
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
  • 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

Antworten (2)

Torsten
Torsten am 22 Okt. 2024
Bearbeitet: Torsten am 22 Okt. 2024
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

Sam Chak
Sam Chak am 22 Okt. 2024
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

Kategorien

Mehr zu Linear Algebra finden Sie in Help Center und File Exchange

Produkte


Version

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by