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

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

Answer 1

[Matt J]

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.

[Nathan]

I want to minimize the moment across the entire bar. So the bounds are bm1~[0, x1), bm2~[x1, x2), bm3~[x2, L] and x is when you go from 0 to L. How would I go about this? Also, P0, L and c are given.

I feel like my class is trying to make me take advantage of fmincon.

Is there a way to make those a function thorugh the sfbm_ function?

And I thought the way fmincon works is it takes an initial guess and then uses that initial guess to find other values?

[Matt J]

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]

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?

Answer 2

[Torsten]

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