octave - lm_feasible 算法的自由度数有限制吗？如果是这样，限制是多少？

Question

我正在开发一个有限元软件，可以最大限度地减少机械结构的能量。使用 Octave 及其优化包，我遇到了一个奇怪的问题:当我使用超过 300 个自由度 (DoF) 时，lm_feasible 算法根本不计算。另一种算法 (sqp) 执行计算，但当我使结构复杂化并且不在我的测试用例中时效果不佳。

lm_feasible 算法的 DoF 数量有限制吗？

如果是这样，最大可能有多少 DoF？

给出代码如何工作的概述和一般概念:

[x,y] = geometryGenerator()

U = zeros(lenght(x)*2,1);
U(1:2:end-1) = x;
U(2:2:end) = y;

%Non geometric argument are not optimised, and fixed during calculation
fct =@(U)complexFunctionOfEnergyIWrap(U(1:2:end-1),U(2:2:end), variousMaterialPropertiesAndOtherArgs)

para = optimset("f_equc_idx",contEq,"lb",lb,"ub",ub,"objf_grad",dEne,"objf_hessian",d2Ene,"MaxIter",1000);
[U,eneFinale,cvg,outp] = nonlin_min(fct,U,para)

完整示例:

clear

pkg load optim

function [x,y] = geometryGenerator(r,elts = 100)
  teta  = linspace(0,pi,elts = 100);
  x = r * cos(teta);
  y = r * sin(teta);
endfunction

function ene  = complexFunctionOfEnergyIWrap (x,y,E,P, X,Y)
  ene = 0;
  for i = 1:length(x)-1
    ene += E*(x(i)/X(i))^4+ E*(y(i)/Y(i))^4- P *(x(i)^2+(x(i+1)^2)-x(i)*x(i+1))*abs(y(i)-y(i+1));
  endfor
endfunction

[x,y] = geometryGenerator(5,100)

%Little distance from axis to avoid division by zero
x +=1e-6;
y +=1e-6;
%Saving initial geometry
X = x;
Y = y;

%Vectorisation of the function
%% Initial vector
U = zeros(length(x)*2,1);
U(1:2:end-1) = linspace(min(x),max(x),length(x));
U(2:2:end) = linspace(min(y),max(y),length(y));

%%Constraints
Aeq = zeros(3,length(U));
%%% Blocked bottom
    Aeq(1,1) = 1;
    Aeq(2,2) = 1;
%%% Sliding top    
    Aeq(3,end-1) = 1;
%%%Initial condition
    beq = zeros(3,1);
    beq(1) = U(1);
    beq(2) = U(2);
    beq(3) = U(end-1);

    contEq = @(U) Aeq * U - beq;

%Parameter
Mat = 0.2e9;
pressure = 50;

%% Vectorized function. Non geometric argument are not optimised, and fixed during calculation
fct =@(U)complexFunctionOfEnergyIWrap(U(1:2:end-1),U(2:2:end), Mat, pressure, X, Y)

para = optimset("Algorithm","lm_feasible","f_equc_idx",contEq,"MaxIter",1000);
[U,eneFinale,cvg,outp] = nonlin_min(fct,U,para)

xFinal = U(1:2:end-1);
yFinal = U(2:2:end);

plot(x,y,';Initial geo;',xFinal,yFinal,'--x;Final geo;')

Answer 1

有限元方法通常被表述为最小化问题的最佳标准，这相当于虚拟工作原理(参见 Hughes of Bathe 等书籍)。 Virtual Work 表示一组线性(或非线性)方程，可以更有效地求解(使用 fsolve)。

如果出于某种原因，您必须将问题作为优化问题来解决，那么，如果您正在考虑线性弹性，则应变能是二次的，因此您可以使用 qp Octave 函数。

使用稀疏矩阵也可能会有所帮助。