function f=NonLinF(x)
% 本函数求解多齿作用下的转矩优化目标函数
% T0为单齿作用产生的转矩
% x=0-2*pi 摩擦齿轮的旋转角度，w是设计变量
% w(1)-切入角
% w(2)-齿高H
% w(3)-基圆半径
% w(4)-齿距P
% w(5)-齿根宽S
w=0:0.01:2*pi;
i=-4:1:10;
c1=433900; %轮胎径向刚度，单位为N/m
%c2=w3; %摩擦齿轮基圆半径
%c3=w2; %齿顶高度
c4=0.371; %轮胎外径
c5=0.35; %损耗因子，来自参考文献“基于粗糙接触理论的橡胶-金属摩擦副的摩擦分析”
c6=0.4997; %橡胶材料泊松比，来自参考文献“轮胎胶料有限元分析的实验基础及计算”
c7=5000000; %橡胶材料弹性模量，单位Pa，来自参考文献“轮胎胶料有限元分析的实验基础及计算”
c8=3/8; %公式M中的常系数
c9=3/4; %公式M中的常系数
c10=2; %公式M中的常系数
c11=0.5; %公式M中的常系数
e1=2; %公式M中的常指数
e2=1/3; %公式M中的常指数
e3=4/3; %公式M中的常指数
e4=2/3; %公式M中的常指数
e5=0.5; %公式M中的常指数
a3=x(4)/(x(2)+x(3));%齿距p对应的弧度
w0=pi/2-acos(((x(3)+c4)^e1+(x(3)+x(2))^e1-c4^e1)/(c10*(x(3)+c4)*(x(3)+x(2))));%单齿旋转初始角
w1=pi-w0;%单齿旋转结束角度
%a1=w1;%轮齿的齿侧面和齿顶面的夹角，也是摩擦轮胎的切入角
%y1=zeros(size(x));
%y2=zeros(size(x));
%y3=zeros(size(x));
z=zeros(size(w));
N=length(w);
M=length(i);
A=zeros(M,N);
for J=1:M
for k=1:N
if w(k)<w0+i(J)*a3&w(k)>=0;
z(k)=0;
elseif w(k)>=w0+i(J)*a3&w(k)<=w1+i(J)*a3;
a2=atan(((x(3)+x(2))*cos(w(k)-i(J)*a3))./((x(3)+c4)-(x(3)+x(2))*sin(w(k)-i(J)*a3)));%论文图5-10中的beta
c12=((c10*(c4+x(3))-(x(3)+x(2))*sin(w(k)-i(J)*a3)-c4*cos(a2)).^e1+((x(3)+x(2))*cos(w(k)-i(J)*a3)-c4*sin(a2)).^e1).^e5;%轮胎的径向变形量
p1=c1*c11*c12.*sin(c10*(x(1)-w(k)+i(J)*a3+a2));
p2=c8*c5*((c9*((1-c6)^e1)/c7)^e2);
p3=((c1*c12.*sin(w(k)-i(J)*a3-a2)).^e3)/((x(3)+x(2))^e4);
p4=cos(pi/2-w(k)+i(J)*a3+a2);
p5=c4-c12;
z(k)=-((p1+p2*p3.*p4).*p5);
else w(k)<=(2*pi)&w(k)>w1+i(J)*a3;
z(k)=0;
end
end
B(J,
=z;
end
y=B;
sumT=zeros(size(w));
for J=1:M
sumT=sumT+y(J,
;
end
WsumT=sumT.*W(w);
indxx=find(w>w0&w<w1);
Tsum=WsumT(indxx);
minT=min(Tsum);
maxT=max(Tsum);
f=maxT-minT;