LaGrange Multiplier constrain with discretised functions

Hello!

I am a little stuck with my optimization assignment. It's, no surprise, the catenary but the approach we are looking for is a little different than whatever I could track down on different forums across Google. Or maybe I just don't get something and can't incorporate it correctly into my code.

Background

The first part was solving the unconstrained problem (soap film minimizing area), and the second part is the fixed length chain case. I am not able to put the whole script in here as it involves a few different functions, and also I need help with only this one specific idea and not sorting out the whole project for, so I will try to explain this as well as possible based on only part of the code - based on threads I've enountered in here already (https://se.mathworks.com/matlabcentral/answers/531298-finding-minimum-maximum-of-function-using-lagrange-multipliers and https://se.mathworks.com/matlabcentral/answers/436044-is-there-anyone-to-help-me-for-lagrange-multipliers-method) I see people who really know that stuff and I also have this violent presumption that I am VERY close to figuring that out, I'm just missing something important.

Intro

Let's start from the unconstrained input so it's easier to point out the differences:

function y=F_Discrete(U, N)
 
%Input
X = 3;              %Separation between rings
delX = X/(N-1);     %Segments definition
U(1)   = 2;         %Start height = first ring radius
U(N)   = 3;         %End height = second ring radius
 
f_1 = (U(1)*sqrt(1+((U(2)-U(1))/(delX))^2));    %First sum component
f_N = (U(N)*sqrt(1+((U(N)-U(N-1))/(delX))^2));  %Last sum component
f = 0;                                          %Rest of sum preallocation
 
for i=2:(N-1)
f = f + (U(i)*sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2));    %Second to (N-1)th component
end
 
%Sum
y = delX * (1/2) * (f_1 + 2*f + f_N);

Nothing fancy here, it's basically taking the Area integral and using the trapezoidal rule on it:

I even have a slide from my report to show steps that were behind this one. It works perfectly fine

Problem

Now is the part I am stuck with. The new code I came up with looks like that:

function y=F_Discrete(U, N)

%Input
X = 4;              %Separation
delX = X/(N-1);     %Segments definition
U(1)   = 2;         %Start height 
U(N)   = 2;         %End height 
L = 6;              %Cable length

%Energy
f_1 = (U(1)*sqrt(1+((U(2)-U(1))/(delX))^2));    %First sum component
f_N = (U(N)*sqrt(1+((U(N)-U(N-1))/(delX))^2));  %Last sum component
f = 0;                                          %Rest of sum preallocation

for i=2:(N-1)
f = f + (U(i)*sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2));    %Second to (N-1)th component
end

E = delX * (1/2) * (f_1 + 2*f + f_N);

%Length
g_1 = (sqrt(1+((U(2)-U(1))/(delX))^2));    %First sum component
g_N = (sqrt(1+((U(N)-U(N-1))/(delX))^2));  %Last sum component
g = 0;                                     %Rest of sum preallocation

for i=2:(N-1)
g = g + (sqrt(1+((U(i+1)-U(i-1))/(2*delX))^2));    %Second to (N-1)th component
end

syms lambda
Len = lambda*(delX * (1/2) * (g_1 + 2*g + g_N)) - L == 0

F = E + Len == 0;

%Sum
% y = E + lambda*Len
% y = E + Len
% y = E + Len;
y = E;

Don't mind this y = E; in the end, this is only so the whole program runs through (this in only the input for the bigger optimization 'calculator' with different gradient-like methods, in the input N is the number of segments, and U is the value for that specific point that is calculated outside through BFGS or other Conjugated Gradient function), and so I can check how Len or lambda, or whatever I want to check, looks like.

As You can see, I already tried to incorporate the Multiplier in some way here, I also tried a million other things too, deriving it here and there, multiplying different things in different places, but nothing got me even remotely close to getting the values close to the analytical ones. I am sure I just don't understand the concept correctly and my restriction doesn't apply to the value in the way I want it.

Here's where I need help: since I already have the integrated values from the trapezoidal approach, how do I apply the length constrain there? Where to put it, am I thinking correctly on how to include it in this 'Len'?

My deadline is until the end of the week but I still have a few things to work on there, this is just a bump in the road that totally locks me out! I'm losing my mind a little, please help! :)

Answer 1

[Alan Weiss]

I think that you are making a mistake by trying to use symbolic variables here.

If I understand what you are trying to do, you have a bunch of fixed-length segments that you want to place so that the segments are all connected, the end of one to the beginning of the next, and you want to minimize the total gravitational potential energy of the weights, where I suppose that the weights are at segment ends or in the middle or something.

I think that you should use 2-D variables to represent the x- and y-coordinates of each segment end. Or, if you like, use an angle that represents the angle of each segment from the hoorizontal or from the preceding segment. Anything you like. Then write the equations that represent the constraints and that represent the total potential energy that you are trying to minimize. Then minimize the objective function subject to the constraints using a standard soolver such as fmincon.

Good luck,

[Kamil Pospiech]

Hello Alan!

Thanks for taking a look!

Interesting, I didn't think of it that way, I need to figure out if changing the frame of reference like that is going to work. Up until now, it was a little different approach:

As You can see in the unconstrained example, my 'F_Discrete' takes two input arguments: number of segments N and U, which is another function (one that is being minimized externally by gradient like algorithms written step by step - we don't use the toolbox in here). The code then calculates the integral of the Area by using trapezoid sum and gives its value as y. This goes into the optimization alogorithm that tries to minimize it and comes back with the new U etc. etc.

I am 'discretizing the grid' sort of - delX is dividing the distance between the start and end point so the points are always in the same disctance from each other in X but the length of the segments themselves are not equal. Based on this X coordinate I am making guesses for the Y values and minimizing the function based on this output.

Now You made me think if switching into working on fixed segment length could work better - I'll see if I can up with something bright but nothing pops out on me straight away. Good idea, I need to process it now :)