vpa solving 5 equations with 5 unknowns, but result is an empty structure

1 回表示 (過去 30 日間)
Raine Wong
Raine Wong 2023 年 5 月 30 日
コメント済み: John D'Errico 2023 年 5 月 31 日
I have 5 unknowns and 5 equations. There is no error or warning when matlab runs it, but for some reason the resulting structure from vpasolve is an empty structure, what should I do? Many parameters were precalculated and set so I made sure there were only 5 unkowns, R alpha t L n. Are the equations too complex? Thanks in advance!
syms R alpha t L n
m_p = 45/1000 %mass of panel in kg, 90/2 for small trap
m_p = 0.0450
I_p = I_cg*1000^2 %moment of inertia abt centroid in kgm2, try chg to mm
Unrecognized function or variable 'I_cg'.
a_p = y_centroid * 1000 %vertical distance from base to centroid in m, try change to mm
% Copper Beryllium from p.2
E = 131000; %N/mm2
poisson_ratio = 0.3;
yield_stress = 1175; %N/mm2
D = E*t^3 / (12*(1-poisson_ratio^2)); %plate bending stiffness
l = L/(R*alpha);
%2.32
M_max_pos = (D*R/t) * (1.152/1000 - 2.21/(1000*l) + (-2.061/10^9 + 7.096/(1000000*l^4))*(R/t)^2)^0.5 * alpha^(2.84 + 18.17/l^2 + (-2.281/1000 + 6.809/(100*l) - 0.245/l^2)*R/t);
%2.33
theta_max_pos = L/R * (1.52/10000 - 4.365/(1000*l) + 1.993/(100*l^2) + (13.45 - 564.5/l^4 + 9.929/(100000*l^8))*(t/R)^2)^0.5 * alpha^(-2.86 + 26.69/l - 81.49/l^2 + (-1.071*10000/l^2 + 4.134*100000/l^4 - 3.74*1000000/l^6)*t/R);
%2.35
M_max_neg = -(D*(2.6*R/(100*t) - 2.143/10^5 * (R/t)^2) * alpha^(2.224+1.338/1000 *R/t));
%2.36
theta_ramp_neg = -(alpha/8 + 18*t/R);
%2.37
M_const_pos = (1+poisson_ratio)*D*alpha;
%2.38
M_const_neg = -(1-poisson_ratio)*D*alpha;
%2.34
theta_heel_pos = 3.673/100*alpha - 4.932/100*alpha^2 + 1.08/100*alpha^3 +t/R*exp(pi) * alpha^(-0.8); %in rad
%6.12
M_heel_pos = M_max_pos * theta_heel_pos / theta_max_pos;
%6.11
theta_max_neg = -abs(M_max_neg)/M_max_pos * theta_max_pos;
%6.23
% M_crit_pos = yield_stress * R^2 * t * (alpha^2 - 4*(1-cos(alpha)) + alpha*sin(alpha))/(2*(2*sin(alpha/2)-alpha*cos(alpha/2)));
%choose arbitraliy
time = 1; %s Harry agreed as well
fn = 1; %Hz
theta_max_0 = 80 *pi/180; %change to rad
%the 5 simulataneous eqs
eq_10 = (0.5*M_max_pos + abs(M_const_neg))*theta_max_pos + (-0.5*M_heel_pos + M_const_pos)*theta_heel_pos + 0.5*((abs(M_max_neg)-abs(M_const_neg))/(abs(theta_max_neg)-abs(theta_ramp_neg)))*(theta_max_pos - abs(theta_ramp_neg))^2 == theta_max_0 * (M_const_pos + abs(M_const_neg));
eq_17 = sqrt(((m_p * (L-L/3+a_p)^2+I_p)*theta_max_0)/(n*D*alpha)) == time; % p1 = 1/3L
eq_18 = sqrt((2*n*M_max_pos)/((m_p*(L+a_p)^2+I_p)*theta_max_pos))/(2*pi) == fn;
%material for 'running check'
%please unpercentage below for second method
eq_23 = yield_stress * R^2 * t * (alpha^2 - 4*(1-cos(alpha)) + alpha*sin(alpha))/(2*(2*sin(alpha/2)-alpha*cos(alpha/2))) >= M_max_pos; %material
eq_25 = E/(yield_stress*(1+poisson_ratio)) <= R/t;
% assume(n,'integer')
sol = vpasolve([eq_10 eq_17 eq_18 eq_23 eq_25],[R alpha t L n])
R_sol = sol.R %in mm
alphaSol = sol.alpha; % in rad
alpha_result = double(alphaSol * 180 / pi) %in deg
tSol = sol.t %in mm
LSol = sol.L %in mm
nSol = sol.n % integer
  2 件のコメント
Alex Sha
Alex Sha 2023 年 5 月 31 日
what are the values of "I_cg" and "y_centroid"?
Raine Wong
Raine Wong 2023 年 5 月 31 日
I calculated them in above sections, y_centroid = 0.0494 and I_cg = 4.0317e-5.

サインインしてコメントする。

回答 (3 件)

