Error : Not Enough Input Arguments at u1 = z(1)

Hello, I am new to matlab. I wrote this script by going through basics and by looking up answers in this community.

So, I am trying to solve six simultaneous non linear equations with six unknowns. I am using the fsolve method. I am working on parametric surface equations. I am getting the following error, Not Enough Input Arguments at u1 = z(1);

My Questions are:

1) How to solve this issue ?

2) How to get the values of the angles u1, v1, u2 and v2 in either radians or degrees?

function F=myFunction(z)
    syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
    
    %Six Unknowns
    u1 = z(1);
    v1 = z(2);
    u2 = z(3);
    v2 = z(4);
    lambda = z(5);
    t = z(6);
    
    %Define Known Quantities

    r=3; %Radius of the probe tip
    R=5; %Radius of the Hemisphere
    L=50; %Length of the Stylus Stem
    k=1; %Number of the subprocess
    df=0.001; %Deflection along the feed direction

    %Define Surface 1 (Hemispherical Surface) in Parametric Form
    S1=[R*cos(u1)*sin(v1);
        R*sin(u1)*sin(v1);
        R*cos(v1)]; 

    %Define Surface 2(Probe tip) in Parametric Form
    S2=[Cx+r*cos(u2)*sin(v2);
        Cy+r*sin(u2)*sin(v2);
        Cz+r*cos(v2)];

    %Surface 1 Normal

    n1 = cross(diff(S1,u1),diff(S1,v1));

    %Normal Vector at Initial Point of Contact

    n0 = subs(n1,[u1,v1],[deg2rad(0),deg2rad(50)]);

    %Vector Vk in XY Plane
    v0 = n0-[0;0;n0(3)];


    %Initial Position of the Spindle
    S0 = [(r+R)*cos(ui)*sin(vi);
        (r+R)*sin(ui)*sin(vi);
        (r+R)*cos(vi)+L];

    %ui = deg2rad(ui);
    %vi = deg2rad(vi);

    S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(50)]);

    %Position of Spindle at any subprocess k

    Sk = S0 + k*df*n0;


    %Cx = Sk(1)+v0(1)*t;
    %Cy = Sk(2)+v0(2)*t;
    %Cz = Sk(3)-L;


    S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*t,Sk(2)+v0(2)*t,Sk(3)-L]);

    %Surface 2 Normal
    n2 = cross(diff(S2,u2),diff(S2,v2));

    %Six Simultaneous Equations
    F(1) = S1(1)-S2(1);
    F(2) = S1(2)-S2(2);
    F(3) = S1(3)-S2(3);
    F(4) = n1(1)-lambda*n2(1);
    F(5) = n1(2)-lambda*n2(2);
    F(6) = n1(3)-lambda*n2(3);
    
    
    F = [F(1),F(2),F(3),F(4),F(5),F(6)];
    
end

To run this function, I use the following code,

%guess values of six unknowns
zg = [0;0;0;0;1;1];
options = optimoptions('fmincon',...
    'Algorithm','sqp','Display','iter','ConstraintTolerance',1e-12,'StepTolerance',1e-10);
z = fsolve(@myFunction, zg,options);

Thank you, If you need additional information or you think there is something missing please feel free to comment.

Answer 1

[Walter Roberson]

%guess values of six unknowns
zg = [0;0;0;0;1;1];
options = optimoptions('fsolve',...
    'Algorithm','trust-region-dogleg','Display','iter','StepTolerance',1e-10);
z = fsolve(@myFunction, zg,options);

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          7         83.0177                           126               1
     1         14         29.6986       0.395074           83.8               1
     2         21         15.0704              1           33.8               1
     3         22         15.0704            2.5           33.8             2.5
     4         29         6.91011          0.625           11.3           0.625
     5         30         6.91011        1.34722           11.3            1.56
     6         37          3.0559       0.336805           6.06           0.337
     7         44        0.483268       0.678309           9.36           0.842
     8         51     0.000156153      0.0857199          0.267           0.842
     9         58     1.21879e-11    0.000685997       3.04e-05           0.842
    10         65     7.00571e-23    5.24866e-05       2.51e-11           0.842

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

format long g
disp(z)

     -5.24927677027809e-05
      1.55798767520773e-14
     -2.78954536863992e-12
          1.61414015729979
     -7.88445049427304e-15
         0.621028858616472

