时间戳

首页 > 解决方案 > lm_feasible 算法的自由度数量是否有限制？如果是这样，限制是多少？

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

问题描述

我正在开发一种将机械结构的能量降至最低的有限元软件。使用 octave 及其 optim 包，我遇到了一个奇怪的问题：当我使用超过 300 个自由度 (DoF) 时，lm_feasible 算法根本不计算。另一种算法（sqp）执行计算，但当我复杂化结构并且超出我的测试用例时效果不佳。

lm_feasible 算法的自由度数量有限制吗？

如果是这样，最大可能有多少自由度？

概述和大致了解代码的工作原理：

[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;')

标签： octavenonlinear-optimizationfinite-element-analysis

解决方案

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

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

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

推荐阅读