clearvars;close all;clc
clf;close all;


Nyy = 446; % Number of nodes to discretize the bridge structure. We need a spatial resolution of 1 m in our case.
Nmodes = 6; % number of modes for each degree of freedoms, i.e. a total of 3x6 =18 nodes
[Bridge] = LysefjordBridge(Nyy); % bridge structural properties


[wn,phi] = eigenBridge(Bridge,Nmodes); % get the eigen frequencies and mode shapes [2,3]
fn=wn/(2*pi); % get the frequencies in Hz
zetaStruct = 5e-3*ones(3,Nmodes); % modal damping ratios




% Define the structure bridge
Bridge.zetaStruct = zetaStruct;
Bridge.wn = 2*pi*fn;
Bridge.phi = phi;

Sampling parameters

fs = 15; % sampling frequency
M = 15; %  M must be a power of 2
[t,f] = getSamplingPara(M,fs);
N=numel(t);

Wind turbulence

% No wind in this example
tmax = t(end);
Wind.u = zeros(Nyy,N);
Wind.w = zeros(Nyy,N);
Wind.t = t;

Computation of the vehicle-induced response at midspan (target data)

yTarget = Bridge.L/2; % midspan
[~,indY]=min(abs(Bridge.x.*Bridge.L-yTarget)); % get indice where the midspan is located 
Nvehicle = 4; % Number of vehicle selected
tic
%Create the structure variable vehicle_Target that contains the aprameters
%of each vehicle
rng(1) % To ensure reproducibility of this example
clear vehicleTarget
for ii=1:Nvehicle
    vehicleTarget(ii).mass = randi([1e3 50e3],1); % mass of the vehicle (kg)
    vehicleTarget(ii).speed = randi([30 80],1)/3.6; % speed of the vehicle (m/s)
    vehicleTarget(ii).direction = 1; % All the vehicle have the same direction. If "-1", the vehicle is moving in the opposite direction
    vehicleTarget(ii).tStart = 55+ii.*400; % arrival time (sec). t = 0s is the initial time of the simulation.  
end


% Compute the wind and vehicle-induced bridge response
% Newmark-beta algorithm is used for the solver
[Do1,Ao1] = dynaResp_vehicle_TD(Bridge,Wind,'vehicle',vehicleTarget); 
% Get the vertical bridge displacement response at midspan
rz = squeeze(Do1(2,indY,:));
% Some white noise is added to challenge the identification algorithm
rz = rz + 0.03*std(rz).*randn(size(rz));
toc

Visualization of the vehicle-induced response at midspan

gcf;close all;
figure
plot(t,rz,'r')
xlabel('Time (s)')
ylabel('r_z (m)')
set(gcf,'color','w')
axis tight

Outlier-detection algorithm and clustering algorithm

tSim = t;
minNp = 30; % minimum number of cluster points
CO = 30*fs; % cutoff value for the clusters, constructed using a hierarchical cluster tree (cf Matlab function "cluster")
fMin = 0.04; % cut-off frequency (high pass filter) -> must be above the lower frequency correctly measured by the accelerometer
fMax = 0.15; % cut-off frequency (low pass filter) -> must be lower than the first eigen frequency. here the lowest vertical eigenfrequency is 0.21 Hz
[rz_filtered,t_filtered,tImpact_guess,tCluster,rzCluster] = findVehicleID(rz,t,...
    'minNp',minNp,'CuttOff',CO,'deltaTmax',inf,'fcLow',fMin,'fcUp',fMax);
% tImpact_guess is the Initial guess estimate based on the cluster analysis
% rz_filtered,t_filtered, tCluster and rzCluster are used only for visualization purpose!


Ncluster = numel(tCluster);


clf;close all;
figure
COLOR = {'r','b','g','m','c','y'};
if Ncluster>numel(COLOR), COLOR = repmat(COLOR,1,Ncluster); end
for ii=1:Ncluster
    subplot(Ncluster,1,ii)
    plot(t_filtered,rz_filtered,'k')
    hold on
    box on
    plot(tCluster{ii},rzCluster{ii},'o','color',COLOR{ii},'markersize',2)
    xlabel('time (s)')
    ylabel('r_z (m)')
    axis tight
