I have an optimization problem where I want to optimize the thicknesses t1, t2, t3 by solving them symbolically. The solutions should be approximately:
t1 = 1.72053;
t2 = 1.91395;
t3 = 0.781725;
%Given
M_y = 2000000; %[Nmm]
b = 100; %[mm]
E = 70000; %[N/mm^2]
k_Druck = 3.62;
k_Biegung = 21.7;
sigma_zul = 100; %[N/mm^2]
syms t1 t2 t3 positive
% center of gravity
zGes =(b^2*t1+b^2*t3)/(t1*b+t2*b+2*t3*b);
zSeite = zGes - b/2;
zUnten = b - zGes;
zOben= -zGes;
%area moment of inertia
I_y = 2*((t3*b^3/12)+zSeite^2*t3*b)+ ((b*t1^3/12)+zUnten^2*t1*b) + ((b*t2^3/12)+zOben^2*t2*b);
sigma_Oben = -M_y / I_y * zOben; %Stabilitätsversagen (druckbelastete Oberseite)
sigma_Unten = M_y / I_y * zUnten; %Festigkeitsversagen (zugbelasteten Unterseite)
sigma_krit_Oben =k_Druck * E * (t2/b)^2; %kritische Spannung Oberseite
sigma_krit_Seiten =k_Biegung * E * (t3/b)^2; %kritische Spannung Seitenwand
%% condition
Beulen_Oben = sigma_Oben == sigma_krit_Oben;
Beulen_Seiten= sigma_Oben == sigma_krit_Seiten;
Festig_Unten = sigma_Unten==sigma_zul;
%Solver
eqns = [Beulen_Oben,Beulen_Seiten,Festig_Unten];
vars = [t1, t2, t3];
[t1, t2, t3]=vpasolve(eqns,vars)
I tried to solve the problem using the vpasolve function in Matlab, but with the positive condition, the function always outputs Empty sym: 0-by-1. When I input the given solutions for the thicknesses the given conditions are nearly satisfied. I assume them to be a bit off due to the number of fractional digits that I input.
There may or may not be an actual solution (root) for this system in the mathematical sense. vpasolve can return an empty solution when it fails to find a solution, as is described in the documentation. Sometimes this can be resolved by providing bounds on an initial guess as an optional third argument. This does not always work, as in this case. vpasolve will also have trouble in cases where the equation has a minimum that just touches zero, but does not cross or almost crosses zero, but not quite. In your case, this may just be due to the specific combination or parameter values that you have selected.
I would suggest solving this numerically via fsolve instead. You can use your existing symbolic math code and convert your equations to a numerical function using matlabFunction. The first step is to move the righthand side of your equality equations to the lefthand side using the lhs and rhs functions:
eqns0 = lhs(eqns)-rhs(eqns);
Then you can convert this to a numeric anonymous function with vector input as expected by fsolve:
f = matlabFunction(eqns0, 'Vars', {vars});
Then select an initial guess and call fsolve:
x0 = [1.7 1.9 0.8]; % Initial guess
[x, fval] = fsolve(f, x0)
On my R2023a system, this returns:
x =
1.7205 1.9139 0.7817
fval =
1.0e-12 *
0.0284 -0.1137 0.0426
For this system of equations, the solution does not appear to be very sensitive to the initial guess, but it is possible that there may be more than one solution. See the documentation for how to adjust the convergence criteria to reduce the residual fval closer to zero.