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

4 views (last 30 days)
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

Answers (2)

Matt J
Matt J on 28 Sep 2022
Edited: Matt J on 28 Sep 2022
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 Comments
Nathan
Nathan on 30 Sep 2022
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?

Sign in to comment.


Torsten
Torsten on 30 Sep 2022
Edited: Torsten on 2 Oct 2022
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

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by