matlab用fsolve求解非线性方程组为何得到复数的结果

function F=hw(X)
dr1_=X(1);dr2_=X(2);da1_=X(3);da2_=X(4);Fa1=X(5);Fa2=X(6);
W=935;H=480;g=9.81;T=1565;ag=-5.3955;
Fzr=W*g*(1-2*(H*ag)/(T*g))/2%轮胎径向载荷
Fyr=-W*g*(ag/g-2*(H/T)*(ag/g)^2)/2%轮胎轴向载荷
e=-3;Rs=296;Sbc=20;Dpw=56.5;a=35;
S=Sbc+Dpw*tand(a);Dpw%表示钢球节圆直径，Sbc表示钢球中心间距
Fr1=(Fzr*(S/2+e)-Fyr*Rs)/S;%内圈径向载荷
Fr2=(Fzr*(S/2-e)+Fyr*Rs)/S;%外圈径向载荷
fi=0.515;fe=0.525;Dw=11.112;
A=(fi+fe-1)*Dw;%内外沟曲率中心距
D=74;d=40;
dm=(D+d)/2;
r=Dw*cosd(a)/dm;
p11i=2/Dw;p12i=2/Dw;p21i=-1/(fi*Dw);p22i=2*r/(Dw*(1-r));
sumpi=p11i+p12i+p21i+p22i;%sumpi表示接触点主曲率
p11e=2/Dw;p12e=2/Dw;p21e=-1/(fe*Dw);p22e=-2*r/(Dw*(1+r));
sumpe=p11e+p12e+p21e+p22e;%sumpi表示接触点主曲率
ni=0.626;ne=0.688;%与接触点主曲率有关的系数，查表确定
Ki=2.15*10^5*(sumpi)^(-1/2)*(ni)^(-3/2);%负荷变形常量
Ke=2.15*10^5*(sumpe)^(-1/2)*(ne)^(-3/2);%负荷变形常量
Kn=1/(((1/Ki)^(2/3)+(1/Ke)^(2/3))^(3/2));
d0=-0.02;

F(1)=A*da1_+A*da2_+d0;
F(2)=Fa1-Fa2+Fyr;
f1=0;
for  p=0:pi/15:pi
    f1=f1+(cosd(a)+dr1_*cos(p))*cos(p)*max(0,(((sind(a)+da1_)^2+(cosd(a)+dr1_*cos(p))^2)^(0.5)-1)^(1.5))/((sind(a)+da1_)^2+(cosd(a)+dr1_*cos(p))^2)^(0.5);
end 
f2=0;
for   p=0:pi/15:pi
    f2=f2+(cosd(a)+dr2_*cos(p))*cos(p)*max(0,(((sind(a)+da2_)^2+(cosd(a)+dr2_*cos(p))^2)^(0.5)-1)^(1.5))/((sind(a)+da2_)^2+(cosd(a)+dr2_*cos(p))^2)^(0.5);
end 
f3=0;
for   p=0:pi/15:pi
    f3=f3+(sind(a)+da1_)*max(0,(((sind(a)+da1_)^2+(cosd(a)+dr1_*cos(p))^2)^(0.5)-1)^(1.5))/((sind(a)+da1_)^2+(cosd(a)+dr1_*cos(p))^2)^(0.5);
end 
f4=0;
for  p=0:pi/15:pi
    f4=f4+(sind(a)+da2_)*max(0,(((sind(a)+da2_)^2+(cosd(a)+dr2_*cos(p))^2)^(0.5)-1)^(1.5))/((sind(a)+da2_)^2+(cosd(a)+dr2_*cos(p))^2)^(0.5);
end 
F(3)=Fr1-Kn*A^(1.5)*f1;
F(4)=Fr2-Kn*A^(1.5)*f2;
F(5)=Fa1-Kn*A^(1.5)*f3;
F(6)=Fa2-Kn*A^(1.5)*f4;
end

大家好，我想用fsolve求解图片上的六个非线性方程，上面是建立的函数，下面是我用fsolve求解的程序，但是求解的结果是复数，请问是方程的建立出现了问题吗，谢谢。

h=optimset;
h.Display='off';
h.MaxFunEvals=1e+12;%允许的最大函数
h.MaxIter=1000%;允许的最大迭代次数
h.TolX=1e-8;
[jg,fval,exitflag]=fsolve('hw',[0,0,0,0,5000,5000],h) %求解

可能是只有近似解吧？对复数结果用取模方试一下看看数据离标准答案差的打吗

约束多了，不一定有整数解