function F=myFunction(z)
    syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
    syms U1 U2 V1 V2
    
    %Six Unknowns
    u1 = z(1);
    v1 = z(2);
    u2 = z(3);
    v2 = z(4);
    lambda = z(5);
    t = z(6);
    
    %Define Known Quantities

    r=3; %Radius of the probe tip
    R=5; %Radius of the Hemisphere
    L=50; %Length of the Stylus Stem
    k=1; %Number of the subprocess
    df=0.001; %Deflection along the feed direction

    %Define Surface 1 (Hemispherical Surface) in Parametric Form
    S1=[R*cos(U1)*sin(V1);
        R*sin(U1)*sin(V1);
        R*cos(V1)]; 

    %Define Surface 2(Probe tip) in Parametric Form
    S2=[Cx+r*cos(U2)*sin(V2);
        Cy+r*sin(U2)*sin(V2);
        Cz+r*cos(V2)];

    %Surface 1 Normal

    n1 = cross(diff(S1,U1),diff(S1,V1));

    %Normal Vector at Initial Point of Contact

    n0 = subs(n1,[U1,V1],[deg2rad(0),deg2rad(50)]);

    %Vector Vk in XY Plane
    v0 = n0-[0;0;n0(3)];


    %Initial Position of the Spindle
    S0 = [(r+R)*cos(ui)*sin(vi);
        (r+R)*sin(ui)*sin(vi);
        (r+R)*cos(vi)+L];

    %ui = deg2rad(ui);
    %vi = deg2rad(vi);

    S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(50)]);

    %Position of Spindle at any subprocess k

    Sk = S0 + k*df*n0;


    %Cx = Sk(1)+v0(1)*t;
    %Cy = Sk(2)+v0(2)*t;
    %Cz = Sk(3)-L;


    S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*t,Sk(2)+v0(2)*t,Sk(3)-L]);

    %Surface 2 Normal
    N2 = cross(diff(S2,U2),diff(S2,V2));
    n2 = subs(N2, [U2,V2], [u2,v2]);

    %Six Simultaneous Equations
    F(1) = S1(1)-S2(1);
    F(2) = S1(2)-S2(2);
    F(3) = S1(3)-S2(3);
    F(4) = n1(1)-lambda*n2(1);
    F(5) = n1(2)-lambda*n2(2);
    F(6) = n1(3)-lambda*n2(3);
    
    
    F = double(subs([F(1),F(2),F(3),F(4),F(5),F(6)], [U1,U2,V1,V2], [u1, u2, v1, v2]));
    
end

[Muerus Rodrigues]

@Walter Roberson Thank you so much! Thats great!. Solved the error part. However, The values u1, v1, u2 and v2 are angles. I need this code in degrees, can you help me doing that? Also, the values of these 4 angles are unrealistic as seen in the output, is it something to do with the guess values or options ? Thank you.

[Walter Roberson]

In most cases you should not use symbolic expressions inside an fsolve objective function. You should instead call such a function once with symbolic arguements, and then use matlabFunction() to generate a numeric function for use with fsolve.

Your function involves, for example, taking the derivative of an expression with respect to a point. It is not useful to have to calculate the function and its derivative each time: you can be more efficient by building the function and taking its derivative once outside the objective function, and convert it to numeric anonymous function for use inside the objective.

There are some functions that are not defined except in Symbolic Toolbox, and there are some cases where you need to evaluate an expression to a lot of digits in order to get correct ratios, in which case you might need to duck into symbolic work for a few lines. It is, for example, common to take ratios of Bessel functions, but those explode to 10^1000 or 10^-1000 easily, so it just doesn't work to calculate them numerically, in which case Symbolic might be needed.

So it isn't "never" on using symbolic work within an fsolve objective function, but it is a case of "precompute what you can and turn that into numeric work."

[Muerus Rodrigues]

@Walter Roberson Thank you so much for your input!

[Muerus Rodrigues]

@Walter Roberson Is there a way we can define the ranges for the angles u1, v1, u2 and v2. If yes, where should we define in the code?

I want to restrict the solution between the range defined.

[Walter Roberson]

When you use fsolve() it is not possible to put in constraints. However, you can use this adaptation.

You can change lb and ub to put in vectors of bounds.

syms Z [1 6]
myFunctionSym = myFunction(Z);
residueSym = sum(myFunctionSym.^2);

