Есть ли ограничение на количество степеней свободы с помощью алгоритма lm_feasible?Если так, то каков предел? - PullRequest
Новая
       42

Есть ли ограничение на количество степеней свободы с помощью алгоритма lm_feasible?Если так, то каков предел?

8 голосов
/ 11 июля 2019

Я занимаюсь разработкой программного обеспечения для конечных элементов, которое минимизирует энергию механической конструкции.Используя octave и его пакет optim, я столкнулся со странной проблемой: алгоритм lm_feasible вообще не рассчитывается, когда я использую более 300 степеней свободы (DoF).Другой алгоритм (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;')

1 Ответ

0 голосов
/ 15 июля 2019

Метод конечных элементов обычно формулируется как оптимальный критерий для задачи минимизации, что эквивалентно принципу виртуальной работы (см. Книги, подобные Хьюзу из Бате). Виртуальная работа представляет собой набор линейных (или нелинейных) уравнений, которые могут быть решены более эффективно (с помощью fsolve).

Если по какому-либо мотиву вы должны решить задачу как задачу оптимизации, то, если вы рассматриваете линейную упругость, ваша энергия деформации является квадратичной, поэтому вы можете использовать функцию qp Octave.

Также полезно использовать разреженные матрицы.

...