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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序功能：1.生成并输出骨架曲线点至FrameCurve.txt，提取后需要对骨架曲线点按照位移值排序；%默认格式：txt文件为两列，顺序依次为：位移、荷载%读取的曲线默认：从正向开始加载，若为负向开始加载，应先进行调整。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;clc;fid=fopen('FrameCurve.txt','wt');%%创建输出文件FrameCurve.txt，初次写入骨架曲线点[0,0]，应覆盖内容coordinate=[0.00.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负-零fori=3:LineNumif(outdata(i,1)*outdata(i-1,1)=0)&&(outdata(i-1,1)0);%--默认初始加载从正半周开始EndCirclePoints(EndCircle)=i;%分割点的行号EndCirclePointsEndCircle=EndCircle+1;%中间分割点个数为EndCircle-1end;end;%--------------------检查最后是否存在不完整滞回环----------------------ifEndCirclePoints(EndCircle-1)LineNumLoopNum=length(EndCirclePoints)+1;%LoopNum仅为完整循环的个数，应还有最后一个不完整循环elseLoopNum=length(EndCirclePoints);end;%-------------------------分割各个滞回环---------------------------fork=1:LoopNumifk==1LoopCircles(k)={outdata(1:EndCirclePoints(k),:)};elseifkLoopNumLoopCircles(k)={outdata(EndCirclePoints(k-1):EndCirclePoints(k),:)};elseLoopCircles(k)={outdata(EndCirclePoints(k-1):LineNum,:)};%后续步骤中，最后一段应区分是少于半周、少于整周end;end;%--------------------提取各个滞回环荷载骨架曲线点，储存在矩阵FramePointsPostive和FramePointsNegative中-----------------------------fork=1:LoopNumA=LoopCircles{k};ifkLoopNum[ColMaxValue,LineMax]=max(A);%获取各列最大值所在行号LineMax（1*n向量），ColMaxValue为各列的最大值（1*n向量），n为A的列数[ColMinValue,LineMin]=min(A);%获取各列最小值所在行号LineMinFramePointMax=A(LineMax(2),:);%第2列为荷载列FramePointMin=A(LineMin(2),:);FramePointsPostive(k,:)=FramePointMax;%骨架曲线点储存在矩阵FramePoints中FramePointsNegative(k,:)=FramePointMin;%-------------------区分最后一段是完整周，少于半周，还是多于半周---------------elseifEndCirclePoints(EndCircle-1)==LineNum%完整周----------------------------------------------------------------------[ColMaxValue,LineMax]=max(A);%获取各列最大值所在行号LineMax（1*n向量），ColMaxValue为各列的最大值（1*n向量），n为A的列数[ColMinValue,LineMin]=min(A);%获取各列最小值所在行号LineMinFramePointMax=A(LineMax(2),:);%第2列为荷载列FramePointMin=A(LineMin(2),:);FramePointsPostive(k,:)=FramePointMax;%骨架曲线点储存在矩阵FramePoints中FramePointsNegative(k,:)=FramePointMin;elseifoutdata(LineNum,1)=0%仅正半周，仅一个最大值----------------------------------------------------[ColMaxValue,LineMax]=max(A);%获取最大值行号LineMaxFramePointMax=A(LineMax(2),:);FramePointsPostive(k,:)=FramePointMax;else%有负半周，有最大值和最小值----------------------------------------------------[ColMinValue,LineMin]=min(A);%获取最小值行号LineMinFramePointMin=A(LineMin(2),:);FramePointsPostive(k,:)=FramePointMax;FramePointsNegative(k,:)=FramePointMin;end;end;end;%---------------------------------------比较最大位移，删除骨架曲线中同一荷载级的较小骨架点---------------------------------------------------fork=LoopNum:-1:2A=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);%获取各列最小值所在行号LineMinifabs(ColMaxValueA(1)-ColMaxValueB(1))4;%判断为同一荷载级，容差取为4mmifColMaxValueA(2)ColMaxValueB(2);%比较最大荷载值，较小者删去FramePointsPostive(k,:)=[];elseFramePointsPostive(k-1,:)=[];end;end;ifabs(ColMinValueA(1)-ColMinValueB(1))4;ifColMinValueA(2)ColMinValueB(2);%比较最小荷载值，较小者删去FramePointsNegative(k-1,:)=[];elseFramePointsNegative(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')