residue = matlabFunction(residueSym, 'vars', {Z.'});

%guess values of six unknowns
zg = [0;0;0;0;1;1];
options = optimoptions('fmincon',...
    'Algorithm','sqp','Display','iter','ConstraintTolerance',1e-12,'StepTolerance',1e-10);

A = []; b = []; Aeq = []; beq = [];
lb = [];   %vector of values of length 6 can be given here
ub = [];   %vector of values of length 6 can be given here
nonlcon = [];
z = fmincon(residue, zg, A, b, Aeq, beq, lb, ub, nonlcon, options);

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           7    8.301768e+01     0.000e+00     1.000e+00     0.000e+00     2.511e+02  
    1          31    4.779081e+01     0.000e+00     2.326e-03     6.285e-01     3.086e+02  
    2          56    3.118665e+01     0.000e+00     1.628e-03     3.633e-01     1.348e+02  
    3          63    1.195154e+01     0.000e+00     1.000e+00     6.509e-01     4.462e+01  
    4          79    1.104509e+01     0.000e+00     4.035e-02     5.199e-01     3.304e+01  
    5          86    1.007965e+01     0.000e+00     1.000e+00     6.052e-01     5.530e+01  
    6          97    9.969204e+00     0.000e+00     2.401e-01     3.000e-01     2.195e+01  
    7         104    9.375079e+00     0.000e+00     1.000e+00     6.860e-02     3.953e+00  
    8         111    9.172119e+00     0.000e+00     1.000e+00     1.143e-01     8.759e+00  
    9         121    9.115529e+00     0.000e+00     3.430e-01     3.450e-01     4.875e+01  
   10         128    8.754702e+00     0.000e+00     1.000e+00     1.557e-01     1.449e+01  
   11         135    8.375483e+00     0.000e+00     1.000e+00     1.157e-01     1.588e+01  
   12         143    7.404024e+00     0.000e+00     7.000e-01     4.827e-01     6.026e+01  
   13         150    5.585885e+00     0.000e+00     1.000e+00     9.869e-02     3.471e+01  
   14         159    4.035767e+00     0.000e+00     4.900e-01     7.522e-01     3.925e+01  
   15         169    2.805451e+00     0.000e+00     3.430e-01     2.414e-01     4.410e+01  
   16         177    1.452774e+00     0.000e+00     7.000e-01     4.167e-01     7.920e+00  
   17         184    3.677735e-01     0.000e+00     1.000e+00     2.103e-01     5.113e+00  
   18         192    2.782488e-01     0.000e+00     7.000e-01     1.879e-01     7.619e+00  
   19         204    1.450481e-01     0.000e+00     1.681e-01     1.755e-01     9.255e+00  
   20         211    7.439211e-03     0.000e+00     1.000e+00     7.640e-02     1.771e+00  
   21         218    6.340881e-04     0.000e+00     1.000e+00     1.194e-02     4.104e-01  
   22         225    1.111051e-06     0.000e+00     1.000e+00     5.995e-03     6.802e-03  
   23         232    4.270553e-09     0.000e+00     1.000e+00     3.107e-04     1.151e-03  
   24         239    2.230559e-12     0.000e+00     1.000e+00     1.685e-05     2.877e-05  
   25         246    5.325381e-14     0.000e+00     1.000e+00     3.223e-07     1.937e-05  
   26         253    2.605770e-14     0.000e+00     1.000e+00     3.132e-08     1.354e-05  
   27         260    1.369715e-14     0.000e+00     1.000e+00     2.422e-08     1.131e-05  
   28         270    1.363514e-14     0.000e+00     3.430e-01     1.153e-09     8.997e-06  
   29         280    1.336503e-14     0.000e+00     3.430e-01     4.078e-09     8.944e-06  
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
   30         287    1.285603e-14     0.000e+00     1.000e+00     1.879e-09     1.090e-05  
   31         297    1.284282e-14     0.000e+00     3.430e-01     1.141e-09     1.086e-05  
   32         306    1.284282e-14     0.000e+00     4.035e-02     1.190e-10     1.086e-05  

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.

format long g
disp(z)

      0.000730239126125599
     -7.78834468623277e-11
       3.7102982464465e-10
           1.6141401468119
      3.62894950956184e-10
         0.621028851319083

function F=myFunction(z)
    syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
    syms U1 U2 V1 V2
    
    %Six Unknowns
    u1 = z(1);
    v1 = z(2);
    u2 = z(3);
    v2 = z(4);
    lambda = z(5);
    t = z(6);
    
    %Define Known Quantities

    r=3; %Radius of the probe tip
    R=5; %Radius of the Hemisphere
    L=50; %Length of the Stylus Stem
    k=1; %Number of the subprocess
    df=0.001; %Deflection along the feed direction

    %Define Surface 1 (Hemispherical Surface) in Parametric Form
    S1=[R*cos(U1)*sin(V1);
        R*sin(U1)*sin(V1);
        R*cos(V1)]; 

    %Define Surface 2(Probe tip) in Parametric Form
    S2=[Cx+r*cos(U2)*sin(V2);
        Cy+r*sin(U2)*sin(V2);
        Cz+r*cos(V2)];

    %Surface 1 Normal

    n1 = cross(diff(S1,U1),diff(S1,V1));

    %Normal Vector at Initial Point of Contact

    n0 = subs(n1,[U1,V1],[deg2rad(0),deg2rad(50)]);

    %Vector Vk in XY Plane
    v0 = n0-[0;0;n0(3)];


    %Initial Position of the Spindle
    S0 = [(r+R)*cos(ui)*sin(vi);
        (r+R)*sin(ui)*sin(vi);
        (r+R)*cos(vi)+L];

    %ui = deg2rad(ui);
    %vi = deg2rad(vi);

    S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(50)]);

    %Position of Spindle at any subprocess k

    Sk = S0 + k*df*n0;


    %Cx = Sk(1)+v0(1)*t;
    %Cy = Sk(2)+v0(2)*t;
    %Cz = Sk(3)-L;


    S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*t,Sk(2)+v0(2)*t,Sk(3)-L]);

    %Surface 2 Normal
    N2 = cross(diff(S2,U2),diff(S2,V2));
    n2 = subs(N2, [U2,V2], [u2,v2]);

    %Six Simultaneous Equations
    F(1) = S1(1)-S2(1);
    F(2) = S1(2)-S2(2);
    F(3) = S1(3)-S2(3);
    F(4) = n1(1)-lambda*n2(1);
    F(5) = n1(2)-lambda*n2(2);
    F(6) = n1(3)-lambda*n2(3);
    
    
    F = (subs([F(1),F(2),F(3),F(4),F(5),F(6)], [U1,U2,V1,V2], [u1, u2, v1, v2]));
    
