vpasolve More Unknowns than Equations => No Solutions

Question

I have 3 equations and 4 unknowns that I want to solve for. Of course there are solutions, but vpasolve doesnt return any, if I drop down to 3 eq and 3 unknowns, it works good. I know that with more unknowns I have pretty much infinite number of solutions, so how do I make it work in that case?

Here is my code:

syms x y z theta1 theta2 theta3 phi1

xEquation = x == cos(theta1)*cos(phi1) + cos(theta1 + theta2)*cos(phi1) + cos(theta1 + theta2 + theta3)*cos(phi1)
yEquation = y == cos(theta1)*sin(phi1) + cos(theta1 + theta2)*sin(phi1) + cos(theta1 + theta2 + theta3)*sin(phi1)
zEquation = z == sin(theta1) + sin(theta1 + theta2) + sin(theta1 + theta2 + theta3)

x = 2;
y = 1.5;
z = 0;
sol = vpasolve([eval(xEquation), eval(yEquation), eval(zEquation)], [theta1, theta2, theta3, phi1], [-pi pi; -pi pi; -pi pi; -pi pi;]);

That produces sol struct with 4 fields but they are empty, no solutions.

Answer 1

Solving m equations with n unknowns such as m<n, means some variables will be parameters that's depending on other variables.

For example

x-3y+z = 2
x-y+5z = 5

Let suppose z to be the parameter
To solve this in Matlab here is the code

syms x y z

eq1 = x-3*y+z ==2;
eq2 = x-y+5*z ==5;
sol = solve(eq1, eq2, x, y);
sol.x
sol.y

As you can see z been omitting in solve expression, which means it will be considered as the parameter
The solution is

sol.x = 13/2 - 7*z
sol.y = 3/2 - 2*z

• From the solution we can see that x, y and z are not numerical values, so that's the reason you can not use vpasolvewhich stands for Variable-precision arithmetic. vpasolve gives a precision to numerical solution
• Also since x, y and z are not numerical values you can not predefine a range for them, unless you fix z first

---

You may use solve to have a view of the solution set, here I consider phi1 to the parameter, so omitted it in solve expression

syms x y z theta1 theta2 theta3 phi1
    
xEquation = 2 == cos(theta1)*cos(phi1) + cos(theta1 + theta2)*cos(phi1) + cos(theta1 + theta2 + theta3)*cos(phi1);
yEquation = 1.5 == cos(theta1)*sin(phi1) + cos(theta1 + theta2)*sin(phi1) + cos(theta1 + theta2 + theta3)*sin(phi1);
zEquation = 0 == sin(theta1) + sin(theta1 + theta2) + sin(theta1 + theta2 + theta3);
  
    
sol = solve(xEquation, yEquation, zEquation, theta1, theta2, theta3);

Solution set

sol.theta1 = [-pi*(2*n - 1); pi*(2*m + 1); pi*k; pi*k]
sol.theta2 = [pi*(2*m + 1);  -pi*(2*n - 1);  -pi*(2*n - 1); z]
sol.theta3 = [pi*k; pi*k; pi*(2*m + 1); pi*(2*m + 1)] 
phi1 is the parameter 

First solution set

X = [sol.theta1(1); sol.theta2(1); sol.theta3(1); phi1]
X = [-pi*(2*n - 1); pi*(2*m + 1); pi*k; phi1]

There are 4 sets written accordingly to the above deductions

• z is the parameter, k, m, n are integer numbers which are mainly used for trigonometric functions periodicity
  

• If you set z in range [-pi, pi] you can adjust k, m and n to get valid solution in range [-pi, pi].

As you can see it's time costing

---

Use fmincon

• Alternatively you can use fmincon to solve your problem
• I mainly minimize a constant function, fmincon will just give solutions that satisfy the constraints, here ceq = 0
• The interval of search [-pi pi] is transformed to lb = -pi and ub = pi
• I set the initial guess to 0

The code is as follow

t = 0:0.1:1;
x = 1.5 + 0.5 .* cos(8 .* pi .* t);
y = 1.5 + 0.5 .* sin(8 .* pi .* t); 
z = 1 .* t .* ones(size(x));

lb = -pi*ones(1, 4);
ub = -lb;

p0 = zeros(1,4);
sol = cell(1,length(t));

for i = 1:length(t)
    sol{i} = fmincon(@(p)0,p0,[],[],[],[],lb,ub,@(p)nonlincon(x(i),y(i), z(i), p(1), p(2), p(3), p(4)));


end


function [c, ceq] = nonlincon(x,y, z, theta1, theta2, theta3, phi1)

    c = [];
    ceq(1) = cos(theta1)*cos(phi1) + cos(theta1 + theta2)*cos(phi1) + cos(theta1 + theta2 + theta3)*cos(phi1)-x;
    ceq(2) = cos(theta1)*sin(phi1) + cos(theta1 + theta2)*sin(phi1) + cos(theta1 + theta2 + theta3)*sin(phi1)-y;
    ceq(3) = sin(theta1) + sin(theta1 + theta2) + sin(theta1 + theta2 + theta3)-z;
end

Solution

Second set of solution when t = 0.1 is sol{2}

sol{2}.(1) = pheta1
sol{2}.(2) = pheta2
sol{2}.(3) = pheta3
sol{2}.(4) = phi1

You can follow same logic to find solution at different time t

The Entire solution