Help using fmincon to find the 2 support locations (x1,x2) that minimize maximum bending moment

1 次查看(过去 30 天)
I am pretty new to Matlab and I feel like this is more of a problem understanding fmincon and function handles than anything. I have found the shear force and bending moment equations for an overhung beam and I need to minimize the maximum bending moment from 0 to L. The equations are shown below. The beam is subject to a uniform load of P0 (and theres a c component too but im not worried about this at the moment) throughout its whole length. L and P0 are given.
My first thoughts are to load the bending moments into the object functions but this is not working. I am looking for some guidance in order to get in the right direction.
syms L P0 c x x1 x2
[sf1,sf2,sf3,bm1,bm2,bm3] = sfbm_(L,P0,c,x,x1,x2);
L = 100;
xl = [0 0];
xu = [L L];
x0 = (xu+xl)/2;
[t,a,w,v] = fmincon([bm1 bm2 bm3],x0,[],[],[],[],xl,xu)
>> [sf1,sf2,sf3,bm1,bm2,bm3] = sfbm_(L,P0,c,x,x1,x2)
sf1 =
- P0*x - (c*x^2)/2
sf2 =
(3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2)/(6*(x1 - x2)) - (c*x^2)/2 - P0*x
sf3 =
((L - x)*(2*P0 + L*c + c*x))/2
bm1 =
x*((c*x^2)/2 + P0*x) - (x^2*(3*P0 + 2*c*x))/6
bm2 =
x*((c*x^2)/2 + P0*x - (3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2)/(6*x1 - 6*x2)) - (x^2*(3*P0 + 2*c*x))/6 + (x1*(3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2))/(6*(x1 - x2))
bm3 =
(x1*(3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2))/(6*(x1 - x2)) - (x2*(3*L^2*P0 + 2*L^3*c - 6*L*P0*x1 - 3*L^2*c*x1))/(6*(x1 - x2)) - (x^2*(3*P0 + 2*c*x))/6 - (x*(L - x)*(2*P0 + L*c + c*x))/2

回答(2 个)

Matt J
Matt J 2022-9-28
编辑:Matt J 2022-9-28
I need to minimize the maximum bending moment from 0 to L
Do you mean the maximum from among [b1 b2 b3]? If so, you probably want to use fminimax instead.
Also, you need non-symbolic fuctions and to provide actual numeric values for the known parameters x,P0, c and L in order to use Optimization Toolbox solvers. So, get rid of syms or use matlabFunction to convert b1,b2,b3 to non-symbolic form.
Also, you should provide a better initial guess than (xl+xu)/2, seeing as how complicated your functions are. Because your function has only two unknowns, it is tractable to do a vectorized sweep over different combinations of xl and xu values from 0 to L to search for an approximate optimum.
  3 个评论
Matt J
Matt J 2022-9-29
编辑:Matt J 2022-9-29
I want to minimize the moment across the entire bar.
But in your post you have a vector quantity [bm1 bm2 bm3]. How is this related to the scalar quantity that you are trying to minimize?
Is there a way to make those a function thorugh the sfbm_ function?
Rewrite sfbm_ to return a numeric value for the moment (for a given x1,x2) rather than symbolic expressions for other stuff.
And I thought the way fmincon works is it takes an initial guess and then uses that initial guess to find other values?
That is true, but your current choice of initial guess looks very arbitrary. Do you have reason to believe that x0 (xu+xl)/2 lies near the optimum? As I mentioned above, it should be easy to do a grid search over x1,x2 to find a better initial guess.
Nathan
Nathan 2022-9-30
After looking more into this I am starting to understand. The function I am trying to pass thorugh is bm1 from 0 to x1, bm2 from x1 to x2 and bm3 from x2 to L. Would it be a good idea to turn this into a piecewise funciton and return that in the funciton?
I understand that I am trying to optimize the values for x1 and x2 so I would need to get values for all of the other paramters. All of them except for x are given, so how would I go about making this a value?

请先登录,再进行评论。


Torsten
Torsten 2022-9-30
编辑:Torsten 2022-10-2
Don't know if the bending moment is correctly implemented. From how I naively interpreted the problem, x1 and x2 should lie symmetrically to L/2. But that's not the result of the optimization.
I took the maximum abs-values for bm1, bm2 and bm3 in the respective intervals.
Maybe you could state the problem more properly for a non-physician.
L = 1;
c = 1;
P0 = 1.5;
A = [1 -1];
b = 0;
lb = [0 0];
ub = [L L];
X0 = [L/3 2*L/3];
[X,fval] = fmincon(@(X)fun(X,L,c,P0),X0,A,b,[],[],lb,ub)
Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.
X = 1×2
0.2325 0.8130
fval = 0.0426
x1 = X(1)
x1 = 0.2325
x2 = X(2)
x2 = 0.8130
function obj = fun(X,L,c,P0)
x1 = X(1);
x2 = X(2);
syms x
bm1 = x*((c*x^2)/2 + P0*x) - (x^2*(3*P0 + 2*c*x))/6;
bm2 = x*((c*x^2)/2 + P0*x - (3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2)/(6*x1 - 6*x2)) - (x^2*(3*P0 + 2*c*x))/6 + (x1*(3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2))/(6*(x1 - x2));
bm3 = (x1*(3*L^2*P0 + 2*L^3*c - 6*L*P0*x2 - 3*L^2*c*x2))/(6*(x1 - x2)) - (x2*(3*L^2*P0 + 2*L^3*c - 6*L*P0*x1 - 3*L^2*c*x1))/(6*(x1 - x2)) - (x^2*(3*P0 + 2*c*x))/6 - (x*(L - x)*(2*P0 + L*c + c*x))/2;
dbm1dx = diff(bm1,x);
dbm2dx = diff(bm2,x);
dbm3dx = diff(bm3,x);
crbm1x = solve(dbm1dx==0);
crbm2x = solve(dbm2dx==0);
crbm3x = solve(dbm3dx==0);
max1 = -Inf;
if double(crbm1x(1)) >=0 && double(crbm1x(1)) <= x1
max1 = max(max1,abs(double(subs(bm1,crbm1x(1)))));
end
if double(crbm1x(2))>=0 && double(crbm1x(2)) <= x1
max1 = max(max1,abs(double(subs(bm1,crbm1x(2)))));
end
max1 = max([max1,abs(double(subs(bm1,x,0))),abs(double(subs(bm1,x,x1)))]);
max2 = -Inf;
if double(crbm2x(1)) >=x1 && double(crbm2x(1)) <= x2
max2 = max(max2,abs(double(subs(bm2,crbm2x(1)))));
end
if double(crbm2x(2))>=x1 && double(crbm2x(2)) <= x2
max2 = max(max2,abs(double(subs(bm2,crbm2x(2)))));
end
max2 = max([max2,abs(double(subs(bm2,x,x1))),abs(double(subs(bm2,x,x2)))]);
max3 = -Inf;
if double(crbm3x(1)) >=x2 && double(crbm3x(1)) <= L
max3 = max(max3,abs(double(subs(bm3,crbm3x(1)))));
end
if double(crbm3x(2))>=x2 && double(crbm3x(2)) <= L
max3 = max(max3,abs(double(subs(bm3,crbm3x(2)))));
end
max3 = max([max3,abs(double(subs(bm3,x,x2))),abs(double(subs(bm3,x,L)))]);
obj = max([max1,max2,max3]);
end

类别

Help CenterFile Exchange 中查找有关 Logical 的更多信息

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by