end

[Muerus Rodrigues]

This is awesome! I was going through one of the video about nonlinear optimization problem using fmincon. Thank you so much again @Walter Roberson

[Muerus Rodrigues]

Thank you for the updated code. Cant thank you enough. Really appreciate your time and valuable feedback.

[Walter Roberson]

Notice that here I used the technique of calling the function once with symbolic arguments to get a symbolic expression, and then turning that into an anonymous handle with matlabFunction() .

[Muerus Rodrigues]

Yes gotcha. Also, I forgot to mention the function F that needs to be minimzed. Overall problem : minimize the distance (I already have the objective function equation in terms of unknowns, which is basically a distance formula) subject to F(1) upto F(6) constraints, and ranges of u1,v1,u2 and v2.

[Walter Roberson]

Ah, you would proceed a bit differently in that case. You would still create

myFunctionSym = myFunction(Z);

but you would then

myFunction_numeric = matlabFunction(myFunctionSym, 'vars',  {Z.'});
nonlcon = @(z) deal(myFunction_numeric(z), []);

z = fmincon(Function_to_minimize, zg, A, b, Aeq, beq, lb, ub, nonlcon, options);

[Muerus Rodrigues]

Perfect. I will try this out. Thanks for all the help.

Answer 2

[Muerus Rodrigues]

@Walter RobersonI made some changes in my code. I have defined Constraints , Objective Function and the Function calls in three different files (Code below). The code executes but the results are incorrect. Is it something to do with the initial guess z0 ? Have I defined lb and ub (angle ranges) correctly ? Thanks!

%Constraint Fuction
function [c,ceq]=constraint(z)
    syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
    syms U1 U2 V1 V2 T;

    %Define Known Quantities
   
    r=3; %Radius of the probe tip
    R=5; %Radius of the Hemisphere
    L=50; %Length of the Stylus Stem
    k=1; %Number of the subprocess
    df=0.001; %Deflection along the feed direction

    %Hemispherical Surface Parametric Equation
    S1=[R*cos(U1)*sin(V1);
    R*sin(U1)*sin(V1);
    R*cos(V1)]; 

    %Normal Vector to Surface 1 at point (U1,V1)
    n1 = cross(diff(S1,U1),diff(S1,V1));

    %Normal Vector to Surface 1 at Initial Known Point given by Azimuthal
    %and Polar Angles
    n0 = subs(n1,[U1,V1],[deg2rad(0),deg2rad(30)]);

    %Vector Vk in XY Plane
    v0 = n0-[0;0;n0(3)];

    %Initial Position of the Spindle at Initial Known Point given by Azimuthal
    %and Polar Angles
    S0 = [(r+R)*cos(ui)*sin(vi);
    (r+R)*sin(ui)*sin(vi);
    ((r+R)*cos(vi))+L];

    S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(30)]);
  
    %Position of Spindle at any subprocess k
 
    Sk = S0 + k*df*n0;
  
    %Define Surface 2(Probe tip) in Parametric Form
    S2=[(Cx+r*cos(U2)*sin(V2));
    (Cy+r*sin(U2)*sin(V2));
    (Cz+r*cos(V2))];

    S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*T,Sk(2)+v0(2)*T,Sk(3)-L]);

    %Surface 2 Normal
    n2 = cross(diff(S2,U2),diff(S2,V2));
          
    u1 = z(1);
    v1 = z(2);
    u2 = z(3);
    v2= z(4);
    lambda = z(5);
    t = z(6);
    
    S1 = double(subs(S1,[U1,V1],[u1,v1]));
    S2 = double(subs(S2,[U2,V2,T],[u2,v2,t]));
    n1 = double(subs(n1,[U1,V1],[u1,v1]));
    n2 = double(subs(n2,[U2,V2,T],[u2,v2,t]));
      
    %inequality Constraint
    c=[];
    %Equality Constraint
    ceq=[S1(1)-S2(1);
        S1(2)-S2(2);
        S1(3)-S2(3);
        n1(1)-lambda*n2(1);
        n1(2)-lambda*n2(2);
        n1(3)-lambda*n2(3)];
