Optimize function while meeting a condition

How can i optimize the I function, i want to find the values of h(j) that minimize I, meetentig the condition h(j+1)>h(j)?

Thanks for the help

function [h] = peralte(v,rt,Lc,Lt,s)

lt=0:0.1:Lt;
lc=0:0.1:Lc;

    for i=1:length(lt)
        for j=1:length(lc)
            ht(i,j)=(((v(i,j)).^2)*s)/(rt(j)*9.81);
        end
    end

    for i=1:length(lt)
        for j=1:length(lc)
            I(i,j)=ht(i,j)-h(j);
        end
    end


end

Answer 1

[William Rose]

@Jon Bilbao,

Please explain your function in more detail. It is good practice to include comments in a function eplaining the meaning and size (if not scalars) of all input and output arguments.Is it true that function peralte(v,rt,Lc,Lt,s) returns a vector of length Lt+1? Is v an array? What is its physical meaning? Why does g=9.81 appear?

Your function never assigns a value to h. Therefore it is not complete.

If you want to adjust h to minimize I, then you should define a function I(h), which takes h as an argument and returns I. By the way, I try to avoid using upper case I (eye) or lower case l (ell) as names of functions or variables, since they look the same.

When you have defined I(h), you can try one of matlab's minimization routines. fmincon() is a good place to start. It is not yet clear to me if fmincon() can implement the constraint you specified. But first let's figure out h() and I(h) more clearly.

[William Rose]

@Jon Bilbao,

[Edit function for ht(i,j). I forgot to include the scalras s and 9.81, whcih I have now added.]

You said you want to find the value of vector h that will minimize I. But in the function h=peralte(), I is an array. What part of array I, or what property of I, do you want to minimize?

It would also be helpful if you would provide sample values of the arguments s, Lt, Lc, array v, vector rt, for experimentation.

It appears from the code for [h]= peralte(...) that you expect v to be an array of size lt by lc, and that you expect r to be a vector of length lc. At least those are the ranges of elements read by the code. But lt and lc are funcitons of the arguments Lt and Lc. What happaens if the values of Lt and Lc do not result in values of lt and lc that match the sizes of v and r? An error may occur. Therefore I encourage you to think carefully about the arguments to [h]=peralte(...).

Why do you define vectors lt and lc? You never use their values. You only use their lengths.

You can compute array ht() from array v and vector rt without using Lt or Lc or lt or lc, and without for loops:

v=[10,11; 12, 13];
r=[2 3];
s=9.81;
ht=(s/9.81)*(1./r).*(v.*v)

ht = 2×2

50.0000 40.3333 72.0000 56.3333

Compare to your nested for loops that compute ht(i,j).

[Jon Bilbao]

v is a speed and is calculated with this function, the thin is that Lt is the longitude of a train and Lc the logitude of the curve and with that i calculate the speed of each point of the train in each point.

function [v] = VelTren2(v0,a,Lt,Lc)
%Velocidades de entrada a la curva de los ejes del tren
lt=0:0.1:Lt;
lc=0:0.1:Lc;
v(1)=v0;

    for i=2:length(lt)
        vi(i)=((v0.^2)+2*a*lt(i)).^(1/2);
        for j=1:length(lc)
            v(i,j)=((vi(i).^2)+2*a*lc(j)).^(1/2);
        end
    end

end

[Jon Bilbao]

rt are the curvature radius of the curve o each point lc=0:0.1:Lc;

The final objective of the function is tu calculae the cant (Inclination of the railway to prevent the train from derailing) needed in the curve and its formula is h=v^2*s/r*gravity, there must be only one value of h for each point on the curve, but us you can see there are as many as differents points of the train, so "I" is the diference between the theorical valoue(same sice as v) and the practical(lc=1:0.1:Lc) one.There will be as many I as points of the curve (lc=1:0.1:Lc) so i need a function tha calculates the h tha minimize the "I" values.

I hope that this may help you understanding the question

[Jon Bilbao]

And im sorry but a dont have sample values beacaus this function is for any train and curve

[William Rose]

@Jon Bilbao,

Thank you for the additional explanation.

You said "Lt is the longitude of a train and Lc the logitude of the curve and with that i calculate the speed of each point of the train in each point."

