支持向量机截齿疲劳寿命预测 MATLAB 程序
%************************************************************%%** 子函数——核函数 **
%***********************************************************%
function =ker(x,y,sig)
kerRow=size(x,2);
kerColumn=size(y,2);
omega=zeros(kerRow,kerColumn);
for i=1:kerRow
for j=1:kerColumn
omega(i,j)=norm(x(:,i)-y(:,j));
end
end
K=exp(-omega.^2/(2*sig^2));
%***************************************************************
%** 截割三向力——疲劳——内插预测 **
%***************************************************************
clear;
clc;
close all;
stress_x=;
stress_y=;
stress_z=;
fatigue=;
%数据归一化
stress_x=(stress_x-min(stress_x))/(max(stress_x)-min(stress_x))*0.85+0.15;
stress_y=(stress_y-min(stress_y))/(max(stress_y)-min(stress_y))*0.85+0.15;
stress_z=(stress_z-min(stress_z))/(max(stress_z)-min(stress_z))*0.85+0.15;
E_x=stress_x;
E_y=stress_y;
E_z=stress_z;
E_f=fatigue;
%********学习样本*******************
x1=stress_x(1:2:48);
x2=stress_y(1:2:48);
x3=stress_z(1:2:48);
xInput=;
yOutput=fatigue(1:2:48);
%*********预测样本******************
xs1=stress_x(2:4:48);
xs2=stress_y(2:4:48);
xs3=stress_z(2:4:48);
xsInput=;
yrOutput=fatigue(2:4:48);
sig=sqrt(4.5021);
gamma=1.6236e+008;
alpha0=0.2;
%********构造 RBF 核矩阵***************
ker1=ker(xInput,xInput,sig);
for i=1:size(xInput,2)
ker1(i,i)=ker1(i,i)+1/gamma;
end
%*********LSSVM 模型训练**************
I=ones(1,size(xInput,2));
b=I*inv(ker1)*yOutput'/(I*inv(ker1)*I');
alpha=inv(ker1)*(yOutput'-b*I');
for i=1:size(alpha,1);
if abs(alpha(i))<alpha0;
alpha(i)=0;
end
end
ker2=ker(xInput,xsInput,sig);
ysOutput=alpha'*ker2+b*(ones(1,size(xsInput,2)));
%***********画图*********************
hold on;
plot(1:length(ysOutput),ysOutput,'b-');
plot(1:length(yrOutput),yrOutput,'r-.');
grid;
xlabel('序列');
ylabel('循环载荷次数/次');
title('截齿疲劳寿命曲线');
gtext('—real output');
gtext('.-predict output');
hold off;
%***********计算整体误差*********************
ave=sum(ysOutput)/length(ysOutput);
error1=sqrt(sum((yrOutput-ysOutput).^2)/sum((yrOutput-ave).^2))
%***********计算局部误差*********************
error2=sqrt(((ysOutput-yrOutput)./yrOutput).^2);
=max(error2')
=min(error2')
谢谢你的分享,希望以后今后能见到能多的好资料和大家分享
页:
[1]