end	
-----------------------------------------------------------------------------------------------------
%Objective Function
function f=objfun(z)
    syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
    syms U1 U2 V1 V2 T;

    r=3; %Radius of the probe tip
    R=5; %Radius of the Hemisphere
    L=50; %Length of the Stylus Stem
    k=1; %Number of the subprocess
    df=0.001; %Deflection along the feed direction

    S1=[R*cos(U1)*sin(V1);
    R*sin(U1)*sin(V1);
    R*cos(V1)]; 

    n1 = cross(diff(S1,U1),diff(S1,V1));

    n0 = subs(n1,[U1,V1],[deg2rad(0),deg2rad(30)]);

    %Vector Vk in XY Plane
    v0 = n0-[0;0;n0(3)];

    %Initial Position of the Spindle
    S0 = [(r+R)*cos(ui)*sin(vi);
    (r+R)*sin(ui)*sin(vi);
    ((r+R)*cos(vi))+L];

    %ui = deg2rad(ui);
    %vi = deg2rad(vi);

    S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(30)]);
  
    %Position of Spindle at any subprocess k
 
    Sk = S0 + k*df*n0;
  
    %Define Surface 2(Probe tip) in Parametric Form
    S2=[(Cx+r*cos(U2)*sin(V2));
    (Cy+r*sin(U2)*sin(V2));
    (Cz+r*cos(V2))];

    S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*T,Sk(2)+v0(2)*T,Sk(3)-L]);
 
    u1 = z(1);
    v1 = z(2);
    u2 = z(3);
    v2= z(4);
    %lambda = z(5);
    t = z(6);
    
    S1 = double(subs(S1,[U1,V1],[u1,v1]));
    S2 = double(subs(S2,[U2,V2,T],[u2,v2,t]));      
    %objective function  
    f=sqrt(((S1(1)-S2(1))^2)+((S1(2)-S2(2))^2)+((S1(3)-S2(3))^2));
    %f= 0;
end

-------------------------------------------------------------------------------------------------------
%Function Call

z0=[deg2rad(0),deg2rad(30),deg2rad(0),deg2rad(30),0,0];
lb= [deg2rad(0),deg2rad(0),deg2rad(0),deg2rad(0),-Inf,-Inf];
ub= [deg2rad(360),deg2rad(90),deg2rad(360),deg2rad(180),+Inf,+Inf];
options = optimoptions('fmincon',...
    'Algorithm','sqp','Display','iter','ConstraintTolerance',1e-12,'StepTolerance',1e-10);
