Is there a limit in the number of degrees of freedom with the lm_feasible algorithm? If so, what is the limit?

Question

I am developing a finite element software that minimizes the energy of a mechanical structure. Using octave and its optim package, I run into a strange issue: The lm_feasible algorithm doesn't calculate at all when I use more than 300 degrees of freedom (DoF). Another algorithm (sqp) performs the calculation but doesn't work well when I complexify the structure and are out of my test case.

Is there a limit in the number of DoF with lm_feasible algorithm?

If so, how many DoF are maximally possible?

To give an overview and general idea of how the code works:

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

Full example:

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

Fulcrum

can you create a MCVE? – Artur

It probable depends on how much RAM your machine has available, how much is allocated to Octave, etc. I know FEA well, and I understand what optimization means, but I don't know what "minimize the energy of a mechanical structure" means. Is this a static or dynamic problem? What non-geometric variables are allowed to float? Over what range? What do your constraints look like? 300 dof is not a large problem at all. I doubt very much that this is your problem. – Andrews

@Artur I will ASAP. But it may not replicate the issue as the function to optimize will be simpler – Fulcrum

@Andrews This is a static axisymmetric problem. By minimize energy, I mean like my elements are kind of little spring with mass affected by several forces, and all this forces have potential energy (and the spring as well). I search the geometric conformation where the sum of all forces is minimal. For constraints, I have a point pinned to a location (X and Y of first equal initial location), and a sliding point on Y axis (X equal 0 on last point). My geometry is a simple chain of elements – Fulcrum

I know the physics, thank you. I thought that FEA gave you that minimum solution. You start the formulation with a functional integral or weighted residuals. The resulting matrix expression is the minimum energy solution. Optimization suggests something else entirely to me. – Andrews

Optimization in your case would be a step beyond energy minimization, something like "should I use steel, titanium, or bubble gum sugar to fabricate this structure given constraints on max displacement, stress, strain, fatigue life, etc.?" I don't see that in your problem statement. – Andrews

@Andrews Sorry, mean no offense. I optimize the shape and not material because I'm not using a FEA software, I'm doing one. I want to represent how the system act, here: taking the minimal energy shape, not want to design it. The problem is not about how FEA work, but I believe that the issues I go through may be linked to the context and not to an actual limitation of the algorithm/package function, that's why I introduce the context. As there is no "hardwritten" specification described in documentation nor in the script, I think the problem is due to a software problem I don't know of. – Fulcrum

No, the problem is how you are posing the problem. If you were using FEA you'd have to tell the optimizer how to vary the geometry and how to formulate the gradient matrix. Optimization usually means picking a starting point, calculating a "downhill" direction to move in the solution space, and incrementing until you achieve min/max optimum. I see no hint of that in your problem statement or proposed solution. – Andrews

@dufymo That's because I've not put my entire code here. Just enough to showcase the issue. I provide gradient, hessian and starting points to the optimization algorithm in my actual cases, based on potential energies of the forces acting on the mechanical system. And for instance, this particular algorithm can approximate gradient/hessian numericaly. This then go "downhill" toward the minimal energy configuration, by modifying the shape over incrementation. However, as I don't really now how lm_feasible act, I've trouble understanding its limitations – Fulcrum

I don't know the answer to this question, but some ideas might be finding other problems for lm_feasible and testing if they also fail at >300 DoF. Debugging lm_feasible might also help you (if debugging is feasible in lm_feasible). Did you try it? – Aeolotropic

@Trilarion I don't and but I'll try. I begin to think it may not be a problem with the algorithm itself but more with Octave, that may not have enough material ressources while the number of DoF grow, as duffymo point it out on his first comment. I figure out a workaround too using first Sequential Quadratic Programming to get a close initial point, and then using the Levenberg-Marquardt to get the final solution, based on jorgepz answers. Will not answer the question but hopefully solving the actual problem – Fulcrum

Answer 1

Finite Element Method is typically formulated as the optimal criteria for the minimization problem, which is equivalent to the Virtual Work Principle (see books like Hughes of Bathe). The Virtual Work, represents a set of linear (or nonlinear) equations which can be solved more efficiently (with fsolve).

If for some motive you must solve the problem as an optimization problem, then, if you are considering linear elasticity, your strain energy is quadratic, thus you could use the qp Octave function.

To use sparse matrices could also be helpful.

Jacquijacquie

I know that I can use qp to solve the problem (cause I do), and this doesn't really address the question, which is: is there a limit on the number of DoF with lm_feasible. Because I have some non-convergence issue with sqp algorithm, even with many iterations, when my condition begin to be "extreme". With lm_feasible, wich is derivated from Gauss-Newton, I have better and quicker convergence. However, when I use many elements, it doesn't process at all (with lm_feasible). – Fulcrum

Now that you added the functional to minimize I can see that the problem is very nonlinear. Optimization convergence is related to the initial point and the search space size. You should try with an initial point very close to the solution first and see if it converges. If not, then your functional or gradient are wrong. Just to know... Which is the problem u are solving? – Jacquijacquie

I do inflatable structures resolution. My best initial guess is to give the "out of the box"/factory shape and then apply pressure condition. The problem is that in certain pressure cases this initial shape is really different from the final minimal energy one. – Fulcrum

then, you should start validating your algorithm using very small loads. in that case your initial guess is very similar to the optimal, and (if the code is ok) it should converge. if not, you have a problem somewhere. in that case you could then check if you are solving the right problem by minimizing a function with known solution.... solving this is probably something your advisor shall help you to do... – Jacquijacquie