基于 ANN 进行非线性系统识别附matlab代码

clear; clc; close all;%Input%--------------------------------------------------------------------------%Nodes Coordinates and Degrees of Freedoms%------------------------------------------Node{1}.x=0; Node{1}.y=0;                                 %Node CoordinatesNode{2}.x=5; Node{2}.y=0;Node{3}.x=5; Node{3}.y=5;Node{4}.x=0; Node{4}.y=5;Node{1}.DOFx=1;  Node{1}.DOFy=2;  Node{1}.DOFtheta=3;     %Indices for x, y, and Rotation Degree of Freedoms.  Node{2}.DOFx=4;  Node{2}.DOFy=5;  Node{2}.DOFtheta=6;Node{3}.DOFx=7;  Node{3}.DOFy=8;  Node{3}.DOFtheta=9;Node{4}.DOFx=10; Node{4}.DOFy=11; Node{4}.DOFtheta=12;%Material%-----------------Material{1}.E=200e9;                                      %Modulus of ElasticityMaterial{1}.f_y=400e6;                                    %Stress Beyond Which the Stiffness DecreasesMaterial{1}.Ep=2e9;                                       %Reduced Modulus of ElasticityMaterial{1}.rho=7800;                                     %Material Density%Elements Nodes Connectivity, Area, Moment of Inertia, Modulus of Elasticity, Density, and Rayleigh Damping Parameters%---------------------------------------------------------------------------------------------------------------------Element{1}.Nodes=[1 4];                                   %Element Nodes ConnectivityElement{1}.Material=1;                                    %Element MaterialElement{1}.alpha=0.001;                                   %Element Rayleigh Damping Parameter alphaElement{1}.beta=0.001;                                    %Element Rayleigh Damping Parameter betab1=0.1; h1=0.1;for i=1:1:10Element{1}.Section{i}.x=1/sqrt(3)*b1/2*[-1 -1 1 1];       %Element Section integation points x coordinatesElement{1}.Section{i}.y=1/sqrt(3)*h1/2*[-1 1 -1 1];       %Element Section integation points y coordinatesElement{1}.Section{i}.A=b1*h1/4*[1 1 1 1];                %Element Section integation points areasElement{1}.Section{i}.e=[0 0 0 0];                        %Element Section integation points initial stressesElement{1}.Section{i}.s=[0 0 0 0];                        %Element Section integation points initial strainsElement{1}.Section{i}.k=[0 0 0 0];                        %Element Section integation points initial hardeningendElement{2}.Nodes=[3 2];Element{2}.Material=1; Element{2}.alpha=0;Element{2}.beta=0.00001;b2=0.1; h2=0.1;for i=1:1:10Element{2}.Section{i}.x=1/sqrt(3)*b2/2*[-1 -1 1 1];     Element{2}.Section{i}.y=1/sqrt(3)*h2/2*[-1 1 -1 1];    Element{2}.Section{i}.A=b2*h2/4*[1 1 1 1]; Element{2}.Section{i}.e=[0 0 0 0];Element{2}.Section{i}.s=[0 0 0 0]; Element{2}.Section{i}.k=[0 0 0 0]; endElement{3}.Nodes=[4 3];Element{3}.Material=1; Element{3}.alpha=0;Element{3}.beta=0.00001;b3=0.1; h3=0.1;for i=1:1:10Element{3}.Section{i}.x=1/sqrt(3)*b3/2*[-1 -1 1 1];     Element{3}.Section{i}.y=1/sqrt(3)*h3/2*[-1 1 -1 1];    Element{3}.Section{i}.A=b3*h3/4*[1 1 1 1]; Element{3}.Section{i}.e=[0 0 0 0];Element{3}.Section{i}.s=[0 0 0 0]; Element{3}.Section{i}.k=[0 0 0 0]; end%Support, Support Displacement, and Applied Force%-------------------------------------------------steps=1000;                                               %Number of Analysis StepsSupport=[1 2 3 4 5 6];                                    %Indices of Degrees of Freedom of Supports NSupport=length(Support);U_s=zeros(NSupport,steps);                                %Support Displacement Matrix of Size (Number of Supports, Number of Analysis Steps)                           dU_s=zeros(NSupport,steps);                               %Support Velocity Matrix of Size (Number of Supports, Number of Analysis Steps)    ddU_s=zeros(NSupport,steps);                              %Support Acceleration Matrix of Size (Number of Supports, Number of Analysis Steps)    NNode=length(Node); Force=zeros(3*NNode,steps);                               %Force Matrix of Size (Number of Degrees of Freedoms=3*Number of Nodes, Number of Analysis Steps)Force(10,:)=600.*[1:steps].*sin(1*[1:steps]/10);fs=300;                                                   %Sampling Frequencydelta=10e-4;                                              %Convergance criterion for residual force%Processing%--------------------------------------------------------------------------[Node,Element,F_s,Elements_N]=Analysis(Node,Element,Material,Support,U_s,dU_s,ddU_s,Force,fs,delta);shear=zeros(steps,1);for i=1:1:stepsshear(i)=Elements_N{i}{1}.f(1);end%Output%--------------------------------------------------------------------------Mag_Factor=1;                                             %Magnification Factor for Plotstep=steps;                                               %Step Number to be PlottedPlot_Results(Node,Element,Mag_Factor,0,step); title(strcat('Deflected Shape @ Step=',num2str(step)));Plot_Results(Node,Element,Mag_Factor,1,step); title(strcat('Axial Strain @ Step=',num2str(step)));Plot_Results(Node,Element,Mag_Factor,2,step); title(strcat('Curvature @ Step=',num2str(step)));Plot_Results(Node,Element,Mag_Factor,3,step); title(strcat('Axial Force @ Step=',num2str(step))); Plot_Results(Node,Element,Mag_Factor,4,step); title(strcat('Moment @ Step=',num2str(step))); t=[0:1:steps-1]/fs;figure;subplot(4,1,1);plot(t,Force(10,:));xlabel('Time (s)');ylabel('Force (N)');title('Horizontal Force, DISP, VEL, and ACC at Node 4');subplot(4,1,2);plot(t,Node{4}.Ux);xlabel('Time (s)');ylabel('DISP (m)');subplot(4,1,3);plot(t,Node{4}.dUx);xlabel('Time (s)');ylabel('VEL (m/s)');subplot(4,1,4);plot(t,Node{4}.ddUx);xlabel('Time (s)');ylabel('ACC (m/s^2)');%Identification using neural network%--------------------------------------------------------------------------X=Force(10,:);                         %Input of neural network is forcesT=[Node{4}.Ux*100000 ; shear'];        %Output of neural network is displacement and base shear[X,~] = tonndata(X,1,0);               %Change format of input data[T,~] = tonndata(T,1,0);               %Change format of output datanet = narxnet(0,4,4);                  %Create (nonlinear autoregressive with external input) network with 0 input delays and 4 output delays[Xs,Xi,Ai,Ts] = preparets(net,X,{},T); %Prepares data for trainingnet = train(net,Xs,Ts,Xi,Ai);          %Train the networkview(net)                              %View the networkYopen = net(Xs,Xi,Ai);                 %Predict Outputnet = closeloop(net);                  %convert to closed loop neural networkview(net);                             %view closed-loop network[Xs,Xi,Ai,Ts] = preparets(net,X,{},T); %Prepares data for trainingnet = train(net,Xs,Ts,Xi,Ai);          %Train the networkYclosed = net(Xs,Xi,Ai);               %Predict Output%Plot output with predicted output%--------------------------------------------------------------------------T=cell2mat(T);Yopen=cell2mat(Yopen);Yclosed=cell2mat(Yclosed);figure;plot(t(5:end),T(1,5:end)/100000,'b-',t(5:end),Yopen(1,:)/100000,'r--',t(5:end),Yclosed(1,:)/100000,'k-.'); xlabel('Time (sec)');  ylabel('DISP (m)');legend('Actual','Predicted (openloop)','Predicted (closedloop)');figure;plot(T(1,5:end)/100000,T(2,5:end),'b-',Yopen(1,:)/100000,Yopen(2,:),'r--',Yclosed(1,:)/100000,Yclosed(2,:),'k-.'); xlabel('Interstory Drift (m)');ylabel('Interstory Shear (N)');legend('Actual','Predicted (openloop)','Predicted (closedloop)');