I assume "Lc the logitude of the curve" was a mis-spelling of "Lc the loNgitude of the curve".

Perhaps I am not understanding you correctly. Why is the longitude needed to compute the speed? I understand longitude to be the angle between the prime meridian and a plane through the point in question, perpendicular to the equator. The longitude of Bilbao is -3 degrees.

Perhaps "longitude of a train" means the angle from the beginning of a curve of constant radius to the current train location. Can you provide a diagram showing what Lc and Lt are, and where they are measured from? Is the curve radius constant? Please also provide a diagram that shows the different locaiotns referred to by v(1,1), v(1,2), v(2,1), and v(2.2). I hope such a diagram will help me understand the meaning of i and j in v(i,j).

I assume the units for Lc and Lt are radians. Correct me if I am wrong.

What are the units for I?

You wrote, regarding the angle that will prevent derailment: "its formula is h=v^2*s/r*gravity". What is s, and what are the units for s? If h=angle, then s must dimensionless, for the units to be the same on both sides of the equation.

I think the equation above is not correct, because it does not contain terms for track width and for the height of the center of mass of the car, both of which are key determinants of the angle that is needed to prevent overturning.

[William Rose]

@Jon Bilbao,

Of course the bank angle needed, in order to keep the train's ground reaction force,, perpendicular to the train bed is

but your formula for is not the formula above. You have a factor of s, which I do not understand, and you do not use the inverse tanget formula.

Can you work out the formula for , the bank angle needed to keep the force vector pointed inside the outer rail? Assume the center of mass height above the unbanked track is , and assume the rail gauge is . I get

I could be wrong, so you should check my formula to see if it is correct. I use the max() function so that the formula will not give negative values for bank angle when the velocity is low. I assume each rail is strong enough to withstand any force applied to it.

= vector sum of vertical force (mg) and horizontal force (mv^2/r).

Approximate values for are 2.0 to 2.5 m (see here and here). Standard gauge is 1435 mm, so a reasonable value for is 0.72 m.

In 2015, a train derailed near my home, as it went around a curve. Eight people died and over 200 were injured. The speed limit was 80 kph. The train was going 164 kph when it derailed. I estimated the radius of curvature by adjusting a circle to match the curve, on a screenshot from google maps. See below.

[Jon Bilbao]

s is the width of the train.

The thing is that the value that i want to calculate is the cant (Cant (road/rail) - Wikipedia) of the railway so because of that the formula i use is h=v^2*s/gr. With your formula i understand taht the h calculathed is the position of the trains gravity center.

Answer 2

[William Rose]

@Jon Bilbao, the link to the Wikipedia article is helpful. Now I understand that cant has units of length, not angle.

I cannot assist you if you do not answer the quesitons I have already asked. It will help if you include detailed comments in each function. For example, in function v=velTren2(), please fill in the missing information in the comments below. :

function [v] = VelTren2(v0,a,Lt,Lc)
%Velocidades de entrada a la curva de los ejes del tren
% Inputs
%   v0 = [explain the meaning of v0]  (m/s)
%   a = [explain the meaning of a]  (specify units for a)
%   Lt = [explain meaning of Lt]  (specify units for Lt)
%   Lc = [explain meaning of Lc]  (specify units for Lc)
% Output 
%   v = velocity (NxM array, where N~=10*Lt+1, M~=10*Lc+1) (m/s)

lt=0:0.1:Lt;
lc=0:0.1:Lc;
v(1)=v0;
for i=2:length(lt)
    vi(i)=((v0.^2)+2*a*lt(i)).^(1/2);
    for j=1:length(lc)
        v(i,j)=((vi(i).^2)+2*a*lc(j)).^(1/2);
    end
end
end

Please add similar comments to function h=peralte().

Please answer my previous request to explain, in detail, the significance of i and j, in the array v(i,j), in function velTren2(). What does v(1,2) refer to, compared to v(2,1) or v(2,2) or v(1,1)? A diagram would be helpful.

In function v=velTren2(), you assign a value to v(1). This is confusing, since v is a 2-D array. v(1) will be interpreted by Matlab as a reference to v(1,1). The code never assigns a value to v(1,2) or v(1,3) or v(1,4), etc. Therefore these array elements will be zero when the function returns. Is this what you want?