Torsten
Torsten 2023 年 5 月 31 日
編集済み: Torsten 2023 年 5 月 31 日
Did you read somewhere that "vpasolve" can handle inequalities ?
Anyway, it seems the problem is too hard for a symbolic solution. I suggest you use "fmincon" with an arbitrary objective function and your equality and inequality constraints defined in "nonlcon".

John D'Errico
John D'Errico 2023 年 5 月 31 日
編集済み: John D'Errico 2023 年 5 月 31 日
Suppose I had two equations. Two unknowns. Totally trivial, but they are not unlike your question.
x + y = 1
x - y <= 3
What is THE numerical solution? Note that vpasolve will not solve this. Nor would solve, unless I add some additional prodding.
syms x y
vpasolve(x + y == 1,x-y <= 3)
ans = struct with fields:
x: [0×1 sym] y: [0×1 sym]
There are of course infinitely many solutions. And in this case, the problem is so totally trivial, that solve will succeed.
sol = solve(x + y == 1,x-y <= 3,'returnconditions',true)
sol = struct with fields:
x: 1 - u y: u parameters: u conditions: -1 <= u
Your problem is exactly the same, at least in concept. For more complex of course. The problem for solve would be that the solution locus for your probem is too complex that solve will also probably fail.
  2 件のコメント
Raine Wong
Raine Wong 2023 年 5 月 31 日
Hello, thanks for your reply. What method would you suggest then?
John D'Errico
John D'Errico 2023 年 5 月 31 日
First, you need to decide what are your goals in this exercise.
Do you want to find the set of all possible solutions? Or do you want to find A solution, ANY numerical solution?
The latter case seems rather boring to me, since there are likely infinitely many solutions, and surely some of those solutions would be meaningful, yet others strange and non-useful. And then you know what happens, if asked to pick one solution between infinitely many choices, a computer will always give you the most completely useless one. In that case, you would need to define some sort of measure of "goodness". Then the problem becomes a well-posed one, where you use an optimizer like fmincon to find the solution that maximizes your measure of goodness, while also satisfying a set of equality constraints, and inequality constraints. Do you see why this approach is now well-posed? It is so, as long as your measure of goodness was built properly.
The first case is then actually far more difficult to solve. In fact, I can give you a long list of mathematical problems where there is no simple solution possible. Hmm. Let me offer an analogy...
Suppose I asked you to find the set of numbers n, such that 2^n-1 is prime. Any such number is rather famously known as a Mersenne prime. There are lists of them that have been tabulated. And we think there are probably infinitely many Mersenne primes, though we do not know that for a fact. There is no simple formula known to conclude that for some value of n, that 2^n-1 is itself prime. (Yes, we do know that n must also be prime.) The point is, the locus of solutions to the problem is just the set of numbers n, such that 2^n-1 is prime. That is both the definition of the problem, and the only valid description of the solution locus.
The above example is not unlike your problem. If you want the "solution" of your problem, the solution is defined as the set of unknowns that satisfy that list of nonlinear equalities and inequalities. There simply is no neat solution, nor any way to find one.
So what method would I suggest? Definitely the first. You tried to use vpasolve as it is, so it would appear you are interested in a numerical solution. So I would spend some time, if need bem some serious thinking time in advance. At a quick glance, it looks like there are 4 equality constraints, and one inequality in what you have provided. And there are 5 unknowns. (I did not look carefully at the equations, and one of them was lengthy, and went off the right side of the window.) But I would first reduce the problem. Pick one of the unknowns, and fix it at some value. Then solve for the other 4 unknowns. Attempt to find all solutions if possible. Are some of them reasonable? Are some of them clearly bogus and useless? Non-meaningful solutions? Do this multiple times for different values of the fixed parameter. Get a feeling for what is happening. I have no such feeling, so I would need to do the above prep work.
Once I had spent some time in understanding what can happen, and what a meaingful, useful anser would look like, I would then spend some time in formulating a measure of goodness for a solution. Finally, I would then use a solver, perhaps FMINCON or perhaps GA to optimize goodness, finding the solution to your problem that was optimally good by that measure.
That is what I would do. I would turn an impossible problem to solve into one that is well-posed.

サインインしてコメントする。


Walter Roberson
Walter Roberson 2023 年 5 月 31 日
When you are dealing with symbolic inequalities, you have two choices:
  1. You can convert them directly to symbolic equalities, solve, and then start probing more carefully to see which directions the desired inequalities hold. You can use vpasolve() for this, but there might be multiple solutions that vpasolve() fails to find.
  2. You can introduce "slack" variables that are constrained as non-negative (for <= and >= ) or as positive (for > or <) and convert to equalities in the slack variable. A <= B ---> A+slack1 == B for slack1>=0 . Now you can solve() -- but not vpasolve() because you want to see how the solutions vary in terms of the slack variable.
The first option, conversion to equalities and then vpasolve(), can be useful in finding some solution to start exploring from. However there might be reasons why the equalities are contradictory but inequalities were not. For example x+y==10, 2x+2y==6 is a contradiction whereas x+y<=10, 2x*2y==6 is possible.

カテゴリ

Help Center および File ExchangeNumeric Solvers についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by