从滞回曲线提取骨架曲线点Matlab程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能：1.生成并输出骨架曲线点至FrameCurve.txt，提取后需要对骨架曲线点按照位移值排序；
%默认格式：txt文件为两列，顺序依次为：<位移>、<荷载>
%读取的曲线默认：从正向开始加载，若为负向开始加载，应先进行调整。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;
fid=fopen('FrameCurve.txt','wt'); %%创建输出文件FrameCurve.txt，初次写入骨架曲线点[0,0]，应覆盖内容
coordinate=[0.0 0.0];
dlmwrite('FrameCurve.txt',coordinate,'delimiter','\t','newline','pc','precision','%.2f')
sta=fclose(fid);
outdata=load('jgz3.txt'); %%读取文件数据，默认是为<两列>,<位移，荷载>，outdata矩阵,
LineNum=size(outdata,1);
EndCircle=1;
%--------------------找各级循环分割点位置----------------------负->正or负->零
for i=3:LineNum
if (outdata(i,1)*outdata(i-1,1)<=0)&&(outdata(i-1,1)<0);%--默认初始加载从正半周开始
EndCirclePoints(EndCircle)=i; %分割点的行号EndCirclePoints
EndCircle=EndCircle+1; %中间分割点个数为EndCircle-1
end;
end;

%--------------------检查最后是否存在不完整滞回环----------------------
if EndCirclePoints(EndCircle-1)LoopNum=length(EndCirclePoints)+1; %LoopNum仅为完整循环的个数，应还有最后一个不完整循环
else
LoopNum=length(EndCirclePoints);
end;
%-------------------------分割各个滞回环---------------------------
for k=1:LoopNum
if k==1
LoopCircles(k)={outdata(1:EndCirclePoints(k),:)};
elseif kLoopCircles(k)={outdata(EndCirclePoints(k-1):EndCirclePoints(k),:)};
else
LoopCircles(k)={outdata(EndCirclePoints(k-1):LineNum,:)};%后续步骤中，最后一段应区分是少于半周、少于整周
end;
end;
%--------------------提取各个滞回环荷载骨架曲线点，储存在矩阵FramePointsPostive和FramePointsNegative中-----------------------------
for k=1:LoopNum
A=LoopCircles{k};
if k[ColMaxValue,LineMax]=max(A); %获取各列最大值所在行号LineMax（1*n向量），ColMaxValue为各列的最大值（1*n向量），n为A的列数
[ColMinValue,LineMin]=min(A); %获取各列最小值所在行号LineMin
FramePointMax=A(LineMax(2),:); %第2列为荷载列
FramePointMin=A(LineMin(2),:);
FramePointsPostive(k,:)=FramePointMax; %骨架曲线点储存在矩阵FramePoints中
FramePointsNegative(k,:)=FramePointMin;
%-------------------区分最后一段是完整周，少于半周，还是多于半周---------------
else
if EndCirclePoints(EndCircle-1)==LineNum %完整周-----------------------------------------

-----------------------------
[ColMaxValue,LineMax]=max(A); %获取各列最大值所在行号LineMax（1*n向量），ColMaxValue为各列的最大值（1*n向量），n为A的列数
[ColMinValue,LineMin]=min(A); %获取各列最小值所在行号LineMin
FramePointMax=A(LineMax(2),:); %第2列为荷载列
FramePointMin=A(LineMin(2),:);
FramePointsPostive(k,:)=FramePointMax; %骨架曲线点储存在矩阵FramePoints中
FramePointsNegative(k,:)=FramePointMin;
elseif outdata(LineNum,1)>=0 %仅正半周，仅一个最大值----------------------------------------------------
[ColMaxValue,LineMax]=max(A); %获取最大值行号LineMax
FramePointMax=A(LineMax(2),:);
FramePointsPostive(k,:)=FramePointMax;
else %有负半周，有最大值和最小值----------------------------------------------------
[ColMinValue,LineMin]=min(A); %获取最小值行号LineMin
FramePointMin=A(LineMin(2),:);
FramePointsPostive(k,:)=FramePointMax;
FramePointsNegative(k,:)=FramePointMin;
end;
end;
end;
%---------------------------------------比较最大位移，删除骨架曲线中同一荷载级的较小骨架点---------------------------------------------------
for k=LoopNum:-1:2
A=LoopCircles{k};
B=LoopCircles{k-1};
[ColMaxValueA,LineMaxA]=max(A); %获取各列最大值所在行号LineMax（1*n向量），ColMaxValue为各列的最大值（1*n向量），n为A的列数
[ColMinValueA,LineMinA]=min(A); %获取各列最小值所在行号LineMin
[ColMaxValueB,LineMaxB]=max(B); %获取各列最大值所在行号LineMax（1*n向量），ColMaxValue为各列的最大值（1*n向量），n为A的列数
[ColMinValueB,LineMinB]=min(B); %获取各列最小值所在行号LineMin
if abs(ColMaxValueA(1)-ColMaxValueB(1))<4; %判断为同一荷载级，容差取为4mm
if ColMaxValueA(2)FramePointsPostive(k,:)=[];
else
FramePointsPostive(k-1,:)=[];
end;
end;
if abs(ColMinValueA(1)-ColMinValueB(1))<4;
if ColMinValueA(2)FramePointsNegative(k-1,:)=[];
else
FramePointsNegative(k,:)=[];
end;
end;
end;
dlmwrite('FrameCurve.txt',FramePointsPostive,'-append','delimiter','\t','newline','pc','precision','%.2f')
dlmwrite('FrameCurve.txt',FramePointsNegative,'-append','delimiter','\t','newline','pc','precision','%.2f')