[z,fval] = fmincon(@objfun,z0,[],[],[],[],lb,ub,@constraint,options);
disp(z)
% 

[Walter Roberson]

How can you tell that the results are incorrect?

When I execute, I get

4.92017133636611e-12 3.51764035409061e-25 1.05521301865511e-32 2.26415762710886 -1.98782626363788e-24 1.00816899335998

which look reasonable. They are within the ub and lb, and the nonlinear equality constraints are at most 1e-17 . The function value is 4.55219293546628e-17 which seems pretty low.

rad2deg(z(1:4)) says that z(4) is 129.726676185693 degree, which is well within the limit of 180 degrees for it.

[Muerus Rodrigues]

Is there a possibility that the function is minimum at multiple combinations of unknowns ? If yes, how do we find this set of unknowns ? The results I am getting is not what I am looking for even though the solution is correct.

[Walter Roberson]

There are multiple solutions. I used a different software package to find several exact solutions, but along the way I encounted bugs in that other package, so I have to go back and validate the solutions more carefully.

Two solutions that look legitimate are

Z1 = 0
Z2 = atan(1/1917*413437^(1/2)*3^(1/2))
Z3 = pi,
Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi
Z5 = -25/9
Z6 = -1/1000*413437^(1/2)+639/1000

Z1 = pi
Z2 = atan(1/1917*413437^(1/2)*3^(1/2))
Z3 = 0
Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi
Z5 = -25/9
Z6 = 1/1000*413437^(1/2)+639/1000

There is a third solution that looks self-consistent but which I need to cross-check.

The software package found another 7 solutions (10 total), some of which I might be able to modify into something workable, but others look just wrong so I am going to have to go back and work step-by-step to see what I can find.

[Walter Roberson]

Other solutions include

Z1 = anything within the bounds
Z2 = 0
Z3 = 0
Z4 = -atan((-1634563+1022400*3^(1/2))^(1/2)/(-800+639*3^(1/2))) + pi
Z5 = 0
Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2) + 639/1000

Z1 = anything within the bounds
Z2 = 0
Z3 = pi
Z4 = pi - acos(-5/3 + (213*sqrt(3))/160)
Z5 = 0
Z6 = -sqrt(-1634563 + 1022400*sqrt(3))/1000 + 639/1000

[Walter Roberson]

All ten solution families that the other software package found, lead to an objective function value of exactly 0. But I have reason to believe that not all 10 solutions are consistent with the bounds.

[Muerus Rodrigues]

Thank you Walter Robertson, I figured it out why the solutions were incorrect although the objective function was satisfied. It had something to do with the parameters of the surfaces. I assumed the wrong values altogether. It seems to be working now, Thanks much for your constant help.

[Walter Roberson]

Combining the above with some more results:

Pi = pi; arccos = @acos; arctan = @atan;

Z1 = 0, Z2 = 0, Z3 = 0, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000

Z1 = 0, Z2 = 0, Z3 = Pi, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = -1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000

Z1 = 0, Z2 = arctan(1/1917*413437^(1/2)*3^(1/2)), Z3 = pi, Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi, Z5 = -25/9, Z6 = -1/1000*413437^(1/2)+639/1000

Z1 = Pi, Z2 = 0, Z3 = 0, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000

Z1 = Pi, Z2 = 0, Z3 = Pi, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = -1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000

Z1 = pi, Z2 = arctan(1/1917*413437^(1/2)*3^(1/2)), Z3 = 0, Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi, Z5 = -25/9, Z6 = 1/1000*413437^(1/2)+639/1000

Z1 = Z1, Z2 = 0, Z3 = 0, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000

Z1 = Z1, Z2 = 0, Z3 = 0, Z4 = -arctan((-1634563+1022400*3^(1/2))^(1/2)/(-800+639*3^(1/2))) + pi,  Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2) + 639/1000

Z1 = Z1, Z2 = 0, Z3 = Pi, Z4= Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = -1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000

Z1 = Z1, Z2 = 0, Z3 = pi, Z4 = pi - arccos(-5/3 + (213*sqrt(3))/160), Z5 = 0, Z6 = -sqrt(-1634563 + 1022400*sqrt(3))/1000 + 639/1000

All of these lead to an f value of exactly 0, which is the minimum possible for that function.

The notation Z1 = Z1 indicates that Z1 can be any value within its permitted range.