end
set(gcf,'color','w')

Identification of the vehicles' arrival time and speed

gradS = 0.3; % termination tolerance for the speed estimate, i.e. the speed is identified at +/- gradS (in m/s)
gradT = 0.5; % termination tolerance for the arrival time estimate, i.e. it is identified at +/- gradT (in s)
newN = 150; % number of data point for the time-domain simulation: less points means faster algorithm but also less accurate one
posAcc = yTarget./Bridge.L; % relative position of the accelerometer
tic
[vehicleSpeed,tImpact] = findSpeed(Bridge,t,rz,posAcc,tImpact_guess,...
    'gradS',gradS,'gradS',gradT,'newN',newN,'fcLow',fMin,'fcUp',fMax);
toc
fprintf('Target speed (km/h): \n')
fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[vehicleTarget(:).speed]*3.6)
fprintf('identified speed (km/h): \n')
fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[vehicleSpeed]*3.6)


fprintf('Target arrival time (s): \n')
fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[vehicleTarget(:).tStart])
fprintf('identified arrival time (s): \n')
fprintf([repmat('%2.1f \t',1,Nvehicle),' \n'],[tImpact])

Identification of the vehicles' mass

[vehicleMass] = findMass(Bridge,t,rz,posAcc,tImpact,vehicleSpeed,'fcLow',fMin,'fcUp',fMax);
fprintf('Target mass (kg): \n')
fprintf([repmat('%4.0f \t',1,Nvehicle),' \n'],[vehicleTarget(:).mass])
fprintf('identified mass (kg): \n')
fprintf([repmat('%4.0f \t',1,Nvehicle),' \n'],[vehicleMass])

Comparison between the target and simulated bridge response histories

% Create a structure vehicle with the identified parameters of the vehicles.
clear vehicle
Nvehicle = numel(vehicleSpeed);
for ii=1:Nvehicle
    vehicle(ii).mass = abs(vehicleMass(ii));
    vehicle(ii).speed = vehicleSpeed(ii);
    vehicle(ii).direction = 1;
    vehicle(ii).tStart = tImpact(ii);
    vehicle(ii).wn = 0;
end


% COmpute the relative erros on the identified mas, speed and arrival time
errSpeed = nanmean(100.*([vehicle(:).speed]./[vehicleTarget(:).speed]-1));
errMass = nanmean(100.*([vehicle(:).mass]./[vehicleTarget(:).mass]-1));
errTime = nanmean(100.*([vehicle(:).tStart]./[vehicleTarget(:).tStart]-1));
fprintf(['Relative average error on speed: ',num2str(errSpeed,2),' percent \n'])
fprintf(['Relative average error on speed: ',num2str(errMass,2),' percent \n'])
fprintf(['Relative average error on speed: ',num2str(errTime,2),' percent \n'])


% Compute the bridge response with the identified vehicle parameters
[Do2] = dynaResp_vehicle_TD(Bridge,Wind,'vehicle',vehicle); % both wind and traffic
rzFit = squeeze(Do2(2,indY,:));


clf;close all;
figure
plot(t,rz,'k',t,rzFit,'r');
legend('target','fitted')
axis tight
xlabel('Time (s)')
ylabel('r_z (m)')
set(gcf,'color','w')
axis tight


% PSD calculation
[S1,f1] = pmtm(rz,7,numel(rz),fs);
[S2,f2] = pmtm(rzFit,7,numel(rz),fs);
figure
loglog(f1,S1,'k',f2,S2,'r');
legend('target','fitted')
axis tight
ylim([1e-15,1])
xlabel('Frequencies (Hz)')
ylabel('PSD vertical displacement response (m^2/Hz)')
set(gcf,'color','w')

知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1