clc
clear

%
%

%%肘部法则
%n是样本数，p为特征维数，k为分类数 (列)
[n,p]=size(X);
K=5;D=zeros(K,2);
for k=1:5
    [lable,c,sumd,d]=kmeans(X,k,'dist','sqeuclidean');
    % data，n×p原始数据向量
    % lable，n×1向量，聚类结果标签；
    % c，k×p向量，k个聚类质心的位置
    % sumd，k×1向量，类间所有点与该类质心点距离之和
    % d，n×k向量，每个点与聚类质心的距离
    sse1 = sum(sumd.^2);
    D(k,1) = k;
    D(k,2) = sse1;
end
plot(D(2:end,1),D(2:end,2))
hold on;
plot(D(2:end,1),D(2:end,2),'or');
title('不同K值聚类偏差图')
xlabel('分类数(K值)')
ylabel('簇内误差平方和')
hold off

%

%均值聚类  （行）
x=X.'
rng('default') % For reproducibility
eva1 = evalclusters(x,'kmeans','CalinskiHarabasz','KList',2:4);
eva2 =evalclusters(x,'kmeans','DaviesBouldin','KList',2:4);
eva3 =evalclusters(x,'kmeans','gap','KList',2:4);
eva4 =evalclusters(x,'kmeans','silhouette','KList',2:4);
figure(2)
plot (eva1)  %轮廓系数图
plot (eva2)  %轮廓系数图
plot (eva3)  %轮廓系数图
plot (eva4)  %轮廓系数图
idx = kmeans(x,2);
[s,h] = silhouette(x,idx,'Euclidean')

%

clc
clear

%

%48h 基于混合变量
load('C:\tongbu\BaiduSyncdisk\02PHD\03study\03实验\第一篇论文实验\1.第一次48h\matlab\y1.mat')
load('C:\tongbu\BaiduSyncdisk\02PHD\03study\03实验\第一篇论文实验\1.第一次48h\matlab\u1.mat')
load('C:\tongbu\BaiduSyncdisk\02PHD\03study\03实验\第一篇论文实验\1.第一次48h\matlab\z5.mat')
Ts=5;%采样时间5min
ze1=iddata(y1,u1,Ts);%训练
zv1=iddata(y1,u1,Ts);%验证

V1= arxstruc(zv1(:,:,1),ze1(:,:,1),struc(1:4,1:4,1:4));
nn1 = selstruc(V1,'AIC');
% %第二种输入-输出组合的模型顺序
V2= arxstruc(zv1(:,:,2),ze1(:,:,2),struc(1:4,1:4,1:4));
nn2 = selstruc(V2,'AIC');
% 第3种输入-输出组合的模型顺序
V3= arxstruc(zv1(:,:,3),ze1(:,:,3),struc(1:4,1:4,1:4));
nn3= selstruc(V3,'AIC');

M = arx(ze1,'na',4,'nb',[ 1 2 2 ],'nk',[2 3 1 ]);
compare(ze1,M)%将ARX模型的输出与测量输出进行比较。
present(M);

%


load('C:\tongbu\BaiduSyncdisk\02PHD\03study\03实验\第一篇论文实验\1.第一次48h\matlab\x4.mat')
x1=u1(:,1)
x2=u1(:,2)
x3=u1(:,3)
E=y1(1:4,:)
Et=y1(1:4,:)

%

%输出三个模型对比图
for  t=5:281
    E(t)=1.547*y1(t-1)-0.3602*y1(t-2)-0.3331*y1(t-3)+0.1285*y1(t-4)-0.02217*x1(t-2)+0.4526*x2(t-3)-0.3914*x2(t-4)-0.2327*x3(t-1)+0.1741*x3(t-2);
    Et(t)=1.523*y1(t-1)-0.341*y1(t-2)-0.3603*y1(t-3)+0.152*y1(t-4)-0.004719*x1(t-2)+0.1866*x2(t-3)-0.0988*x2(t-4)+0.6532*x4(t-1)-0.6769*x4(t-2);
end
for  t=1:281
    E_MDL(t)=30.4267+ 1.7547.* x1(t) -1.55.* x2(t) + 5.1559.* x4(t)
end
E_MDL=E_MDL'
figure(1)
plot(z5,y1,z5,E,':',z5,Et,':.',z5,E_MDL,'--','linewidth',1.5)
ylabel('Error[μm]','FontSize',10)
xlabel('Time[min]','FontSize',10)
legend('Measured Values','MV-ARX Model Fitted Values','TV-ARX Model Fitted Values','MLR Model Fitted Values')

%


%输出三个模型残差图
y5=abs(E-y1);
y6=abs(Et-y1);
y7=abs(E_MDL-y1);
figure(2)
plot(z5,y5,z5,y6,':',z5,y7,'--','linewidth',1.5)
ylabel('Residual Error[μm]','FontSize',10)
xlabel('Time[min]','FontSize',10)
legend('MV-ARX Model','TV-ARX Model','MLR Model')

%

R1 = rmse(E,y1);
R2 = rmse(Et,y1);
R3 = rmse(E_MDL,y1);
mae1 = mean(abs(E- y1));
mae2 = mean(abs(Et- y1));
mae3 = mean(abs(E_MDL- y1));

%

ymax=max(y1);
ymin=min(y1);
yc=abs(ymax- ymin);
fit1=(1-R1/yc)*100;
fit2=(1-R2/yc)*100;
fit3=(1-R3/yc)*100;

%

%拟合效果好,非最终使用的
for  t=5:281
    E(t)=0.548*y1(t-1)+0.4606*y1(t-2)-0.32243*x1(t-2)+0.4521*x2(t-3)-0.4912*x2(t-4)-0.2271*x3(t-1)+0.1703*x3(t-2);
end
y5=abs(E-y1)
E_MDL=9.3859 + 2.8510 .* x1 + 1.3213 .* x2 - 1.4585 .* x3
figure(1)
plot(z5,y1,z5,E,':',z5,y5,'--',z5,E_MDL,'linewidth',1.5)
ylabel('Error[μm]','FontSize',10)
xlabel('Time[min]','FontSize',10)
legend('E_arx','Model Fit','Residual','E_MDL')
hold on;

%


%设计成用预测值来拟合，当超过5μm残差时，再用测量值来拟合。
for  t=5:281
    E(t)=0.548*E(t-1)+0.4606*E(t-2)-0.32243*x1(t-2)+0.4521*x2(t-3)-0.4912*x2(t-4)-0.2271*x3(t-1)+0.1703*x3(t-2);
 y5(t)=abs(E(t)-y1(t))
 if y5(t)>=5
     E(t)=0.548*y1(t-1)+0.4606*y1(t-2)-0.32243*x1(t-2)+0.4521*x2(t-3)-0.4912*x2(t-4)-0.2271*x3(t-1)+0.1703*x3(t-2);
 end
end
y5=abs(E-y1)
figure(1)
plot(z5,y1,z5,E,':',z5,y5,'--','linewidth',1.5)
ylabel('Error[μm]','FontSize',10)
xlabel('Time[min]','FontSize',10)
legend('Ex','Model Fit','Residual')
hold off

%

